Bioconductor 分析基因晶片資料(轉帖）

原文地址：

（程式碼太多，請看原貼）

還有 系統的使用affy處理晶片，看原貼：

1. 快速入門

安裝並載入所需R包

source("");

biocLite(“CLL”)

# 載入CLL包，CLL包會自動呼叫affy包，該包含有一系列處理函式

library(CLL)# read example dataset，（CLL包附帶的示例資料集）data("CLLbatch")# pre-process using RMA methodCLLrma<- rma(CLLbatch)# read expression value after pre-processinge <- exprs(CLLrma)# 檢視部分資料e[1:5,1:5]

image.png

2. 基因晶片資料預處理

2.1 資料輸入

# 載入CLL包，CLL包會自動呼叫affy包，該包含有一系列處理函式library(CLL)# read example dataset，（CLL包附帶的示例資料集）data("CLLbatch")# 檢視資料型別data.class(CLLbatch)# 讀入所有樣本的狀態資訊data(disease)# 檢視所有樣本的狀態資訊disease

image.png

# 檢視"AffyBatch"的詳細介紹help(AffyBatch)phenoData(CLLbatch)featureData(CLLbatch)protocolData(CLLbatch)annotation(CLLbatch)exprs_matrix <- assayData(CLLbatch)[[1]]exprs_matrix[1:5,1:5]exprs(CLLbatch)[1:5,1:5]

image.png

2.2 質量控制

質量控制分三步，直觀觀察，平均值方法，資料擬合方法。這三個層次的質量控制分別由image 函式、simpleaffy 包和affyPLM包實現

2.2.1 用image包對晶片資料進行質量評估

# 檢視第一張晶片的灰度影像image(CLLbatch[,1])

如果影像特別黑，說明型號強度低；如果影像特別亮，說明訊號強度很有可能過飽和

image.png

2.2.2 用simpleaffy包對晶片資料進行質量評估

# 安裝simpleaffy包source('')biocLite("simpleaffy")library(BiocInstaller)biocLite("simpleaffy")library(simpleaffy)library(CLL)data(CLLbatch)# 獲取質量分析報告Data.qc <- qc(CLLbatch)# 圖型化顯示報告plot(Data.qc)

image.png

第一列是所有樣品的名稱

第二列檢出率和平均背景噪聲（往往較高的平均背景噪聲都伴隨著較低的檢出率）

第三列

藍色實現為尺度因子，取值（-3，3）

"o" 不能超過1.25，否則資料質量不好

“三角型“不能超過3，否則資料質量不好

bioB 說明晶片檢測沒有達標

2.2.3 用affyPLM包對晶片資料進行質量評估

source('')biocLite('affyPLM')library(affyPLM)data(CLLbatch)# 對資料集做迴歸分析，結果為一個PLMset型別的物件Pset <- fitPLM(CLLbatch)image(CLLbatch[,1])# 根據計算結果，畫權重圖image(Pset,type="weights",which=1, main="Weights")# 根據計算結果，畫殘差圖image(Pset,type="resids",which=1, main="Residuals")# 根據計算結果，畫殘差符號圖image(Pset,type="sign.resids",which=1, main="Residuals.sign")

image.png

image.png

2.2.4 檢視總體質量

一個樣品的所有探針組的RLE分佈可以用一個統計學中常用的箱型圖形表示

source('')biocLite("RColorBrewer")library(affyPLM)library(RColorBrewer)library(CLL)# 讀入資料data("CLLbatch")# 對資料集做迴歸計算Psel<-fitPLM(CLLbatch)# 載入一組顏色colors<-brewer.pal(12,"Set3")# 繪製RLE箱線圖Mbox(Pset, ylim=c(-1,1),col=colors,main="RLE",las=3);# 繪製NUSE箱線圖boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3)

從以下兩幅圖中可以看出CLL1，CLL10 的質量明顯有別與其他樣品，需要去除

image.png

image.png

2.2.5 檢視RNA降解曲線

source('')biocLite("affy")library(affy)library(RColorBrewer)library(CLL)data("CLLbatch")# 獲取降解資料data.deg <- AffyRNAdeg(CLLbatch)# 載入一組顏色colors <- brewer.pal(12,"Set3")# 繪製RNA降解圖plotAffyRNAdeg(data.deg, col=colors)# 在左上部位加註圖注legend("topleft", rownames(pData(CLLbatch)), col=colors, lwd=1, inset=0.05, cex=0.5)

image.png

去掉三個質量差的

CLLbatch<-CLLbatch[, -match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),                              sampleNames(CLLbatch))]

