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

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

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

• 小白菊内酯处理人胆管癌细胞的基因变化

• 药物处理：小白菊内酯

• Platform：GPL6102

• Series：GSE22633

• https://www.ncbi.nlm.nih.gov/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平台，所以下载到的是一个含有一个元素的list
a=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]=1
boxplot(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.bio-info-trainee.com/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     16
g=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]=2
n[n< -2]= -2
n[1:4,1:4]
#两个分组
batch=pd$characteristics_ch1
batch=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 )
g
g=relevel(g,'normal')
design=model.matrix(~g)
#影响表达量的批次
batch=pd$characteristics_ch1
batch=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 heatmap
if(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 heatmap
if(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热图比较，也发现弱化了组内差别，凸显出组间
• 这样，就可用新的矩阵和差异基因进行下一步分析了

总结

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

如果你也需要免费的数据分析

我们推文里面提到的各种各样的数据分析环节都是我非常有经验的，比如我在lncRNA的一些基础知识 ，和lncRNA芯片的一般分析流程 介绍过的那些图表，以及下面的目录的分析内容 对我来说是举手之劳，希望可以帮助到你！

• 转录组数据分析的4个维度认识（数据分析继续免费哦） RNA-seq数据的2个分组差异分析，热图，PCA图，火山图等等
• 根据感兴趣基因看肝癌免疫微环境的T细胞亚群差异 条形图或者箱线图
• 查看感兴趣基因的甲基化水平和RNA表达水平(数据分析免费做)相关性 散点图或者箱线图

你只需要发送数据分析要求，以及简短的项目描述到我的邮箱 [e-mail]

邮件正文最好是加上你是啥时候认识生信技能树的哦，或者其它一些寒暄的话，自我介绍也行。主要是考虑到可能想免费分析数据的朋友很多，所以会根据你的来信，*我主观判定一个优先级哦。目前我有20多个愿意长期在我的指导下进行数据探索的学徒，等我的团队扩大到200人，我们应该是可以做到数据分析全部免费，敬请期待哈！

如果你想亲自分析自己的数据

那就需要考虑花一个月到半年时间学习了，我们力推生信技能树官方举办的学习班：

• 数据挖掘学习班第6期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。
• 生信爆款入门-第8期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。