r語言代寫使用Bioconductor 分析芯片數據

 

原文鏈接：http://tecdat.cn/?p=4764

 

介紹

芯片數據分析流程有些復雜，但使用 R 和 Bioconductor 包進行分析就簡單多了。本教程將一步一步的展示如何安裝 R 和 Bioconductor，通過 GEO 數據庫下載芯片數據， 對數據進行標准化，然后對數據進行質控檢查，最后查找差異表達的基因。教程示例安裝的各種依賴包和運行命令均是是在 Ubuntu 環境中運行的（版本： Ubuntu 10.04， R 2.121）,教程的示例代碼和圖片在這里。

安裝 R 和 Bioconductor 包

打開命令終端，先安裝 R 和 Bioconductor 的依賴包，然后安裝 R.

1.  
   $ sudo apt-get install r-base-core libxml2-dev libcurl4-openssl-dev curl$ R
2.  
    

之后在 R 環境中安裝 Bioconductor 包

> # 下載 Bioconductor 的安裝程序> source("http://bioconductor.org/biocLite.R")> # 安裝 Bioconductor 的核心包> biocLite()> # 安裝 GEO 包> biocLite("GEOquery")

如果你沒有管理員權限，你需要將這些包安裝到你個人庫目錄中。安裝 Bioconductor 需要一段時間，GEOquery 包也需要安裝，GEOquery 是 NCBI 存儲標准化的轉錄組數據的基因表達綜合數據庫 GEO 的接口程序。

下載芯片數據

本教程中我們使用 Dr Andrew Browning 發表的數據集 GSE20986。HUVECs（人臍靜脈內皮細胞）是從人胎兒臍帶血中提取出來的，通常用來研究內皮細胞的病理生理機制。這個實驗設計中，從捐獻者的眼組織中提取的虹膜、視網膜、脈絡膜微血管內皮細胞用來和 HUVEC 細胞進行比對，以便考查 HUVEC 細胞能否替代其他細胞作為研究眼科疾病的樣本。 GEO 頁面同時也包含了關於本實驗研究的其他信息。對於每組樣本有三次測量，樣品分成四組 iris, retina, HUVEC, choroidal 。

實驗平台是 GPL570 -- GEO 數據庫對人類轉錄組芯片 Affymetrix Human Genome U133 Plus 2.0 Array 的縮寫。通過 GSM 鏈接 GSM524662，我們可以得到各個芯片的更多實驗條件信息。同一個實驗設計中的不同芯片的實驗條件應該是相同的。不同細胞系的提取條件，細胞生長條件，RNA 提取程序，RNA 處理程序，RNA 與探針雜交條件，用哪種儀器掃描芯片，這些數據都以符合 MIAME 標准的形式存儲着。

對於每一個芯片，數據表中存儲着探針組和對應探針組標准化之后的基因表達量值。表頭中提供了表達量值標准化所用的算法。本教程中使用 GC-RMA 算法進行數據標准化，關於 GC-RMA 算法的詳細細節可以參考這篇文獻。

首先，我們從 GEO 數據庫下載原始數據，導入 GEOquery 包，用它下載原始數據，下載的原始數據約 53MB。

> library(GEOquery)> getGEOSuppFiles("GSE20986")

下載完成后，打開文件管理器，在啟動 R 程序的文件夾里可以看到當前文件夾下生成了一個 GSE20986 文件夾，可以直接查看文件夾里的內容：

$ ls GSE20986/filelist.txt GSE20986_RAW.tar

GSE20986_RAW.tar 文件是壓縮打包的 CEL（Affymetrix array 數據原始格式）文件。先解包數據，再解壓數據。這些操作可以直接在 R 命令終端中進行：

> untar("GSE20986/GSE20986_RAW.tar", exdir="data")> cels <- pattern="[gz]"> sapply(paste("data", cels, sep="/"), gunzip)> cels

 

芯片實驗信息整理

在對數據進行分析之前，我們需要先整理好實驗設計信息。這其實就是一個文本文件，包含芯片名字、此芯片上雜交的樣本名字。為了方便在 R 中 使用 simpleaffy 包讀取實驗信息文本文件，需要先編輯好格式：

