AGS测序下游分析一条龙

1. 表达矩阵

ensembl_matrix.Rdata,原始矩阵

rm(list=ls())
options(stringsAsFactors = F)
library(stringr)  

  a1=read.table('/mnt/AGS_RNA-seq/featureCounts/all.counts.txt',header = T)
  dim(a1)
  a1[1:4,1:4]
  mat= a1[,7:ncol(a1)] 
  rownames(mat)=a1$Geneid
  mat[1:4,1:4]
  keep_feature <- rowSums (mat > 1) > 1
  table(keep_feature)
  mat <- mat[keep_feature, ]
  mat[1:4,1:4]
  dim(mat) 
  colnames(mat)
  #删去62，194组内相关性较差的两组
  #mat <- mat[,-2]
  #mat <- mat[,-5]
  colnames(mat)=c("61","62","63","191","193","194")
  colnames(mat)
  ensembl_matrix=mat
  colnames(ensembl_matrix)
  
  #添加分组信息
  b=read.csv('~/AGS_RNAseq/AGS_4/logs/samples.txt')
  group_list=b$status
  table(group_list)
  save(ensembl_matrix,group_list,file='ensembl_matrix.Rdata')

#samples.txt
Sample,status
61,high
62,high
63,high
191,low
193,low
194,low

2. 数据check，PCA，热图，相关性分析

rm(list = ls()) 
options(stringsAsFactors = F)
load(file = 'ensembl_matrix.Rdata')

# 每次都要检测数据
exprSet=ensembl_matrix
dat=log2(edgeR::cpm(exprSet)+1)
dat[1:4,1:4] 
colnames(exprSet)
table(group_list)

## PCA
dat=t(dat) 
library("FactoMineR")
library("factoextra")  
dat.pca <- PCA(dat , graph = FALSE)#现在dat最后列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind =  group_list, # color by groups
             #palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
ggsave('all_samples_PCA_by_type.png')
 

exprSet=ensembl_matrix
dat=log2(edgeR::cpm(exprSet)+1)
dat[1:4,1:4] 
table(group_list)

#表达数据log后较之前均一，看相关性
View(dat)
colnames(dat)
[1] "61"  "62"  "63"  "191" "193" "194"
pheatmap::pheatmap(cor(dat))
colD=data.frame(group=group_list)
rownames(colD)=colnames(dat)
pheatmap::pheatmap(cor(dat),annotation_col = colD, 
                   show_rownames = T, 
                   color = colorRampPalette(brewer.pal(9, "OrRd"))(50),
                   cluster_rows = FALSE, 
                   cluster_cols = FALSE, 
                   display_numbers = T, 
                   number_format = "%.3f", ##显示相关系数，三位小数
                   filename = 'cor_log2_all.png')


##热图

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n) 
pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=ac)
pheatmap(n,show_colnames =T,show_rownames = F,
         annotation_col=ac,filename = 'heatmap_top1000_sd.png')


# 相似性，原始矩阵
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'ensembl_matrix.Rdata')

exprSet=ensembl_matrix
colnames(exprSet)
colD=data.frame(group=group_list)
rownames(colD)=colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),annotation_col = colD, 
                   show_rownames = T, 
                   color = colorRampPalette(brewer.pal(9, "OrRd"))(50),
                   cluster_rows = FALSE, 
                   cluster_cols = FALSE, 
                   display_numbers = T, 
                   number_f