2.2.6 從聚類分析的角度看資料質量

source('')biocLite("gcrma")biocLite("graph")biocLite("affycoretools")library(CLL)library(gcrma)library(graph)library(affycoretools)data(CLLbatch)data("disease")# 使用gcrma演算法來預處理資料CLLgcrma<-gcrma(CLLbatch)# 提取基因表達矩陣eset<-exprs(CLLgcrma)# 計算樣品兩兩之間的Pearson相關係數pearson_cor<-cor(eset)# 得到Pearson距離的下三角矩陣dist.lower<-as.dist(1-pearson_cor)# 聚類分析hc<-hclust(dist.lower,"ave")plot(hc)# PCA分析samplenames<-sub(pattern="\.CEL", replacement ="",colnames(eset))groups<-factor(disease[,2])plotPCA(eset,addtext=samplenames,groups=groups,groupnames=levels(groups))

從聚類分析的結果來看，穩定組（紅框）和惡化組分不開

image.png

從主成成份分析來看，也分不開

image.png

2.3 背景校正、標準化和彙總

晶片資料透過質量控制，剔除不合格的樣品，留下的樣品資料往往要透過三步處理（背景校正、標準化和彙總）才能得到下一步分析所需要的基因表達矩陣

去除背景噪聲的過程叫背景校正

標準化目的是使各次/組測量或各種實驗條件下的測量可以相互比較

使用一定的統計方法將前面得到的熒光強度值從探針水平彙總到探針組水平

> bgcorrect.methods()[1]"bg.correct""mas""none""rma"> normalize.methods(CLLbatch) [1]"constant""contrasts""invariantset""loess""methods""qspline"[7]"quantiles""quantiles.robust""quantiles.probeset""scaling"> pmcorrect.methods()[1]"mas""methods""pmonly""subtractmm"> express.summary.stat.methods()[1]"avgdiff""liwong""mas""medianpolish""playerout"

引數說明

afbatch輸入資料的型別必須是AffyBatch物件

bgcorrect.method背景校正方法

bgcorrect.param指定背景校正方法的引數

normalize.method標準化方法

normalize.param指定標準化方法的引數

pmcorrect.methodPM調整方法

pmcorrect.param指定PM方法引數

summary.method彙總方法

summary.param指定彙總方法的引數

晶片內標準化（Loess）前後MA圖

# 使用mas方法做背景校正>CLLmas5<- bg.correct(CLLbatch, method ="mas")# 使用constant方法標準化> data_mas5 <- normalize(CLLmas5, method="constant")#　檢視每個樣品的縮放係數> head(pm(data_mas5)/pm(CLLmas5),5)CLL11.CELCLL12.CELCLL14.CELCLL15.CELCLL16.CELCLL17.CELCLL18.CELCLL19.CELCLL20.CELCLL21.CELCLL22.CELCLL23.CEL17521811.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236535668911.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236522769611.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236523791911.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.43236527517311.1558491.0238731.4931931.5493692.0002991.4515761.7765010.98251080.70708280.99587331.432365CLL24.CELCLL2.CELCLL3.CELCLL4.CELCLL5.CELCLL6.CELCLL7.CELCLL8.CELCLL9.CEL1752181.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392483566891.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482276961.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482379191.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.89392482751731.7060261.1563780.84254190.97750820.98165731.1829631.1149761.136390.8939248# 檢視第二個樣品的縮放係數是怎麼計算來的> mean(intensity(CLLmas5)[,1])/mean(intensity(CLLmas5)[,2])[1]1.155849

2.4 預處理的一體化演算法

預處理方法

方法背景校正方法標準化方法彙總方法

MASSmasconstantmas

dChipmasinvariantsetliwong

RMArmaquantilemedianpolish

MAS5 每個晶片可以單獨進行標準化；RMA 採用多晶片模型，需要對所有晶片一起進行標準化。

MAS5 利用MM探針的資訊去除背景噪聲，基本思路是MP-MM；RMA 不使用MM資訊，而是基於PM訊號分佈採用的隨機模型來估計表達值。

RMA處理後的資料是經過2為底的對數轉換的，而MAS5不是，這一點非常重要，因為大多數晶片分析軟體和函式的輸入資料必須經過對數變換。

比較不同演算法的處理效果

