怎么样才能正确的学习生信分析呢？—从学徒做起

昨天的我可以为你做些什么好像阅读量不高，不过效果还是蛮显著的，跟部分粉丝聊了聊，希望对他们有帮助吧！

我们的生信故事会一直在征稿，详见；生信故事会栏目征稿启事

下面分享一下最新投稿，讲述跨专业转行学习生物信息学的最正确打开姿势！

(PS : 需要排除部分根本就无心向学的，也许是没有时间，也许是畏难。总之，没有学习条件哪怕是找到我，我也无能为力哈)

前情提要

人生际遇各不同。在机缘巧合的情况下，2019年6月份我从纯实验咖进入到了纯生信的办公室。

为什么会把做生信分析人员所处之地称之为办公室了，那是因为我之前的教育背景，我所认为的实验室，有实验台、细胞间、各种分子生物和细胞生物所需的耗材和仪器，还有实验人员的枪（这里特指移液枪）等，但是生信办公室好像并没有这些，更多的是电脑、计算机和生物理论书籍、白板、文具、各种茶和咖啡。我欣然的接受际遇给我带来的惊喜，从容的进入我的生信生涯:smirk:。

有时间人总是会对周遭的事情有所误解，比如觉得自己收藏了资料，就是学会了资料里面的知识，或者进入的一个圈子，就会认为自己是哪个圈子里面的人了。然而，残酷的现实总会一遍一遍又一遍的告诉世人，这是一种严重的自欺欺人的心理。当然，我也是怀着这种心理很幸运的遇到了生信菜鸟团，从而结识了生信技能树，我还很清楚的记得当时是2019年的8月份。技能树是个为生信小白提供学习的场所，我就是从一片空白在这个场所开始学习。从不会看代码，到学会认识数据结构，在到学会自己敲代码，甚至学习曾老师的单细胞分析流程（第一次看对没有R语言基础的我来说简直就是天书）。就这样吭哧吭哧自己走了大半年。

古语有云：学而不思则罔。学了大半年之后我对我自己做了这样的总结，因为我只会敲代码，而且这些都是流程化的，网上有的，简单来说就是照搬。但是，为什么这样分析？这样分析的意义在哪里呢？该怎样解释这样的分析呢？真是人生三大问呀

在这样的情况下我看到了技能树招学徒的微信短文，上面有句话非常吸引我，好的学习习惯、方法以及名师的指导对个人的成长帮助是毋庸置疑的。

基于以上的缘由，我整理好了我在生信分析上遇到的困难，我的迷茫，以及我的教育背景一起发送给了曾老师，希望能做学徒。

PS: 如果你也有需求，通过学徒培养的方式来学习生物信息学，那么请看好上面的要求哦！找到我招学徒的宣传，里面有我的微信二维码哈，如果你还不知道如何搜索我们生信技能树公众号历史教程，自行点击教程学会在技能树公众号历史教程里面根据关键词查询，基本上初学者问题都有解决方案！

学徒初试

1.练习数据挖掘能力，重复GSE7621、GSE4733、GSE6613三个数据集韦恩图的分析

第一次的作业曾老师就诊断出了我的问题，就是只知道分析的流程，并不知道这其中的缘由，生信的知识并没有系统的建立。看到这样的回复，当时我就激动了，太准确了，一针见血。这更加坚定了我当学徒的初心，我需要知其然并知其所以然。

曾老师很有章法，对症下药简直是一流，这样就有了第二次的学徒作业。

2.绘制GSE2513数据集的火山图及热图
这次的数据集很酷😎，其中大有故事可讲，我已经在接下来的实战演练中详细讲解了。

`实战演练`

曾老师布置的学徒作业，针对GSE2513数据集绘制火山图和热图。按照GEO分析流程绘制完后发送给曾老师，结果得到的回复是全部都是错的，一万点暴击呀:sob:！！！

但是那肿么办了，自己硬着头皮看了自己的代码，还写了自己为什么这么做的原因给曾老师，结果还是没有找到其问题的本质。当然作为学徒有着比较好的优势，就是有大佬在一旁指点。经过曾老师的指点，发现这个数据集是经过zscore了的，导致了我们需要去下载原始的cel文件。那么接下来我要一点点讲解这故事经过。

`错误的表达矩阵`

GEO数据的表达矩阵应该用GEOquery这个函数去下载，代码如下：

