用R分析晶片實戰（上）

劉小澤寫於18.9.10
實戰上從資料下載到差異基因的獲得、初步作圖
實戰下進行富集分析，使用資料庫進行註釋
在此感謝jimmy的limma包等教程

1 下載晶片資料

使用Fabbrini E於2015年發表在J Clin Invest上的文章：Metabolically normal obese people are protected from adverse effects following weight gain

文章研究了體重適度增加對代謝正常(MNO)和代謝異常(MAO)受試者脂肪組織基因表達的影響，選取11個代謝正常的人、7個代謝異常的人增重前後的皮下脂肪組織，得到36個資料

文章下載連結：https://pdfs.semanticscholar.org/6052/f035dfb296559e19314352506f048e13f02e.pdf

使用的資料集是GSE62832

實驗平臺是GPL6244， Affymetrix Human Gene 1.0 ST Array

• 方法一：開啟網頁https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62832 【可能會有網速限制】
  
  

  GEO資料GSE62832

• 方法二：使用GEOquery包

  if(T){
  source("http://bioconductor.org/biocLite.R")
  options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
  biocLite("GEOquery")
  
  library(GEOquery)
  eSet <- getGEO('GSE62832', destdir = '.', getGPL = F, AnnotGPL = F)#一般晶片的註釋檔案比較大，直接下載不容易成功，後兩個選項可以避免下載GPL、Annotation
  save(eSet, file = 'GSE62832.eSet.Rdata')
  }
  

下載完，看一下eSet這個物件都包含什麼

