【原】PCA图显示分组无差异，怎么办？

健明 2021-07-14

展开全文

最近接到粉丝提问，感兴趣的数据集做差异分析，发现很勉强，不好把握。因为我以前在生信技能树写过教程：PCA都分不开的两个组强行找差异是为何，所以征求我的意见。但是我很忙啊，所以就把这个数据集安排给了实习生和学徒。我一直强调，做表达矩阵分析一定要有三张图，见：你确定你的差异基因找对了吗？，基本上看到这3张图，就明白问题出在哪里了。

粉丝求助的数据集介绍如下：

小白菊内酯处理人胆管癌细胞的基因变化
药物处理：小白菊内酯
Platform：GPL6102
Series：GSE22633
https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE22633

第一步是下载数据

rm(list = ls()) ## 魔幻操作，一键清空~

options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型

# 注意查看下载文件的大小，检查数据f='GSE22633_eSet.Rdata'library(GEOquery)# 这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FALSE)if(!file.exists(f)){ gset <- getGEO('GSE22633', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件 save(gset,file=f) ## 保存到本地}load('GSE22633_eSet.Rdata') ## 载入数据class(gset) #查看数据类型length(gset) #class(gset[[1]])gset# assayData: 48701 features, 63 samples

# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]] #dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数dim(dat)#看一下dat这个矩阵的维度#dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列boxplot(dat,las=2)dat=dat[apply(dat,1,sd)>0,]dat[dat<0]=1boxplot(dat,las=2)#经观察，表达量比较小，推测已经log了

#提取临床数据pd=pData(a)#因关注小白菊内酯处理细胞系cell line的变化，所以要去除cell type#且cell line: CK系列没有处理组，所以也需要去掉#删除部分数据pd=pd[grepl('cell line',pd$characteristics_ch1),]pd=pd[-(33:38),]table(pd$characteristics_ch1)dat=dat[,rownames(pd)]

上面的代码重点和难点都是在临床信息整理上面。

第二步是分组&注释

##分组 32个group_list=ifelse(grepl('parthenolide',pd$characteristics_ch1.2),'treat','normal')#设置参考水平，对照在前，处理在后group_list = factor(group_list, levels = c("normal","treat"))

#ID转换 GPL6102#http://www./1399.html#BiocManager::install('illuminaHumanv2.db')library(illuminaHumanv2.db)## 获取表达矩阵的行名，就是探针名称probeset <- rownames(dat)## 使用lookup函数，找到探针在illuminaHumanv2.db中的对应基因名称## 如果分析别的芯片数据，把illuminaHumanv2.db更换即可SYMBOL <- annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")## 转换为向量symbol = as.vector(unlist(SYMBOL))#制作ids转换文件ids = data.frame(probeset,symbol,stringsAsFactors = F)

colnames(ids)=c('probe_id','symbol') ids=ids[ids$symbol != 'NA',]

ids=ids[ids$probe_id %in% rownames(dat),]dat[1:4,1:4] dat=dat[ids$probe_id,]

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行

ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids

ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s

dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat

rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名dat[1:4,1:4] #保留每个基因ID第一次出现的信息

save(pd,dat,group_list,file = 'step1-output.Rdata')

上面代码的重点是芯片探针的注释，因为我们对基因才具备理解力。

第三步是检测数据

就是群主一直强调的做表达矩阵分析一定要有三张图，见：你确定你的差异基因找对了吗？

rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)

