PCA图显示分组无差异,怎么办?
时间:2022-07-24
本文章向大家介绍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热图比较,也发现弱化了组内差别,凸显出组间
- 这样,就可用新的矩阵和差异基因进行下一步分析了
总结
- 挖掘数据集前,务必做好PCA图与热图的检查,观察组间是否有差异,以此确定分组是否正确
- 要注意去除批次效应后矩阵表达量是会发生变化的 (ex_b_limma),所以在画差异基因热图进行检验的时候,要留意修改矩阵
- 不太确定利用统计学算法去除批次效应是否能否得到真正的生物学差异,目前还没有一个统一的定论。可能需要进一步了解去除批次效应的内部算法才行。而且,并不是所有的批次效应都是可以去除的,见:并不是所有的批次效应都可以被矫正
- 通过“四大行为”对WCF的扩展[原理篇]
- WCF客户端运行时架构体系详解[下篇]
- WCF客户端运行时架构体系详解[上篇]
- WCF服务端运行时架构体系详解[续篇]
- [WCF-Discovery] 实例演示:如何利用服务发现机制实现服务的“动态”调用?
- [WCF-Discovery]服务如何能被”发现”
- 我的数据访问函数库的源代码(一)—— 共用部分
- 《WCF服务编程》关于“队列服务”一个值得商榷的地方
- 我的数据访问函数库的源代码(二)—— SQL语句部分
- 来源于WCF的设计模式:可扩展对象模式[上篇]
- 我的数据访问函数库的源代码(三)——返回结构数组
- 我的数据访问函数库的源代码(四)—— 存储过程部分,包括存储过程的参数的封装
- [WCF 4.0新特性] 路由服务[实例篇]
- [WCF 4.0新特性] 默认终结点
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法