library(affy)library(gcrma)library(affyPLM)library(RColorBrewer)library(CLL)data("CLLbatch")colors <- brewer.pal(12,"Set3")# use MAS5CLLmas5<- mas5(CLLbatch)# use rma CLLrma<- rma(CLLbatch)# use gcrmaCLLgcrma<- gcrma(CLLbatch)## hist plothist(CLLbatch, main="orignal", col=colors)legend("topright", rownames(pData(CLLbatch)), col=colors,       lwd=1, inset=0.05, cex=0.5, ncol=3)hist(CLLmas5, main="MAS 5.0", xlim=c(-150,2^10), col=colors)hist(CLLrma, main="RMA", col=colors)hist(CLLgcrma, main="gcRMA", col=colors)## boxplotboxplot(CLLbatch, col=colors, las=3, main="original")boxplot(CLLmas5, col=colors, las=3, ylim=c(0,1000), main="MAS 5.0")boxplot(CLLrma, col=colors, las=3, main="RMA")boxplot(CLLgcrma, col=colors, las=3, main="gcRMA")

RMA演算法將多條曲線重合到了一起，有利於進一步的差異分析，但卻出現了雙峰現象，不符合高斯正態分佈。很顯然gcRMA演算法在這裡表現的更好。

image.png

image.png

image.png

image.png

三個演算法處理後的各樣品的中值都十分接近。MAS5演算法總體而言還不錯，有一定拖尾現象；而gcRMA的拖尾現象比RMA要明顯得多。這說明針對低表達量的基因，RMA演算法比gcRMA演算法表現更好一些。

image.png

image.png

image.png

image.png

透過MA圖來檢視標準化處理的效

library(gcrma)library(RColorBrewer)library(CLL)library(affy)data("CLLbatch")colors<-brewer.pal(12."Set3")CLLgcrma<-gcrma(CLLbatch)MAplot(CLLbatch[,1:4],pairs=TRUE,plot.method="smoothScatter",cex=0.8,main="original MA plot")MAplot(CLLgcrma[,1:4],pairs=TRUE,plot.method="smoothScatter",cex=0.8,main="gcRMA MA plot")

原始資料中，中值（紅線）偏離0，經過gcRMA處理後，中值基本保持在零線上。

image.png

image.png

3 基因晶片資料分析

3.1 選取差異表達基因

######## DEG analysis#　如果沒有安裝limma包，請取消下面兩行註釋# library(BiocInstaller)# biocLite("limma")# 匯入包library(limma)library(gcrma)library(CLL)# 匯入CLL資料data("CLLbatch")data(disease)# 移除 CLL1, CLL10, CLL13CLLbatch<-CLLbatch[, -match(c("CLL10.CEL","CLL1.CEL","CLL13.CEL"),                              sampleNames(CLLbatch))]# 用gcRMA演算法進行預處理CLLgcrma<- gcrma(CLLbatch)# remove .CEL in sample namessampleNames(CLLgcrma) <- gsub(".CEL$","", sampleNames(CLLgcrma))# remove record in data disease about CLL1, 10 and 13.CELdisease <- disease[match(sampleNames(CLLgcrma), disease[,"SampleID"]), ]# 獲得餘下21個樣品的基因表達矩陣eset <- exprs(CLLgcrma)# 提取實驗條件資訊disease <- factor(disease[,"Disease"])# 構建實驗設計矩陣design <- model.matrix(~-1+disease)# 構建對比模型，比較兩個實驗條件下表達資料#　這裡的.是簡寫而不是運算子號contrast.matrix <- makeContrasts(contrasts ="diseaseprogres. - diseasestable",                                 levels=design)# 線性模型擬合fit <- lmFit(eset, design)# 根據對比模型進行差值計算 fit1 <- contrasts.fit(fit, contrast.matrix)# 貝葉斯檢驗fit2 <- eBayes(fit1)# 生成所有基因的檢驗結果報告dif <- topTable(fit2, coef="diseaseprogres. - diseasestable", n=nrow(fit2), lfc=log2(1.5))# 用P.Value進行篩選，得到全部差異表達基因dif <- dif[dif[,"P.Value"]<0.01,]# 顯示一部分報告結果head(dif)

