怎麼樣才能正確的學習生信分析呢？—從學徒做起

  • 2020 年 3 月 10 日
  • 筆記

昨天的我可以為你做些什麼好像閱讀量不高，不過效果還是蠻顯著的，跟部分粉絲聊了聊，希望對他們有幫助吧！ 我們的生信故事會一直在徵稿，詳見；生信故事會欄目徵稿啟事

下面分享一下最新投稿，講述跨專業轉行學習生物資訊學的最正確打開姿勢！

(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，怎麼處理異常數據集，拿到表達矩陣後應該怎麼樣做功能分析，這一步一步走的過程肯定也會遇到不懂不明白的地方，幸運的是，作為學徒有大佬的指點，大佬指一個方向，你是可以節省好多時間。

此時此刻，我想到是杜甫的《春夜喜雨》：好雨知時節，當春乃發生。隨風潛入夜，潤物細無聲。

VirMach 便宜 VPS

QNews

QNews