最近接到粉丝提问,感兴趣的数据集做差异分析,发现很勉强,不好把握。因为我以前在生信技能树写过教程:PCA都分不开的两个组强行找差异是为何,所以征求我的意见。但是我很忙啊,所以就把这个数据集安排给了实习生和学徒。我一直强调,做表达矩阵分析一定要有三张图,见:你确定你的差异基因找对了吗? ,基本上看到这3张图,就明白问题出在哪里了。
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),看起来是比较有编程天赋的,值得大力培养!
all_samples_PCA
top1000_sd
这三张图,见:你确定你的差异基因找对了吗? 非常重要,提升我们这个数据集的质量!
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