>head(dif)logFCAveExprtP.Valueadj.P.ValB39400_at-0.99978505.634004-5.7273291.482860e-050.10345442.44583541303_at-1.34303064.540225-5.5968131.974284e-050.10345442.154635033791_at1.96479626.8379035.4004993.047498e-050.10345441.713558336131_at0.95742149.9453345.3677413.277762e-050.10345441.639622337636_at2.05340935.4786835.1205195.699454e-050.14391121.078831336122_at0.80086047.1462934.7767211.241402e-040.26121180.2922106

下面逐次介紹這個分析過程的六個關鍵步驟：構建基因表達矩陣、構建實驗設計矩陣、構建對比模型（對比矩陣）、線性模型擬合、貝葉斯檢驗和生成結果報表。

design

image.png

3.2 註釋

# 載入註釋工具包library(annotate)# 獲得基因晶片註釋包名稱affydb <- annPkgName(CLLbatch@annotation,type="db")affydb# 如果沒有安裝，先安裝biocLite(affydb)#　載入該包library(affydb, character.only = TRUE)# 根據每個探針組的ID獲取對應基因Gene Symbol，並作為新的一列dif$symbols<- getSYMBOL(rownames(dif), affydb)# 根據探針ID獲取對應基因Entrez IDdif$EntrezID<- getEG(rownames(dif), affydb)# 顯示前幾行head(dif)

image.png

3.3 統計分析及視覺化

# 載入包library(GOstats)# 提取HG_U95Av2晶片中所有探針組對應的EntrezID，注意保證uniqentrezUniverse <- unique(unlist(mget(rownames(eset), hgu95av2ENTREZID)))# 提取差異表達基因及其對應的EntrezIDentrezSelected <- unique(dif[!is.na(dif$EntrezID),"EntrezID"])# 設定GO富集分析的所有引數params <-new("GOHyperGParams", geneIds=entrezSelected, universeGeneIds=entrezUniverse,              annotation=affydb, ontology="BP", pvalueCutoff=0.001, conditional=FALSE,              testDirection="over")# 對所有的GOterm根據params引數做超幾何檢驗hgOver <- hyperGTest(params)# 生成所有GOterm的檢驗結果報表bp <- summary(hgOver)# 同時生成所有GOterm的檢驗結果檔案，每個GOterm都有指向官方網站的連結，以獲得詳細資訊htmlReport(hgOver, file="ALL_go.html")# 顯示結果前幾行head(bp)

image.png

# 安裝並載入所需包biocLite("GeneAnswers")library(GeneAnswers)# 選取dif中的三列資訊構成新的矩陣，新一列必須是EntrezIDhumanGeneInput <- dif[, c("EntrezID","logFC","P.Value")]# 獲得humanGeneInput中基因的表達值humanExpr <- eset[match(rownames(dif), rownames(eset)), ]humanExpr <- cbind(humanGeneInput[,"EntrezID"], humanExpr)# 去除NA資料humanGeneInput <- humanGeneInput[!is.na(humanGeneInput[,1]),]humanExpr <- humanExpr[!is.na(humanExpr[,1]),]# KEGG通路的超幾何檢驗y <- geneAnswersBuilder(humanGeneInput,"org.Hs.eg.db", categoryType ="KEGG",                        testType ="hyperG", pvalueT=0.1, geneExpressionProfile=humanExpr,                        verbose = FALSE)getEnrichmentInfo(y)[1:6]

image.png

biocLite("pheatmap")library(pheatmap)# 從基因表達矩陣中，選取差異表達基因對應的資料selected <- eset[rownames(dif), ]# 將selected矩陣每行的名稱由探針組ID轉換為對應的基因symbolrownames(selected) <- dif$symbols# 考慮顯示問題，這裡只畫前20個基因的熱圖pheatmap(selected[1:20,], color = colorRampPalette(c("green","black","red"))(100),         fontsize_row = 4, scale ="row", border_color = NA)

image.png

biocLite("Rgraphviz")library(Rgraphviz)# 顯著富集的GO term的DAG關係圖ghandle <- goDag(hgOver)# 這個圖太大，這裡只能選一部分資料構建區域性圖subGHandle <- subGraph(snodes=as.character(summary(hgOver)[,1]), graph = ghandle)plot(subGHandle)

image.png

source("")biocLite("KEGG.db")yy<- geneAnswersReadable(y)geneAnswersConceptNet(yy, colorValueColumn ="logFC", centroidSize ="pvalue", output="interactive")

這裡我報錯了，這個網站解決問題

image.png

yyy<- geneAnswersSort(yy,sortBy ="pvalue")geneAnswersHeatmap(yyy)

作者：蘇慕晨楓
連結：