rm(list = ls ())
options(stringAsFactors = F)
f='GSE2513_eSet.Rdata'    #先设置一个Rdata赋值给对象f，让f成为一个Rdata对象。
library(GEOquery)
if(!file.exists(f)){
      gset <- getGEO('GSE2513',destdir='.',
                 AnnotGPL = F,    #注释文件
                 getGPL = F,        #平台文件 ）
         save(gset,file=f)   #保存文件到本地          
    }
a=gset[[1]]  ##因为GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list   
dat=exprs(a) ##a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
boxplot(dat,las=2)  ##原始数据做箱式图看数据分布                                         

箱式图如下：

boxplot是用来查看数据表达量的分布，但是我通过boxplot并不能对数据有正确的认识，经曾老师指导，从图中我们可以知道此表达矩阵经过了zscore的标准化，大概是为了让数据正态分布吧。那肿么办了，这下相当于是错误的表达矩阵呀:confused:，是无法正确的进行下面的GEO分析的。可能就是绝处逢生吧，这时曾老师非常合时宜的分享了一微信原创给我了。顿时豁然开朗呀，眼前一片繁花似锦:laughing:，原来原来原来我要重新去下载原始数据，那么我要看看GSE2513是什么芯片平台了，根据进入GEO网页文件一看GSE2513是GPL96的芯片平台，然而GPL96对应的是hgu 133系列的芯片注释包，这时候我们就要去下载CEL文件了，不要问我为什么，我只会告诉你一个网址自己去查看。

以上的这段话用到了两个参考资料：

一个是关于CEL格式数据的解释：http://www.bio-info-trainee.com/1580.html

`下载CEL文件`

这个时候来下载原始数据：

得到GSE2513的芯片的原始数据，cel文件的下载地址：
  • ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE2nnn/GSE2513/suppl/GSE2513_RAW.tar
怎么下载了，使用代码下载，或者直接在网页中下载都是可行的，看个人习惯，自行选择。（如果网速不够，我们生信技能树以前也分享过解决方案哈）

`处理CEL文件`

BiocManager::install('oligo',ask = F,update = F)
library(oligo)
BiocManager::install('pd.hg.u133.plus.2',ask = F,update = F)
library(pd.hg.u133.plus.2)
dir='GSE2513_RAW/'  ##让这个文件赋值给dir这个变量
od=getwd()  ##当前所处路径
setwd(dir)  ##设置为文件路径，为了处理CEL文件
celFiles <- list.celfiles(listGzipped = T) ##为了列出路径下文件名中包含.cel或.CEL的文件，返回9值是一个包含有文件名的list，主要oligo包中的该函数只能列出被压缩的cel文件或是不被压缩的cel文件，不能同时列出，因此需要设定参数listGzipped。
celFiles
affyRaw <- read.celfiles(celFiles) ##读取CEL文件
setwd(od)
eSet <- rma(affyRaw)#rma函数是对表达矩阵进行normalization(mas5也有异曲同工之妙)
eSet#其实严格来说，这个芯片得到表达矩阵后，是需要过滤的，但是我想看看过滤和未过滤的对比，所以这里我就先不过滤了

这样处理后，我就得到了eSet这个对象，这与我一开始用GEOquery包下载是一样的，后续的分析代码不变。

`得到表达矩阵和表型信息`

a=eSet
dat=exprs(a)
dim(dat)
#22283features 12samples
dat[1:4,1:4]#每步都要查看数据结构
par(cex = 0.7)
n.sample=ncol(dat)
if(n.sample>40) par(cex=0.5)
cols <- rainbow(n.sample*1.2)
boxplot(dat,col=cols,main='new expression value',las=2,cex.axis=0.6)
#通过这些boxplot可以         看到各个芯片直接数据是否整齐，看了一下并不是很整齐，所以需要用normalizeBetweenArrays这个函数处理一  下，这步QC也很重要。

我们可以看得出来数据并不整齐，这时就需要一个神奇的函数normalizeBetweenArrays来处理一下，代码如下：

dat=normalizeBetweenArrays(dat)#运用quantile的进行组间归一化，如果上一步boxplot样本间不整齐，这   步肯定要执行的。只能在同一个数据集里面使用。
dat[1:4,1:4]
boxplot(dat,col=cols,main='new expression value',las=2,cex.axis=0.6)#可以看到不整齐的那些已经处理好了

这样数据整齐，接下来就可以使用GEO流程来绘制热图和火山图了。

`PCA主成分分析`

为什么要有PCA呢？如果数据之中的某些维度之间存在较强的线性相关关系，那么样本在这两个维度上提供的信息就会有一定程度上的重复，所以我们希望数据各个维度之间是不相关的（也就是正交的）。此外，出于降低处理数据的计算量或去除噪音等目的，我们也希望能够将数据集中一些不那么重要（方差小）的维度剔除掉。这样将样本从输入空间通过线性或非线性映射到一个低维空间，减少了后续步骤处理的计算量，当降至三维以下时还可用于可视化技术，从而发挥人在低纬空间感知上的优点，发现数据集的空间分析、数据样本之间的差异、聚类性质等结构特征。总之，PCA对于分析样本的相关性具有自己独有的优势。代码如下：

dat[1:4,1:4]    #每步都要检查数据
dat=t(dat)    #转置矩阵，这步非常重要，如果之前的QC没有做好，这步容易出现NA值。转置的目的是符合PCA的输入数据格式：数值型数据，每一列是一个指标(例如基因)，每一行是一个观测(例如样本)
dat=as.data.frame(dat)    #转换成矩阵
dim(dat)    #查看一下维度
library(stringr)
data=data.table::fread('gsm_group.txt',data.table=F,header=F)    #多GSM排序不规则，利用stringr进行分组，做PCA要用到分组信息
names(data)=c('GSM','group')    #给数据的两列命变量名
head(data)
group_list=str_match(data$group,'(.*).*')    #正则表达式的哲学人生，请自行看我给出的网站解释
class(group_list)
head(group_list)
group=data.frame(GSM=data$GSM,group=group_list[,2])
head(group)    
#小白的人生就是每处理完一个数据格式，都要查看一下。
b=group[,2]
class(b)
b=as.character(b)    #一定要变因子
class(b)
colnames(b)='group'
table(b)
b$group[b$group=='Pterygium']='pterygium'    #在原始的txt中本来只有两个处理条件的，但是R是区分大小写的，所以如果你有一个样本是Pterygium，它也会认为这个样本和pterygium是不同的处理条件。
dat=cbind(dat,b)    #加上分组信息
library('FactoMineR')
library('factoextra')    #PCA必要用的两个函数
dim(dat)    #我是小白，每步我都要看看维度，才能放心
dat.pca=PCA(dat[,-ncol(dat)],graph=F)
fviz_pca_ind(dat.pca,
            geom.ind='point',
            col.ind=dat$group,
            palette=c('#00AFBB','#E7B800'),
            addEllipses=T,
            legend.title='Groups_PCA')    #画出PCA图

出图之后，我们来总结一下，用PCA看分组，即检验现在的表达矩阵中的样本信息所对应的分组信息是否有以下情况：

  • 是否有离群样本;
  • 实验组和对照组是否正确（有没有标反）;
  • 有没有批次效应。
    所以跟着PCA分析样本进行聚类，代表样本的点在坐标轴上距离越远，说明样本差异越大。
    上图里面的两个分组，虽然样本有混杂，但是两个组的中心点还是可以分开的，后续分析问题不大。

`差异分析`

两个不同分组之间基因表达值有差异的基因。一般通过两个指标去进行筛选：Fold change（变化倍数，简称FC），以及P value（P值）。常用FC阈值为2或者1，Pvalue的阈值为0.05或者0.01

dat[1:4,1:4]    #我是小白，我检查
library(limma)    #利用limma包进行差异分析
design=model.matrix(~0+factor(group_list))    #这是可以通过limma包做到随心所欲的指定任意两组进行比较的代码
class(design)
dim(design)
colnames(design)=levels(factor(group_list))
exprSet=dat
exprSet[1:4,1:4]
rownames(design)=colnames(exprSet)
design    ##design就是我们制作好的分组矩阵，需要根据我们下载的芯片数据的实验设计方案来，此处例子是GSE2513数据集来探究，12个样本分成了两组，当然读者有自己的数据只需要按照同样方法制作即可！
contrast.matrix=makeContrasts('ptergium-conjunctiva',
                              levels=design）    ##contranst.matrix，这是要把pterygium组合conjunctiva组进行差异分析比较。                              
deg=function(exprSet,design,contrast.matrix){
                                  fit=lmFit(exprSet,design)
                                  fit2=contrasts.fit(fit,contrast.matrix)##这一步很重要
                                  fit2=eBayes(fit2)
                                  tempOutput=topTable(fit2,coef=1,n=Inf)
                                  nrDEG=na.omit(tempOutput)
                                  head(nrDEG)
                                  return(nrDEG)
 }
deg=deg(exprSet,design,contrast.matrix)
dim(deg)
head(deg)
save(deg,file = 'deg.Rdata')
load(file = 'deg.Rdata')
head(deg)
bp(dat[rownames(deg)[1],])
#为deg数据框添加1列
#1.加symbol列，把行名变成一列
library(dplyr)
deg <- mutate(deg,symbol=rownames(deg))
head(deg)
nrDEG=deg
head(nrDEG)
attach(nrDEG)#直接使用数据框中的元素，不用在用$
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
head(df)
df$v=-log10(P.Value)
dim(df)
ggscatter(df,x='logFC',y='v',size=0.5)
df$change=ifelse(df$P.Value > 0.05,'stable',
                 ifelse( df$logFC > 0.5849625,'up',
                        ifelse( df$logFC < -0.5849625,'down','stable') )
               )     ##logFC的条件照着文章处理
head(df)
table(df$change)
write.csv(df,file='DEGs.csv')

在基因表达芯片数据的分析中，差异分析的主角就是limma包，当然也有其它的分析方法，但是limma无疑是其中的佼佼者。

使用这个包需要三个数据，在以上的代码中我都有制作：
  • 表达矩阵
  • 分组矩阵
  • 差异比较矩阵
而且总共也只有三个步骤，在以上的代码中也有体现，现在只是总结
  • lmFit
  • eBayes
  • topTable

`火山图和热图绘制`

差异基因表达出来之后我们应该可视化，通过火山图我们可以方便直观地展示两个样本间基因差异表达的分布情况。
火山图：
  • 通常横坐标用log2(fold change)表示，差异越大的基因分布在两端，
  • 纵坐标用-log10(pvalue)表示，T检验显著性P值的负对数，
  • 通常差异倍数越大的基因T检验越显著，所以往往关注左上角和右上角的值。

ggscatter(df,x='logFC',y='v',size = 0.5,color = 'change')
ggscatter(df,x='logFC',y='v',color = 'change',size = 0.5,
          label = 'symbol',repel = T,
          label.select = c('NDRG1','TNNI2','JUNB','MYC','CYP26A1','ATF3','UPK1B','RAB38','SERPINB1','ATP1B3','RBP1','TYRP1'),
          palette = c('#00AFBB','#D3D3D3','#FC4E07'))
ggsave('volcano.png')
dev.off()

上调和下调基因的热图可视化，让我们看到组与组之间基因的变化趋势

dat[1:4,1:4]
head(df)
table(df$change)
df1=df[df$change=='up'|df$change=='down',]
dim(df1)
head(df1)
gene_up=df1[df1$change=='up',]
length(rownames(gene_up))
gene_down=df1[df1$change=='down',]
length(rownames(gene_down))
x=gene_up$logFC
head(x)
names(x)=rownames(gene_up)
head(x)
y=gene_down$logFC
head(y)
names(y)=rownames(gene_down)
head(y)
cg=c(names(head(sort(x),50)),names(head(sort(y),50)))#上调基因和下调基因各取50画图。
head(cg)
length(cg)
dat[1:4,1:4]
pheatmap(dat[cg,],show_rownames = F,show_colnames = F)
n=t(scale(t(dat[cg,])))
dim(n)
head(n)
max(n)
min(n)
n[n>2]=2
n[n< -2]=-2
n[1:4,1:4]
pheatmap(n,show_rownames = F,show_colnames = F)
ac=data.frame(g=group_list)
ac
rownames(ac)=colnames(n)
ac
pheatmap(n,
         show_colnames = T,
         show_rownames = T,
         cluster_cols = T,
         #cluster_rows = T,
         border_color = NA,
         annotation_col = ac,
         annotation_legend = T,
         cellwidth = 30,
         cellheight = 8,
         fontsize = 8,
         filename = 'heatmap_up50_down50_DEG.png')

以上就是曾老师布置的学徒作业，初试GEO分析中的火山图和热图。看我在这里讲了一大堆，能放在文章里面的也只有火山图和热图了。

可是这个过程走来我是清楚分析流程中的QC，怎么处理异常数据集，拿到表达矩阵后应该怎么样做功能分析，这一步一步走的过程肯定也会遇到不懂不明白的地方，幸运的是，作为学徒有大佬的指点，大佬指一个方向，你是可以节省好多时间。

此时此刻，我想到是杜甫的《春夜喜雨》：好雨知时节，当春乃发生。随风潜入夜，润物细无声。

赞 (0)

相关推荐