> eSet
$GSE62832_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 36 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1534034 GSM1534035 ... GSM1534069 (36 total) #36個樣本
  varLabels: title geo_accession ... tissue:ch1 (35 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: GPL6244 #檢視對應晶片平臺

class(eSet) #先檢視資料種類【列表還是資料框】
str(eSet) #再看資料結構【可以看到第二行涉及到了表達矩陣】
#既然是列表，那麼從列表中提取資訊就是 eSet[[1]]
e <- exprs(eSet[[1]])

2. 將表達矩陣的探針ID轉換為gene ID

• 首先需要知道GSE62832對應平臺是GPL6244

• 然後需要知道GPL6244對應哪個註釋包【閱讀jimmy之前整理的平臺資訊https://www.jianshu.com/p/f6906ba703a0】

• 下載相應的註釋包

  source("https://bioconductor.org/biocLite.R")
  options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
  biocLite("hugene10sttranscriptcluster.db")
  

• 開始探索、過濾、整合

  library(hugene10sttranscriptcluster.db) 
  ls("package:hugene10sttranscriptcluster.db") #檢視所有包含的物件
  s <- toTable(hugene10sttranscriptclusterSYMBOL) #將SYMBOL物件轉為資料框
  #【補充：這個註釋包是別人上傳的，我們只是拿來用。其中的表達矩陣中原來有3w多探針，作者過濾後只留下接近2w的探針。這僅僅是作者過濾後的結果，並非原始資料結果】
  

  簡單的探索

  library(magrittr) 
  s$symbol %>% unique() %>% length() #看一下有多少基因【人類正常蛋白編碼基因就2w左右】
  s$symbol %>% table() %>% sort() %>% tail() #基因使用探針最多的前6名【發現有兩個基因用了10個探針】
  s$symbol %>% table() %>% sort() %>% table()
  

  過濾+整合

  dim(e) #過濾前
  #過濾【只保留和註釋檔案探針id相同的探針】
  efilt <- e[rownames(e)%in%s$probe_id,]
  dim(efilt)#過濾後
  #整合1【目的：保證一個基因對應一個探針；如果基因和探針一一對應很好說，但如果一個基因對應多個探針：每個探針取一行的均值-》對應同一基因的探針取表達量最大的探針-》按照基因名給他們建索引，因為是按照基因來過濾探針（不用s$probe_id構建索引的原因是，看清楚我們的目的是讓註釋包的一個基因對應我們自己表達矩陣的一個探針。如果用s$probe_id那麼結果就成了讓註釋包的一個探針對應我們自己表達矩陣的一個探針，當然這樣也執行不成功，因為自己表達矩陣的探針過濾後的數量和註釋包的探針數量不相等，這樣沒法一一對應。但基因名數量是不變的，什麼是索引？以不變應萬變的就是索引）】
  maxp = by(efilt,s$symbol,function(x) rownames(x)[which.max(rowMeans(x))]) 
  uniprobes = as.character(maxp)
  efilt=efilt[rownames(efilt)%in%uniprobes,]
  #整合2【目的：將我們表達矩陣的行名換成剛才一對一的基因名，並且match這個函式保證了表達矩陣和註釋包的順序是一致的】
  rownames(efilt)=s[match(rownames(efilt),s$probe_id),2]
  

3. 對錶達矩陣進行一些檢驗

表達矩陣不侷限於GEO資料庫的晶片分析，轉錄組及其他涉及基因、樣本、分組關係的都會有一個表達量矩陣，就是一個基因在不同樣本中（對照、處理；是否患病等）的表達差異。
拿到表達矩陣，在進行後續分析之前，首先要檢測這個矩陣是不是合理的，比如看管家基因是否表達量突出、一致；樣本分組是不是和實驗設計一致，用PCA、hclust檢驗

3.1 檢測一些管家基因表達量

boxplot(efilt[,1])#看看第一個樣本中總體基因表達量分佈，可以看到基本為5左右
efilt['ACTB',] #激動蛋白Beta-actin的基因名是ACTB，管家基因
efilt['GAPDH',] #也是管家基因

3.2 看錶達矩陣的整體分佈

先把表達矩陣=》tidy data【四列：基因名、樣本、表達量、表型分組(看文獻按MAO、MNO分組)】

library(reshape2)
pdata=pData(eSet[[1]]) #將樣本表型資訊從資料框中提取出來【取出來的是表型、樣本的資料框】
group_list=as.character(pdata$`metabolic status:ch1`) 
m_efilt = melt(efilt) #先將原來矩陣“融化
colnames(m_efilt)=c('symbol','sample','value') #重新命名三列
m_efilt$group=rep(group_list,each=nrow(efilt)) 

以圖為證

4. 對錶達矩陣進行差異分析

只需要提供表達矩陣efilt、分組資訊group_list，就能使用limma進行分析
詳細內容參考jimmy的http://www.bio-info-trainee.com/bioconductor_China/software/limma.html

suppressMessages(library(limma))
#limma需要三個矩陣：表達矩陣（efilt）、分組矩陣（design）、比較矩陣（contrast）
#先做一個分組矩陣～design，說明MAO是哪幾個樣本，MNO又是哪幾個，其中1代表“是”
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(efilt)
design
#再做一個比較矩陣【一般是case比control】
contrast<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast

4.1 準備就緒，就可以開始差異分析

DEG <- function(efilt,design,contrast){
  ##step1
  fit <- lmFit(efilt,design)
  ##step2
  fit2 <- contrasts.fit(fit, contrast) 
  fit2 <- eBayes(fit2)  
  ##step3
  mtx = topTable(fit2, coef=1, n=Inf)
  deg_mtx = na.omit(mtx) 
  return(deg_mtx)
}
DEG_mtx <- DEG(efilt,design,contrast) #得到全部的差異基因矩陣

4.2 對一個小表達矩陣，如30、50個基因，可以用熱圖

top30_gene=head(rownames(DEG_mtx),30)
top30_matrix=efilt[top30_gene,] #得到top30的表達量矩陣
top30_matrix=t(scale(t(top30_matrix)))
#這裡做個top30 gene heatmap

4.3 對一個大的表達矩陣，如全部的差異基因，可以用火山圖

火山圖實際上就是根據兩列進行作圖：logFC、pvalue

plot(DEG_mtx$logFC, -log10(DEG_mtx$P.Value)) #先畫一個最簡單的圖（能說明原理）
#再根據jimmy的程式碼做一個。網路上也有很多火山圖的程式碼可以參考

歡迎關注我們的公眾號～_～　　
我們是兩個農轉生信的小碩，打造生信星球，想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請後臺留言或傳送郵件到Bioplanet520@outlook.com