不一定正确的多分组差异分析结果热图展现
最近给学徒安排了一个文献数据挖掘任务，但是检查她结果发现跟文献差异有点大，因为时间关系，来不及全盘细致检查，所以先发出来交给读者考量。主要是基于我前些天在生信技能树分享的一个策略：如果你的分组比较多，差异分析策略有哪些？
1.复现的图片
Paper:Machine learning analysis of gene expression data reveals novel diagnostic and prognostic biomarkers and identifies therapeutic targets for soft tissue sarcomas
2.术语：
3.Something to clarify
不是直接从TCGA下载的，而是文章整理的，文章是下面这篇：
Comprehensive and Integrated Genomic Characterization of Adult Soft Tissue Sarcomas
这里的热图不是用表达矩阵作图，而是，用的logFC；
低表达基因作过滤用的是，cpm<2，差异分析用的是logFC>0和Adjust.pvalue <0.05；
GO富集分析用的是enrichR；
Step01-Data download
GDCquery+GDCdownload+GDCprepare进行数据下载, 大家可能需要从头理解TCGA背景知识：悄咪咪的上线了TCGA知识图谱视频教程（B站和YouTube直达）
rm(list=ls())
library(TCGAbiolinks)
library(SummarizedExperiment)
query <- GDCquery(project = "TCGA-SARC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
GDCdownload(query)
SARC_DATA<-GDCprepare(query,save=FALSE)
save(SARC_DATA,file='sarc.Rdata')
Step02-Data_for_analyse
上文有写，临床信息来自一篇文章，NIHMS912614-supplement-10.xlsx为其附表；
rm(list=ls())
load('sarc.Rdata')
library(SummarizedExperiment)
summary(SARC_DATA)
sarc_data<- assay(SARC_DATA)
dim(sarc_data)
####临床信息来源于另一篇paper的附表
library('openxlsx')
clinic <- read.xlsx('NIHMS912614-supplement-10.xlsx',sheet = 2,startRow = 2)
colnames(sarc_data)<- substr(colnames(sarc_data),start = 1,stop = 15)
sarc_choose<- sarc_data[,clinic$TCGA.barcode]
save(sarc_choose,clinic,file='ex_pd.Rdata')
Step03-Differential analysis
阈值设定（cpm）均源于我对文章的理解，差异分析的步骤也是如此；
这里整个差异分析过程，采用了循环的方式进行，里面主要在循环的就是group_list和design，可能需要详细看才能理解；
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'ex_pd.Rdata')
library(edgeR)
#####过滤掉cpm<2的counts
tmp<- cpm(sarc_choose)
tmp1<- tmp[rowSums(tmp<2)>1,]
sarc_choose<- sarc_choose[rownames(tmp1),]
identical(colnames(sarc_choose),clinic$TCGA.barcode)
####用TMM方法进行normalization
dge <- DGEList(counts=sarc_choose)
dge <- calcNormFactors(dge,method = 'TMM')
####构建对应的group_list
subtype <- unique(clinic$short.histo)
group_list <- c()
for( i in 1:length(subtype)){
tmp<- ifelse(clinic$short.histo==subtype[i],subtype[i],'other')
group_list<- cbind(group_list,tmp)
colnames(group_list)[i] <- subtype[i]
}
#####差异分析和PCA
deg_re <- matrix()
for( i in 1:ncol(group_list)){
design <- model.matrix(~0+factor(group_list[,i]))
colnames(design)=levels(factor(group_list[,i]))
rownames(design)=colnames(sarc_choose)
data_v <- voom(dge,design)
###PCA
library("FactoMineR")
library("factoextra")
####
df.pca <- PCA(t(data_v$E), graph = FALSE)
fviz_pca_ind(df.pca,
geom.ind = "point",
col.ind = group_list[,i],
addEllipses = TRUE,
legend.title = "Groups"
)
title<- paste0('./figures/',colnames(group_list)[i],'.png')
ggsave(title)
constrsts<- paste0(colnames(group_list)[i],"-other")
contrast.matrix<-makeContrasts(contrasts=constrsts,levels = design)
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
fit <- lmFit(data_v,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1, 0.01)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
tT_tmp <- topTable(fit2, adjust="fdr", number=nrow(fit2))
tT_tmp<- subset(tT_tmp, select=c('adj.P.Val',"P.Value","logFC"))
colnames(tT_tmp)<- paste0(colnames(group_list)[i],colnames(tT_tmp))
deg_re <- cbind(deg_re,tT_tmp)
}
deg_re <- deg_re[,-1]
save(deg_re,file='deg.Rdata')
Step04-pheatmap
文章是用logFC做热图，这里等于挑选了，差异基因并集的logFC；
之前直接用logFC>0的阈值进行热图绘制，热图很丑，后面把阈值调高了进行绘图（曾老师原话是排序，这里并没有这样）；
看文章的图，像是聚类后的，但又没有展示tree，我就把聚类结果存为变量，再进行的热图绘制；
rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
histo<- unique(clinic$short.histo)
keep <- c()
for(i in 1:length(histo)){
tT<- deg_re[,grepl(histo[i],colnames(deg_re))]
colnames(tT)<- gsub(histo[i],'',colnames(tT))
####这个logFC的阈值是文章选的
tT$change = ifelse(tT$adj.P.Val < 0.05 & tT$logFC>2,TRUE,FALSE)
keep <- cbind(keep,tT$change)
}
#####文章是用logFC做热图，这里等于挑选了，差异基因并集的logFC
logFC_pheatmap<- deg_re[apply(keep,1,function(x){sum(x)>0}),grepl('logFC',colnames(deg_re))]
colnames(logFC_pheatmap) <- gsub('logFC','',colnames(logFC_pheatmap))
nrow(logFC_pheatmap)
library(pheatmap)
tmp<- pheatmap(logFC_pheatmap,silent=T)
####调出热图聚类结果，并存储变量
final_phe<- logFC_pheatmap[tmp$tree_row$order,tmp$tree_col$order]
bk = c(seq(min(final_phe),max(final_phe), length=100))
pheatmap(final_phe,
breaks=bk,
show_rownames = F,
cluster_rows=F,
cluster_cols = F,
color = c(colorRampPalette(c("blue", 'white','red'))(100)))
Step05-GO enrichment
paper中富集所用的database为GO_Biological_Process_2015和R包为enrichR；
最后的富集结果，某些相应的Adjust.pvalue 没有满足条件的，我选了前三条，共7+7*3行，这里截图展示；
rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
histo<- unique(clinic$short.histo)
library(org.Hs.eg.db)
E2S<-select(org.Hs.eg.db,keys=keys(org.Hs.eg.db),columns = c("SYMBOL",'ENSEMBL'))
deg_re$ENSEMBL <- rownames(deg_re)
deg_re<- merge(deg_re,E2S,by='ENSEMBL')
deg_re <- na.omit(deg_re)
GO_histo <- list()
for(i in 1:length(histo)){
tT<- deg_re[,grepl(histo[i],colnames(deg_re))]
colnames(tT)<- gsub(histo[i],'',colnames(tT))
tT$change = ifelse(tT$adj.P.Val < 0.05 & tT$logFC0,TRUE,FALSE)
######paper中富集所用的database为GO_Biological_Process_2015和R包为enrichR
library(enrichR)
dbs <- c("GO_Biological_Process_2015")
enriched <- enrichr(deg_re$SYMBOL[tT$change], dbs)
bp <- enriched[["GO_Biological_Process_2015"]]
GO_histo[[histo[i]]]<- bp[1:3,]
}
GO_final <- c()
for(i in 1:length(GO_histo)){
name <- paste0(names(GO_histo)[i])
tmp<- GO_histo[[i]]
tmp<- c(name,as.character(tmp[,grepl('Term',colnames(tmp))]))
GO_final <- c(GO_final,tmp)
}
GO_final <- as.data.frame(GO_final)
因为并没有看到GO的数据库注释结果，而且在我检查热图里面各个癌症亚型的差异基因数量的时候，发现跟原文不一致！
首先，可以看到原文数量如下:
然后我使用代码检查学徒的差异分析结果：
rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
sarc_choose=sarc_choose[rownames(sarc_choose) %in% rownames(deg_re),]
histo<- unique(clinic$short.histo)
keep=do.call(cbind,lapply(histo, function(x){
print(x)
cc=deg_re[,grep(x,colnames(deg_re))]
return(cc[,1]<0.05 & cc[,3]>2)
}))
colnames(keep)=histo
rownames(keep)=rownames(deg_re)
apply(keep,2,table)
得到的差异基因数量分布如下：
原文是DDLPS的差异基因数量最多，而学徒代码里面是 SS 的最多，所以中间过程还是需要仔细检查，争取找到差异。
不过，今天是国庆节最后一天，我还在路上，就不耗费时间在这上面了，感兴趣的粉丝可以根据推文里面的代码来尝试学习+理解+发给我你的见解！
我还想说一下的是，其实这个项目的样本数量已经超过200，而且分组是7个，已经可以使用单细胞数据分析策略了，比如我这里就使用seurat对表达矩阵进行tSNE可视化：
还可以使用分组，看看与文章的相似性
文章分组如下：
还有更多单细胞探索，大家自行学习，比如 DoHeatmap 可以绘制我给学徒的任务。