pca_plot = function(dddd,ggggg){ library("FactoMineR") library("factoextra") df.pca <- PCA(t(dddd), graph = FALSE) fviz_pca_ind(df.pca, #axes = c(2,3), geom.ind = "point", col.ind = ggggg , addEllipses = TRUE, legend.title = "Groups" )}# 下面的 step1-output.Rdata 文件load(file = 'step1-output.Rdata')# 每次都要检测数据dat[1:4,1:4]table(group_list)#normal treat# 16 16g=factor( group_list )g#对象顺序g=relevel(g,'normal')pca_plot(dat,g)ggsave('all_samples_PCA.png')#热图

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个

library(pheatmap)n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化n[n>2]=2n[n< -2]= -2n[1:4,1:4]#两个分组batch=pd$characteristics_ch1batch=factor(batch)annotation_col=data.frame(group=as.character(group_list), pair = as.character(batch))rownames(annotation_col)=colnames(dat)pheatmap(n,show_colnames =F,show_rownames = F, annotation_col=annotation_col，filename = 'top1000_sd.png')

PS ：我看到实习生还自创了一个函数：pca_plot = function(dddd,ggggg)，看起来是比较有编程天赋的，值得大力培养！

PCA图如下

all_samples_PCA

发现问题：
正常组和处理组间靠得很近，可以看出两者间差别不大。但根据我们常识，处理前后细胞表达量应该会有变化的。

热图如下

top1000_sd

热图也显示细胞系间差异比处理组间要强，所以估计可能是因为细胞系之间的批次效应影响过大，而弱化了处理前后的差异。

这三张图，见：你确定你的差异基因找对了吗？非常重要，提升我们这个数据集的质量！

去除批次效应

定义：不同平台的数据，同个样品不同实验条件，以及同一个样品不同时间的数据等等都会产生一种batch effect 。就像这个数据集，不同细胞系的差异就可能成为干扰基因表达量的因素。
如何检测是否存在批次效应：PCA图或者热图

PCA图：看组间中心点之间的距离，若离得远则说明分组间差异大，否则差异小
热图：每列代表样本，每行代表基因。观察色块间的颜色差别是否明显。

去除方法：使用 limma 的 removeBatchEffect 函数
参考网站：多种批次效应去除的方法比较

rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)load(file = 'step1-output.Rdata')# 每次都要检测数据dat[1:4,1:4]table(group_list)library(limma)g=factor( group_list )gg=relevel(g,'normal')design=model.matrix(~g)#影响表达量的批次batch=pd$characteristics_ch1batch=factor(batch)## 使用 limma 的 removeBatchEffect 函数dat[1:4,1:4]ex_b_limma <- removeBatchEffect(dat, batch = batch, design = design)#ex_b_limma 去除批次效应后的矩阵 dim(ex_b_limma)ex_b_limma[1:4,1:4]##差异分析#去除批次效应前差异分析fit=lmFit(dat,design)fit=eBayes(fit)options(digits = 4)deg1=topTable(fit,coef=2,adjust='BH',number = Inf)save(deg1,file = 'deg1.Rdata')#后fit=lmFit(ex_b_limma,design) fit=eBayes(fit)options(digits = 4)topTable(fit,coef=2,adjust='BH')deg2=topTable(fit,coef=2,adjust='BH',number = Inf)pca_plot(ex_b_limma,g)ggsave('ex_b_limma.png')save(deg2,ex_b_limma,group_list,file = 'deg2.Rdata')

rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)load(file = 'deg1.Rdata')head(deg1)## for heatmapif(T){ load(file = 'step1-output.Rdata') # 每次都要检测数据 dat[1:4,1:4] table(group_list) x=deg1$logFC #deg取logFC这列并将其重新赋值给x names(x)=rownames(deg1) #deg取probe_id这列，并将其作为名字给x

  cg=c(names(head(sort(x),100)),#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg

names(tail(sort(x),100))) library(pheatmap) pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行，所得到的矩阵来画热图

  n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

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) #将ac的行名也就分组信息（是'no TNBC’还是'TNBC’）给到n的列名，即热图中位于上方的分组信息

pheatmap(n,show_colnames =F, show_rownames = F, cluster_cols = T, annotation_col=ac,filename = 'pheatmap_top200_DEG.png') #列名注释信息为ac即分组信息

}rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)load(file = 'deg2.Rdata')head(deg2)## for heatmapif(T){ # 每次都要检测数据 ex_b_limma[1:4,1:4] table(group_list) x=deg2$logFC #deg取logFC这列并将其重新赋值给x names(x)=rownames(deg2) #deg取probe_id这列，并将其作为名字给x

  cg=c(names(head(sort(x),100)),#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg

names(tail(sort(x),100))) library(pheatmap) pheatmap(ex_b_limma[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行，所得到的矩阵来画热图

  n=t(scale(t(ex_b_limma[cg,])))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

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) #将ac的行名也就分组信息（是'no TNBC’还是'TNBC’）给到n的列名，即热图中位于上方的分组信息

pheatmap(n,show_colnames =F, show_rownames = F, cluster_cols = T, annotation_col=ac,filename = 'heatmap_top200_DEG2.png') #列名注释信息为ac即分组信息

}

去除批次效应后的PCA图如下：

ex_b_limma

这个去除批次效应的PCA图。校正之后，可以很明显看出两组的差别，证明去除批次效应是有效的。
校正前后top200_DEG2热图比较，也发现弱化了组内差别，凸显出组间
这样，就可用新的矩阵和差异基因进行下一步分析了

总结

挖掘数据集前，务必做好PCA图与热图的检查，观察组间是否有差异，以此确定分组是否正确
要注意去除批次效应后矩阵表达量是会发生变化的（ex_b_limma），所以在画差异基因热图进行检验的时候，要留意修改矩阵
不太确定利用统计学算法去除批次效应是否能否得到真正的生物学差异，目前还没有一个统一的定论。可能需要进一步了解去除批次效应的内部算法才行。而且，并不是所有的批次效应都是可以去除的，见：并不是所有的批次效应都可以被矫正