1、下載安裝R相關(guān)包
> source("http://bioconductor.org/biocLite.R") # 由于我之前已經(jīng)安裝了，這里提示升級(jí)
# 網(wǎng)上找了下解決方法
> remove.packages("BiocInstaller")
> install.packages("BiocInstaller") # 更新安裝包
# 安裝、更新bioconductor核心安裝包
> biocLite()
# 安裝 GEO 包
> biocLite("GEOquery")
# 出現(xiàn)如下問題
installation path not writeable, unable to update packages: cluster, foreign, MASS, Matrix, nlme, survival
# 解決方法：windows以管理員身份打開，linux建議的解決方法為使用管理者權(quán)限啟動(dòng)R，即 sudo R
GEOquery 是 NCBI 存儲(chǔ)標(biāo)準(zhǔn)化的轉(zhuǎn)錄組數(shù)據(jù)的基因表達(dá)綜合數(shù)據(jù)庫 GEO 的接口程序。
2、下載芯片數(shù)據(jù)
教程中我們使用 Dr Andrew Browning 發(fā)表的數(shù)據(jù)集 GSE20986。HUVECs（人臍靜脈內(nèi)皮細(xì)胞）是從人胎兒臍帶血中提取出來的，通常用來研究?jī)?nèi)皮細(xì)胞的病理生理機(jī)制。這個(gè)實(shí)驗(yàn)設(shè)計(jì)中，從捐獻(xiàn)者的眼組織中提取的虹膜、視網(wǎng)膜、脈絡(luò)膜微血管內(nèi)皮細(xì)胞用來和 HUVEC 細(xì)胞進(jìn)行比對(duì)，以便考查 HUVEC 細(xì)胞能否替代其他細(xì)胞作為研究眼科疾病的樣本。 GEO 頁面同時(shí)也包含了關(guān)于本實(shí)驗(yàn)研究的其他信息。對(duì)于每組樣本有三次測(cè)量，樣品分成四組 iris, retina, HUVEC, choroidal 。
實(shí)驗(yàn)平臺(tái)是 GPL570 -- GEO 數(shù)據(jù)庫對(duì)人類轉(zhuǎn)錄組芯片 Affymetrix Human Genome U133 Plus 2.0 Array 的縮寫。通過 GSM 鏈接 GSM524662，我們可以得到各個(gè)芯片的更多實(shí)驗(yàn)條件信息。同一個(gè)實(shí)驗(yàn)設(shè)計(jì)中的不同芯片的實(shí)驗(yàn)條件應(yīng)該是相同的。不同細(xì)胞系的提取條件，細(xì)胞生長(zhǎng)條件，RNA 提取程序，RNA 處理程序，RNA 與探針雜交條件，用哪種儀器掃描芯片，這些數(shù)據(jù)都以符合 MIAME 標(biāo)準(zhǔn)的形式存儲(chǔ)著。
對(duì)于每一個(gè)芯片，數(shù)據(jù)表中存儲(chǔ)著探針組和對(duì)應(yīng)探針組標(biāo)準(zhǔn)化之后的基因表達(dá)量值。表頭中提供了表達(dá)量值標(biāo)準(zhǔn)化所用的算法。本教程中使用 GC-RMA 算法進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化，關(guān)于 GC-RMA 算法的詳細(xì)細(xì)節(jié)可以參考這篇文獻(xiàn)。
2.1 下載原始數(shù)據(jù)
首先，我們從 GEO 數(shù)據(jù)庫下載原始數(shù)據(jù)，導(dǎo)入 GEOquery 包，用它下載原始數(shù)據(jù)，下載的原始數(shù)據(jù)約 53MB。
# 先設(shè)置我們的工作路徑
> dir.create('D:/TCGA/microarray_analysis',recursive = T)
> setwd('D:/TCGA/microarray_analysis')
> getwd()
[1] "D:/TCGA/microarray_analysis"
# 下載數(shù)據(jù)
> library(GEOquery)
> getGEOSuppFiles("GSE20986")
# 查看下載好的文件
> list.files('GSE20986')
[1] "GSE20986_RAW.tar" # 可以看到下載了一個(gè)打包文件
> setwd('GSE20986')
> getwd()
[1] "D:/TCGA/microarray_analysis/GSE20986"
> untar('GSE20986_RAW.tar',exdir = 'data') # 解包文件到data目錄下
> setwd('data')
> getwd()
[1] "D:/TCGA/microarray_analysis/GSE20986/data" # 當(dāng)前的工作目錄
# 查看解包后的文件，可以看到共12個(gè)壓縮文件，每個(gè)樣本一個(gè)
> list.files()
[1] "GSM524662.CEL.gz" "GSM524663.CEL.gz" "GSM524664.CEL.gz" "GSM524665.CEL.gz" "GSM524666.CEL.gz" "GSM524667.CEL.gz"
[7] "GSM524668.CEL.gz" "GSM524669.CEL.gz" "GSM524670.CEL.gz" "GSM524671.CEL.gz" "GSM524672.CEL.gz" "GSM524673.CEL.gz"
# 將這些文件進(jìn)行解壓縮
> cels <- list.files() ; sapply(cels, gunzip)
GSM524662.CEL.gz GSM524663.CEL.gz GSM524664.CEL.gz GSM524665.CEL.gz GSM524666.CEL.gz GSM524667.CEL.gz GSM524668.CEL.gz
13555726 13555055 13555639 13560122 13555663 13557614 13556090
GSM524669.CEL.gz GSM524670.CEL.gz GSM524671.CEL.gz GSM524672.CEL.gz GSM524673.CEL.gz
13560054 13555971 13554926 13555042 13555290
# 查看解壓后的文件
> list.files()
[1] "GSM524662.CEL" "GSM524663.CEL" "GSM524664.CEL" "GSM524665.CEL" "GSM524666.CEL" "GSM524667.CEL" "GSM524668.CEL"
[8] "GSM524669.CEL" "GSM524670.CEL" "GSM524671.CEL" "GSM524672.CEL" "GSM524673.CEL"
Affymetrix的探針（proble）一般是長(zhǎng)為25堿基的寡聚核苷酸；探針總是以perfect match 和mismatch成對(duì)出現(xiàn)，其信號(hào)值稱為PM和MM，成對(duì)的perfect match 和mismatch有一個(gè)共同的affyID。
CEL文件：信號(hào)值和定位信息。
CDF文件：探針對(duì)在芯片上的定位信息。
當(dāng)然，如果沒有特殊要求，也可直接使用GEOquery包中的getGEO方法直接下載處理好的表達(dá)矩陣：
> gse20986 <- getGEO("GSE20986", destdir=".") # 會(huì)下載處理好的表達(dá)矩陣
3、芯片實(shí)驗(yàn)信息整理
在對(duì)數(shù)據(jù)進(jìn)行分析之前，我們需要先整理好實(shí)驗(yàn)設(shè)計(jì)信息。這其實(shí)就是一個(gè)文本文件，包含芯片名字、此芯片上雜交的樣本名字。為了方便在 R 中 使用 simpleaffy 包讀取實(shí)驗(yàn)信息文本文件，需要先編輯好格式：
# windows下Rstudio里面自帶的terminal具有部分linux的終端功能，可以借用一下
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ ls -1 *.CEL > finame
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ awk 'BEGIN{OFS="\t";print "Name","Filename","Target"}{print $1,$1}' finame > phenodata.txt
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ cat phenodata.txt
Name Filename Target
GSM524662.CEL GSM524662.CEL
GSM524663.CEL GSM524663.CEL
GSM524664.CEL GSM524664.CEL
GSM524665.CEL GSM524665.CEL
GSM524666.CEL GSM524666.CEL
GSM524667.CEL GSM524667.CEL
GSM524668.CEL GSM524668.CEL
GSM524669.CEL GSM524669.CEL
GSM524670.CEL GSM524670.CEL
GSM524671.CEL GSM524671.CEL
GSM524672.CEL GSM524672.CEL
GSM524673.CEL GSM524673.CEL
最終的實(shí)驗(yàn)信息文本需要包含三列數(shù)據(jù)（用 tab 分隔），分別是 Name, FileName, Target。本教程中 Name 和 FileName 這兩欄是相同的，當(dāng)然 Name 這一欄可以用更加容易理解的名字代替。
Target 這一欄數(shù)據(jù)是芯片上的樣本標(biāo)簽，例如 iris, retina, HUVEC, choroidal，這些標(biāo)簽定義了那些數(shù)據(jù)是生物學(xué)重復(fù)以便后面的分析。具體的信息可以在GEO數(shù)據(jù)庫樣本信息中找到，如下圖。
最后的實(shí)驗(yàn)信息文本文件如下：
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ cat phenodata.txt
Name Filename Target
GSM524662.CEL GSM524662.CEL iris
GSM524663.CEL GSM524663.CEL retina
GSM524664.CEL GSM524664.CEL retina
GSM524665.CEL GSM524665.CEL iris
GSM524666.CEL GSM524666.CEL retina
GSM524667.CEL GSM524667.CEL iris
GSM524668.CEL GSM524668.CEL choroid
GSM524669.CEL GSM524669.CEL choroid
GSM524670.CEL GSM524670.CEL choroid
GSM524671.CEL GSM524671.CEL huvec
GSM524672.CEL GSM524672.CEL huvec
GSM524673.CEL GSM524673.CEL huvec
4、載入數(shù)據(jù)并對(duì)其進(jìn)行標(biāo)準(zhǔn)化
需要先安裝 simpleaffy 包，simpleaffy 包提供了處理 CEL 數(shù)據(jù)的程序，可以對(duì) CEL 數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化同時(shí)導(dǎo)入實(shí)驗(yàn)信息（即前一步中整理好的實(shí)驗(yàn)信息文本文件內(nèi)容），導(dǎo)入數(shù)據(jù)到 R 變量 celfiles 中：
# 中間可能會(huì)有報(bào)錯(cuò)，按照?qǐng)?bào)錯(cuò)提示解決即可
> biocLite("simpleaffy")
> library(simpleaffy)
> celfiles <- read.affy(covdesc='phenodata.txt', path='.')
你可以通過輸入 ‘celfiles’ 來確定數(shù)據(jù)導(dǎo)入成功并添加芯片注釋（第一次輸入 ‘celfiles’ 的時(shí)候會(huì)進(jìn)行注釋）：
> celfiles
installing the source package ‘hgu133plus2cdf’
trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/hgu133plus2cdf_2.18.0.tar.gz'
Content type 'application/x-gzip' length 4353300 bytes (4.2 MB)
downloaded 4.2 MB
* installing *source* package 'hgu133plus2cdf' ...
** R
** data
** preparing package for lazy loading
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
** help
*** installing help indices
converting help for package 'hgu133plus2cdf'
finding HTML links ... 好了
geometry html
hgu133plus2cdf html
hgu133plus2dim html
** building package indices
** testing if installed package can be loaded
*** arch - i386
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
*** arch - x64
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
* DONE (hgu133plus2cdf)
In R CMD INSTALL
The downloaded source packages are in
‘C:\Users\sxuan\AppData\Local\Temp\RtmpItQIXU\downloaded_packages’
AffyBatch object
size of arrays=1164x1164 features (21 kb)
cdf=HG-U133_Plus_2 (54675 affyids)
number of samples=12
number of genes=54675
annotation=hgu133plus2
notes=
Warning messages:
1: replacing previous import ‘AnnotationDbi::tail’ by ‘utils::tail’ when loading ‘hgu133plus2cdf’
2: replacing previous import ‘AnnotationDbi::head’ by ‘utils::head’ when loading ‘hgu133plus2cdf’
現(xiàn)在我們需要對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，使用 GC-RMA 算法對(duì) GEO 數(shù)據(jù)庫中的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，第一次運(yùn)行的時(shí)候需要下載一些其他的必要文件：
> celfiles.gcrma <- gcrma(celfiles)
Adjusting for optical effect............Done.
Computing affinities[1] "Checking to see if your internet connection works..."
installing the source package ‘hgu133plus2probe’
trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/hgu133plus2probe_2.18.0.tar.gz'
... ... # 這里省略了一些下載信息
.Done.
Adjusting for non-specific binding............Done.
Normalizing
Calculating Expression
如果你想看標(biāo)準(zhǔn)化之后的數(shù)據(jù)，輸入 celfiles.gcrma， 你會(huì)發(fā)現(xiàn)提示已經(jīng)不是 AffyBatch object 了，而是 ExpressionSet object，是已經(jīng)標(biāo)準(zhǔn)化了的數(shù)據(jù)：
> celfiles.gcrma
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 12 samples
element names: exprs
protocolData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)
varLabels: ScanDate
varMetadata: labelDescription
phenoData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)
varLabels: sample Filename Target
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu133plus2
5、芯片數(shù)據(jù)可視化
再進(jìn)行下一步的數(shù)據(jù)分析之前，我們有必要對(duì)數(shù)據(jù)質(zhì)量進(jìn)行檢查，確保沒有其他的問題。首先，可以通過對(duì)標(biāo)準(zhǔn)化之前和之后的數(shù)據(jù)畫箱線圖來檢查 GC-RMA 標(biāo)準(zhǔn)化的效果：
> # 載入色彩包
> library(RColorBrewer)
> # 設(shè)置調(diào)色板
> cols <- brewer.pal(8, "Set1")
> # 對(duì)標(biāo)準(zhǔn)化之前的探針信號(hào)強(qiáng)度做箱線圖
> boxplot(celfiles, col=cols)
> # 對(duì)標(biāo)準(zhǔn)化之后的探針信號(hào)強(qiáng)度做箱線圖，需要先安裝 affyPLM 包，以便解析 celfiles.gcrma 數(shù)據(jù)
> biocLite("affyPLM")
> library(affyPLM)
> boxplot(celfiles.gcrma, col=cols)
> # 標(biāo)準(zhǔn)化前后的箱線圖會(huì)有些變化
> # 但是密度曲線圖看起來更直觀一些
> # 對(duì)標(biāo)準(zhǔn)化之前的數(shù)據(jù)做密度曲線圖
> hist(celfiles, col=cols)
> # 對(duì)標(biāo)準(zhǔn)化之后的數(shù)據(jù)做密度曲線圖
> hist(celfiles.gcrma, col=cols)
通過這些圖我們可以看出這12張芯片數(shù)據(jù)之間差異不大，標(biāo)準(zhǔn)化處理將所有芯片信號(hào)強(qiáng)度標(biāo)準(zhǔn)化到具有類似分布特征的區(qū)間內(nèi)。為了更詳細(xì)地了解芯片探針信號(hào)強(qiáng)度，我們可以使用 affyPLM 對(duì)單個(gè)芯片 CEL 數(shù)據(jù)進(jìn)行可視化。
> # 從 CEL 文件讀取探針信號(hào)強(qiáng)度:
> celfiles.qc <- fitPLM(celfiles)
> # 對(duì)芯片 GSM24662.CEL 信號(hào)進(jìn)行可視化:
> image(celfiles.qc, which=1, add.legend=TRUE)
> # 對(duì)芯片 GSM524665.CEL 進(jìn)行可視化
> # 這張芯片數(shù)據(jù)有些人為誤差
> image(celfiles.qc, which=4, add.legend=TRUE)
> # affyPLM 包還可以畫箱線圖
> # RLE (Relative Log Expression 相對(duì)表達(dá)量取對(duì)數(shù)) 圖中
> # 所有的值都應(yīng)該接近于零。 GSM524665.CEL 芯片數(shù)據(jù)由于有人為誤差而例外。
> RLE(celfiles.qc, main="RLE")
> # 也可以用 NUSE (Normalised Unscaled Standard Errors)作圖比較.
> # 對(duì)于絕大部分基因，標(biāo)準(zhǔn)差的中位數(shù)應(yīng)該是1。
> # 芯片 GSM524665.CEL 在這個(gè)圖中，同樣是一個(gè)例外
> NUSE(celfiles.qc, main="NUSE")
我們還可以通過層次聚類來查看樣本之間的關(guān)系：
> eset <- exprs(celfiles.gcrma)
> distance <- dist(t(eset),method="maximum")
> clusters <- hclust(distance)
> plot(clusters)
圖形顯示，與其他眼組織相比 HUVEC 樣品是單獨(dú)的一組，表現(xiàn)出組織類型聚集的一些特征，另外 GSM524665.CEL 數(shù)據(jù)在此圖中并不顯示為異常值。
6、數(shù)據(jù)質(zhì)量控制
現(xiàn)在我們可以對(duì)數(shù)據(jù)進(jìn)行分析了，分析的第一步就是要過濾掉數(shù)據(jù)中的無用數(shù)據(jù)，例如作為內(nèi)參的探針數(shù)據(jù)，基因表達(dá)無明顯變化的數(shù)據(jù)（在差異表達(dá)統(tǒng)計(jì)時(shí)也會(huì)被過濾掉），信號(hào)值與背景信號(hào)差不多的探針數(shù)據(jù)。
下面的 nsFilter 參數(shù)是為了不刪除沒有 Entrez Gene ID 的位點(diǎn)，保留有重復(fù) Entrez Gene ID 的位點(diǎn)：
> celfiles.filtered <- nsFilter(celfiles.gcrma, require.entrez=FALSE, remove.dupEntrez=FALSE)
> # 哪些位點(diǎn)被過濾掉了？為什么？
> celfiles.filtered$filter.log
$numLowVar
[1] 27307
$feature.exclude
[1] 62
我們可以看出有 27307 個(gè)探針位點(diǎn)因?yàn)闊o明顯表達(dá)差異（LowVar）被過濾掉，有 62 個(gè)探針位點(diǎn)因?yàn)槭莾?nèi)參而被過濾掉。
7、查找有表達(dá)差異的探針位點(diǎn)
現(xiàn)在有了過濾之后的數(shù)據(jù)，我們就可以用 limma 包進(jìn)行差異表達(dá)分析了。首先，我們要提取樣本的信息：
> samples <- celfiles.gcrma$Target
> # 檢查數(shù)據(jù)的分組信息
> samples
[1] "iris" "retina" "retina" "iris" "retina" "iris" "choroid"
[8] "choroid" "choroid" "huvec" "huvec" "huvec"
> # 將分組數(shù)據(jù)轉(zhuǎn)換為因子類型變量
> samples <- as.factor(samples)
> # 檢查轉(zhuǎn)換的因子變量
> samples
[1] iris retina retina iris retina iris choroid choroid choroid
[10] huvec huvec huvec
Levels: choroid huvec iris retina
> # 設(shè)置實(shí)驗(yàn)分組
> design <- model.matrix(~0 + samples)
> colnames(design) <- c("choroid", "huvec", "iris", "retina")
> # 檢查實(shí)驗(yàn)分組
> design
choroid huvec iris retina
1 0 0 1 0
2 0 0 0 1
3 0 0 0 1
4 0 0 1 0
5 0 0 0 1
6 0 0 1 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
10 0 1 0 0
11 0 1 0 0
12 0 1 0 0
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$samples
[1] "contr.treatment"
現(xiàn)在我們將芯片數(shù)據(jù)進(jìn)行了標(biāo)準(zhǔn)化和過濾，也有樣品分組和實(shí)驗(yàn)分組信息，可以將數(shù)據(jù)導(dǎo)入 limma 包進(jìn)行差異表達(dá)分析了：
> library(limma)
> # 將線性模型擬合到過濾之后的表達(dá)數(shù)據(jù)集上
> fit <- lmFit(exprs(celfiles.filtered$eset), design)
> # 建立對(duì)比矩陣以比較組織和細(xì)胞系
> contrast.matrix <- makeContrasts(huvec_choroid = huvec - choroid, huvec_retina = huvec - retina, huvec_iris <- huvec - iris, levels=design)
> # 檢查對(duì)比矩陣
> contrast.matrix
Contrasts
Levels huvec_choroid huvec_retina huvec_iris
choroid -1 0 0
huvec 1 1 1
iris 0 0 -1
retina 0 -1 0
> # 現(xiàn)在將對(duì)比矩陣與線性模型擬合，比較不同細(xì)胞系的表達(dá)數(shù)據(jù)
> huvec_fits <- contrasts.fit(fit, contrast.matrix)
> # 使用經(jīng)驗(yàn)貝葉斯算法計(jì)算差異表達(dá)基因的顯著性
> huvec_ebFit <- eBayes(huvec_fits)
> # 返回對(duì)應(yīng)比對(duì)矩陣 top 10 的結(jié)果
> # coef=1 是 huvec_choroid 比對(duì)矩陣, coef=2 是 huvec_retina 比對(duì)矩陣
> topTable(huvec_ebFit, number=10, coef=1)
ID logFC AveExpr t P.Value adj.P.Val
6147 204779_s_at 7.367947 4.171874 72.77016 3.290669e-15 8.985500e-11
7292 207016_s_at 6.934796 4.027229 57.37259 3.712060e-14 5.068076e-10
8741 209631_s_at 5.193313 4.003541 51.22423 1.177337e-13 1.071612e-09
26309 242809_at 6.433514 4.167462 48.52518 2.042404e-13 1.394247e-09
6828 205893_at 4.480463 3.544350 40.59376 1.253534e-12 6.845801e-09
20232 227377_at 3.670688 3.209217 36.03427 4.200854e-12 1.911809e-08
6222 204882_at -5.351976 6.512018 -34.70239 6.154318e-12 2.400711e-08
27051 38149_at -5.051906 6.482418 -31.44098 1.672263e-11 5.282592e-08
6663 205576_at 6.586372 4.139236 31.31584 1.741131e-11 5.282592e-08
6589 205453_at 3.623706 3.210306 30.72793 2.109110e-11 5.759136e-08
B
6147 20.25563
7292 19.44681
8741 18.96282
26309 18.70767
6828 17.75331
20232 17.01960
6222 16.77196
27051 16.08854
6663 16.05990
6589 15.92277
分析結(jié)果的各列數(shù)據(jù)含義：
第一列是探針組在表達(dá)矩陣中的行號(hào)； # 有的沒有此列
第二列“ID” 是探針組的 AffymatrixID；
第三列“l(fā)ogFC”是兩組表達(dá)值之間以2為底對(duì)數(shù)化的的變化倍數(shù)（Fold change, FC），由于基因表達(dá)矩陣本身已經(jīng)取了對(duì)數(shù)，這里實(shí)際上只是兩組基因表達(dá)值均值之差；
第四列“AveExpr”是該探針組所在所有樣品中的平均表達(dá)值；
第五列“t”是貝葉斯調(diào)整后的兩組表達(dá)值間 T 檢驗(yàn)中的 t 值；
第六列“P.Value”是貝葉斯檢驗(yàn)得到的 P 值；
第七列“adj.P.Val”是調(diào)整后的 P 值；
第八列“B”是經(jīng)驗(yàn)貝葉斯得到的標(biāo)準(zhǔn)差的對(duì)數(shù)化值。
如果要設(shè)置一個(gè)倍數(shù)變化閾值，并查看不同閾值返回了多少基因，可以使用 topTable 的 lfc 參數(shù)，參數(shù)設(shè)置為 5,4,3,2 時(shí)返回的基因個(gè)數(shù)：
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=5))
[1] 88
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=4))
[1] 194
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=3))
[1] 386
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=2))
[1] 1016
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=1))
[1] 4025
> probeset.list <- topTable(huvec_ebFit, coef=1, number=10000, lfc=4)
> head(probeset.list)
logFC AveExpr t P.Value adj.P.Val B
204779_s_at 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
204882_at -5.353547 6.513690 -34.70646 6.141288e-12 2.395629e-08 16.77396
9、注釋差異分析結(jié)果的基因ID
為了將探針集注釋上基因 ID （轉(zhuǎn)換為geneID）我們需要先安裝一些數(shù)據(jù)庫的包和注釋的包，之后可以提取 topTable 中的探針 ID 并注釋上基因 ID：
> biocLite("hgu133plus2.db")
> library(hgu133plus2.db)
> library(annotate)
> gene.symbols <- getSYMBOL(rownames(probeset.list), "hgu133plus2")
> results <- cbind(gene.symbols, probeset.list)
> head(results)
gene.symbols logFC AveExpr t P.Value adj.P.Val B
204779_s_at HOXB7 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at <NA> 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at <NA> 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at IL1RL1 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at NLGN1 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
204882_at ARHGAP25 -5.353547 6.513690 -34.70646 6.141288e-12 2.395629e-08 16.77396
# 可以看到這里有一些芯片ID已經(jīng)找不到對(duì)應(yīng)的基因ID了，我們將這些失敗的提出來先
> trans.failed <- results[is.na(results$gene.symbols),]
> head(trans.failed)
gene.symbols logFC AveExpr t P.Value adj.P.Val B
207016_s_at <NA> 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.449872
209631_s_at <NA> 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.966602
243154_at <NA> 4.213968 3.476199 28.24818 4.932212e-11 1.035992e-07 15.294507
209994_s_at <NA> -5.949092 7.766907 -19.86274 1.695876e-09 1.403261e-06 12.349999
1556314_a_at <NA> -6.200352 7.680901 -18.81114 2.916279e-09 2.041844e-06 11.859483
240466_at <NA> 4.671497 4.946474 14.78163 3.156679e-08 1.134161e-05 9.610741
# 看下有多少行失敗的
> nrow(trans.failed)
[1] 16
下面我們使用另一種方法來進(jìn)行注釋探針I(yè)D，雖然也是使用annotate包，這里雖然探針都注釋到的比上面多，但有的一個(gè)探針注釋到了多個(gè)ID:
# 下載芯片平臺(tái)對(duì)應(yīng)的數(shù)據(jù)包
> biocLite("hgu133plus2.db")
> library("annotate")
> probes <- as.character(rownames(probeset.list))
> out <- select(hgu133plus2.db, probes,"SYMBOL")
'select()' returned 1:many mapping between keys and columns # 表示有的探針注釋到了多個(gè)基因ID
> head(out)
PROBEID SYMBOL
1 204779_s_at HOXB7
2 207016_s_at ALDH1A2
3 207016_s_at LOC101928635
4 209631_s_at GPR37
5 209631_s_at SEL1L2
6 242809_at IL1RL1
> probeset.list$PROBEID <- rownames(probeset.list)
> results.final <- merge(out, probeset.list, by="PROBEID", sort=F)
> head(results.final)
PROBEID SYMBOL logFC AveExpr t P.Value adj.P.Val B
1 204779_s_at HOXB7 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
2 207016_s_at ALDH1A2 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
3 207016_s_at LOC101928635 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
4 209631_s_at GPR37 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
5 209631_s_at SEL1L2 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
6 242809_at IL1RL1 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
> nrow(probeset.list);nrow(results.final)
[1] 194
[1] 222
查看一下兩種方法注釋到結(jié)果的區(qū)別：
> dupID <- results.final[duplicated(results.final$PROBEID),]; nrow(dupID)
[1] 28
> b<- results.final[is.na(results.final$SYMBOL),]; b
PROBEID SYMBOL logFC AveExpr t P.Value adj.P.Val B
13 243154_at <NA> 4.213968 3.476199 28.24818 4.932212e-11 1.035992e-07 15.294507
32 1556314_a_at <NA> -6.200352 7.680901 -18.81114 2.916279e-09 2.041844e-06 11.859483
49 240466_at <NA> 4.671497 4.946474 14.78163 3.156679e-08 1.134161e-05 9.610741
# 從上面可以看到仍然是有三個(gè)探針沒有注釋到
# 將第二種方法一對(duì)多注釋的和未注釋到的合并起來
> a<-as.factor(dupID$PROBEID); f <- as.factor(dupID$PROBEID)
> c <- c(f, as.character(b$PROBEID))
> d<-data.frame(c);d
c
1 204670_x_at
2 207016_s_at
3 208306_x_at
4 209312_x_at
5 209631_s_at
6 209994_s_at
7 211796_s_at
8 215193_x_at
9 218805_at
10 221974_at
11 226558_at
12 229748_x_at
13 64064_at
14 243154_at
15 1556314_a_at
16 240466_at
比較后發(fā)現(xiàn)這兩種方法出現(xiàn)問題的16個(gè)探針是相同的，因此直接使用第一種方法getSYMBOL()更加簡(jiǎn)單。
最后，將注釋好的差異表達(dá)基因?qū)懭胛募M(jìn)行后續(xù)分析。
# 將為注釋到的探針移除
> nrow(results);results<-na.omit(results);nrow(results)
[1] 194
[1] 178
> write.table(results, file="resluts.txt", sep="\t",row.names = F)
參考：
http://bioinformatics.knowledgeblog.org/2011/06/20/analysing-microarray-data-in-bioconductor/
Bioconductor 分析芯片數(shù)據(jù)教程