$ ls -1 data/*.CEL > data/phenodata.txt

將這個文本文件用編輯器打開，現在其中只有一列 CEL 文件名，最終的實驗信息文本需要包含三列數據（用 tab 分隔），分別是 Name, FileName, Target。本教程中 Name 和 FileName 這兩欄是相同的，當然 Name 這一欄可以用更加容易理解的名字代替。

Target 這一欄數據是芯片上的樣本標簽，例如 iris, retina, HUVEC, choroidal，這些標簽定義了那些數據是生物學重復以便后面的分析。

這個實驗信息文本文件最終格式是這樣的：

Name FileName TargetGSM524662.CEL GSM524662.CEL irisGSM524663.CEL GSM524663.CEL retinaGSM524664.CEL GSM524664.CEL retinaGSM524665.CEL GSM524665.CEL irisGSM524666.CEL GSM524666.CEL retinaGSM524667.CEL GSM524667.CEL irisGSM524668.CEL GSM524668.CEL choroidGSM524669.CEL GSM524669.CEL choroidGSM524670.CEL GSM524670.CEL choroidGSM524671.CEL GSM524671.CEL huvecGSM524672.CEL GSM524672.CEL huvecGSM524673.CEL GSM524673.CEL huvec

注意：每欄之間是使用 Tab 進行分隔的，而不是空格！

載入數據並對其進行標准化

需要先安裝 simpleaffy 包，simpleaffy 包提供了處理 CEL 數據的程序，可以對 CEL 數據進行標准化同時導入實驗信息（即前一步中整理好的實驗信息文本文件內容），導入數據到 R 變量 celfiles 中：

> biocLite("simpleaffy")> library(simpleaffy)> celfiles <- read.affy(covdesc="phenodata.txt", path="data")

 

你可以通過輸入 ‘celfiles’ 來確定數據導入成功並添加芯片注釋（第一次輸入 ‘celfiles’ 的時候會進行注釋）：

> celfilesAffyBatch objectsize of arrays=1164x1164 features (12 kb)cdf=HG-U133_Plus_2 (54675 affyids)number of samples=12number of genes=54675annotation=hgu133plus2notes=

 

現在我們需要對數據進行標准化，使用 GC-RMA 算法對 GEO 數據庫中的數據進行標准化，第一次運行的時候需要下載一些其他的必要文件：

> celfiles.gcrma <- gcrma(celfiles)Adjusting for optical effect............Done.Computing affinitiesLoading required package: AnnotationDbi.Done.Adjusting for non-specific binding............Done.NormalizingCalculating Expression

 

如果你想看標准化之后的數據，輸入 celfiles.gcrma， 你會發現提示已經不是 AffyBatch object 了，而是 ExpressionSet object，是已經標准化了的數據：

> celfiles.gcrmaExpressionSet (storageMode: lockedEnvironment)assayData: 54675 features, 12 sampleselement names: exprsprotocolDatasampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)varLabels: ScanDatevarMetadata: labelDeionphenoDatasampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)varLabels: sample FileName TargetvarMetadata: labelDeionfeatureData: noneexperimentData: use 'experimentData(object)'Annotation: hgu133plus2

 

數據質量控制

再進行下一步的數據分析之前，我們有必要對數據質量進行檢查，確保沒有其他的問題。首先，可以通過對標准化之前和之后的數據畫箱線圖來檢查 GC-RMA 標准化的效果：

> # 載入色彩包> library(RColorBrewer)> # 設置調色板> cols <-> # 對標准化之前的探針信號強度做箱線圖> boxplot(celfiles, col=cols)> # 對標准化之后的探針信號強度做箱線圖，需要先安裝 affyPLM 包，以便解析 celfiles.gcrma 數據> biocLite("affyPLM")> library(affyPLM)> boxplot(celfiles.gcrma, col=cols)> # 標准化前后的箱線圖會有些變化> # 但是密度曲線圖看起來更直觀一些> # 對標准化之前的數據做密度曲線圖> hist(celfiles, col=cols)> # 對標准化之后的數據做密度曲線圖> hist(celfiles.gcrma, col=cols)

 

數據標准化之前的箱線圖

標准化之前的箱線圖

數據標准化之后的箱線圖

標准化之后的箱線圖

數據標准化之前的密度曲線圖

標准化之前的密度曲線圖

數據標准化之后的密度曲線圖

標准化之后的密度曲線圖

通過這些圖我們可以看出這12張芯片數據之間差異不大，標准化處理將所有芯片信號強度標准化到具有類似分布特征的區間內。為了更詳細地了解芯片探針信號強度，我們可以使用 affyPLM 對單個芯片 CEL 數據進行可視化。

> # 從 CEL 文件讀取探針信號強度:> celfiles.qc <-> # 對芯片 GSM24662.CEL 信號進行可視化:> image(celfiles.qc, which=1, add.legend=TRUE)> # 對芯片 GSM524665.CEL 進行可視化> # 這張芯片數據有些人為誤差> image(celfiles.qc, which=4, add.legend=TRUE)> # affyPLM 包還可以畫箱線圖> # RLE (Relative Log Expression 相對表達量取對數) 圖中> # 所有的值都應該接近於零。 GSM524665.CEL 芯片數據由於有人為誤差而例外。> RLE(celfiles.qc, main="RLE")> # 也可以用 NUSE (Normalised Unscaled Standard Errors)作圖比較.> # 對於絕大部分基因，標准差的中位數應該是1。> # 芯片 GSM524665.CEL 在這個圖中，同樣是一個例外> NUSE(celfiles.qc, main="NUSE")

 

標准的芯片 AffyPLM 信號圖

標准的芯片 AffyPLM 信號圖

存在人工誤差的芯片 AffyPLM 信號圖

存在人工誤差的芯片 AffyPLM 信號圖

CEL 數據的 RLE 圖

RLE (Relative Log Expression) 圖

CEL 數據的 NUSE 圖

NUSE (Normalised Unscaled Standard Errors) 圖

我們還可以通過層次聚類來查看樣本之間的關系：

1.  
   > eset <-> distance <- method="maximum"> clusters <-> plot(clusters)
2.  
    

CEL 數據的層次聚類圖

層次聚類圖

圖形顯示，與其他眼組織相比 HUVEC 樣品是單獨的一組，表現出組織類型聚集的一些特征，另外 GSM524665.CEL 數據在此圖中並不顯示為異常值。

數據過濾

現在我們可以對數據進行分析了，分析的第一步就是要過濾掉數據中的無用數據，例如作為內參的探針數據，基因表達無明顯變化的數據（在差異表達統計時也會被過濾掉），信號值與背景信號差不多的探針數據。下面的 nsFilter 參數是為了不刪除沒有 Entrez Gene ID 的位點，保留有重復 Entrez Gene ID 的位點：

> celfiles.filtered <- require.entrez="FALSE," remove.dupentrez="FALSE)"> # 哪些位點被過濾掉了？為什么？> celfiles.filtered$filter.log$numLowVar[1] 27307$feature.exclude[1] 62

 

我們可以看出有 27307 個探針位點因為無明顯表達差異（LowVar）被過濾掉，有 62 個探針位點因為是內參而被過濾掉。

 

如果您有任何疑問，請在下面發表評論。