使用oligo軟件包處理芯片數據

使用oligo軟件包處理芯片數據
原創金子哦 最后發布於2017-03-27 18:15:10 閱讀數 11376 收藏
展開
本博客介紹過 Affy芯片的處理方法 ，其中所使用的軟件包有一定的局限性，無法讀取和分析一些新版Affy芯片。本文介紹oligo軟件包的處理方法以解決這些問題。oligo軟件包並不是新出現的軟件包，只因新類型芯片的不斷推出，關注它的用戶越來越多。而且，除了用於Affy芯片處理外，oligo軟件包還可處理NimbleGen芯片。

oligo處理芯片的原理和其他方法相同，難點在最后一步：從探針到基因注釋的轉換。

 

1 准備工作
1.1 數據下載
本文還是以Affy芯片為例，介紹oligo軟件包處理芯片數據的方法。所用數據集為NCBI GSE72154，可通過HTTP或FTP下載原始數據文件。建議使用FTP下載，登錄ftp.ncbi.nlm.nih.gov站點，在/geo/datasets下找到GSE72154子目錄后全部下載。下載的目錄中包含matrix、miniml、soft和suppl四個子目錄。

1.2 軟件包安裝
oligo是bioconductor軟件包，可以通過下面語句安裝：

source("https://bioconductor.org/biocLite.R")
biocLite("oligo")
2 文件讀取與信息修改
oligo讀取芯片原始數據文件的函數為 read.celfiles 或 read.xysfiles (NimbleGen芯片)，可以直接讀取壓縮文件，一張芯片對應一個文件。使用FTP下載后CEL文件位於suppl子目錄下，把GSE72154_RAW.tar解壓縮為gz文件或CEL文件即可。

library(oligo)
data.dir <- "GSE72154/suppl/"
(celfiles <- list.files(data.dir, "\\.gz$"))
## [1] "GSM1856253_Bohmer_685-3_Ler-0_Rep1_AraGene-1-1-st.CEL.gz"
## [2] "GSM1856254_Bohmer_685-4_Ler-0_Rep2_AraGene-1-1-st.CEL.gz"
## [3] "GSM1856255_Bohmer_685-1_wrky43_Rep1_AraGene-1-1-st.CEL.gz"
## [4] "GSM1856256_Bohmer_685-2_wrky43_Rep2_AraGene-1-1-st.CEL.gz"
data.raw <- read.celfiles(filenames = file.path(data.dir, celfiles))
## Reading in : GSE72154/suppl//GSM1856253_Bohmer_685-3_Ler-0_Rep1_AraGene-1-1-st.CEL.gz
## Reading in : GSE72154/suppl//GSM1856254_Bohmer_685-4_Ler-0_Rep2_AraGene-1-1-st.CEL.gz
## Reading in : GSE72154/suppl//GSM1856255_Bohmer_685-1_wrky43_Rep1_AraGene-1-1-st.CEL.gz
## Reading in : GSE72154/suppl//GSM1856256_Bohmer_685-2_wrky43_Rep2_AraGene-1-1-st.CEL.gz

read.celfiles 函數讀取數據的過程中會檢查系統中是否已安裝相應的芯片注釋軟件包，沒有安裝的話將有警告。如果在讀取的過程中沒有自動安裝相應注釋軟件包，請手動安裝后再重新讀取數據。

本例含4張芯片，文件名上面代碼已列出，含野生型和突變體各兩個重復。樣品名稱默認為文件名，需簡單修改：

treats <- strsplit("Ler Ler wrky wrky", " ")[[1]]
(snames <- paste(treats, 1:2, sep = ""))
## [1] "Ler1" "Ler2" "wrky1" "wrky2"
sampleNames(data.raw) <- snames
pData(data.raw)$index <- treats
sampleNames(data.raw)
## [1] "Ler1" "Ler2" "wrky1" "wrky2"
pData(data.raw)
## index
## Ler1 Ler
## Ler2 Ler
## wrky1 wrky
## wrky2 wrky
樣品名稱與芯片文件必需對應，不能有重復；pData的index數據表示處理類型。

3 探針水平分析
原始芯片數據讀入后應簡單評估一下芯片的質量，這在 其他文章 中介紹過，oligo中也提供了很多方法，如果樣品重復量足夠，boxplot是比較直觀的方法：

fit1 <- fitProbeLevelModel(data.raw)
boxplot(fit1, names = NA, col = as.factor(treats))
legend("topright", legend = unique(treats), fill = as.factor(unique(treats)),
box.col = NA, bg = "white", inset = 0.01)

 

當然也可以使用MAplot：

par(mfrow = c(2, 2))
MAplot(data.raw[, 1:4], pairs = FALSE)

4 表達量計算
含背景處理、歸一化和表達量計算三個步驟，可選組合方法很多。下面就用通用的三合一處理函數rma一步完成：

data.eset <- rma(data.raw)
## Background correcting
## Normalizing
## Calculating Expression

得到eset類型的數據，用exprs函數可以提取其中的表達量矩陣：

data.exprs <- exprs(data.eset)
str(data.exprs)
## num [1:38408, 1:4] 6.52 6.67 6.52 6.21 7.07 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:38408] "13320001" "13320003" "13320005" "13320007" ...
## ..$ : chr [1:4] "Ler1" "Ler2" "wrky1" "wrky2"
head(data.exprs)
## Ler1 Ler2 wrky1 wrky2
## 13320001 6.515240 6.975158 6.163749 6.121993
## 13320003 6.670580 6.816606 5.808836 5.653102
## 13320005 6.522566 7.396850 7.011411 6.604062
## 13320007 6.206310 6.538738 5.477303 6.572429
## 13320009 7.071985 7.015605 6.426396 6.932439
## 13320011 6.558686 6.333987 6.670292 6.277999
5 P/A過濾
目的是去除“不表達”的基因/探針數據，使用paCalls函數，選取p值小於0.05的探針：

xpa <- paCalls(data.raw)
head(xpa)
## Ler1 Ler2 wrky1 wrky2
## 7 0.477083333 0.80416667 0.327083333 0.689583333
## 8 0.550619835 0.73347107 0.786157025 0.361570248
## 11 0.327731092 0.01890756 0.231092437 0.189075630
## 14 0.005138746 0.01130524 0.005138746 0.005138746
## 17 0.057851240 0.07231405 0.045454545 0.027892562
## 20 0.013360740 0.01027749 0.051387461 0.006166495
AP <- apply(xpa, 1, function(x) any(x < 0.05))
xids <- as.numeric(names(AP[AP]))
head(xids)
## [1] 11 14 17 20 28 42

計算后發現P/A結果和表達量分析結果所用的探針名稱是不一樣的，需要使用探針信息進行轉換。getProbeInfo函數可以獲取相關信息：

pinfo <- getProbeInfo(data.raw)
head(pinfo)
## fid man_fsetid
## 1 7 13507007
## 2 8 13379105
## 3 11 13402401
## 4 14 13363588
## 5 17 13492405
## 6 20 13488138

兩類id轉換后進行篩選：

fids <- pinfo[pinfo$fid %in% xids, 2]
head(fids)
## [1] "13402401" "13363588" "13492405" "13488138" "13349578" "13529865"
nrow(data.exprs)
## [1] 38408
data.exprs <- data.exprs[rownames(data.exprs) %in% fids, ]
nrow(data.exprs)
## [1] 35697

第一個語句篩選出表達基因的探針ID，第四個語句用篩選出的ID過濾掉不表達的基因/探針。

6 獲取差異表達探針
下面使用的limma方法，更多方法請參考本博客 其他文章 。先根據樣品名稱、類型和順序構建試驗設計矩陣：

library(limma)
(design <- model.matrix(~0 + treats))
## treatsLer treatswrky
## 1 1 0
## 2 1 0
## 3 0 1
## 4 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$treats
## [1] "contr.treatment"
colnames(design) <- gsub("treats", "", colnames(design))
rownames(design) <- snames
design
## Ler wrky
## Ler1 1 0
## Ler2 1 0
## wrky1 0 1
## wrky2 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$treats
## [1] "contr.treatment"

然后將試驗設計矩陣應用於表達量數據（exprs）或eset數據，獲取擬合模型：

fit1 <- lmFit(data.exprs, design)
## fit1 <- lmFit(data.eset, design)

接下來進行試驗處理組間的比較。因為只有兩組樣品，只可能是突變體和野生型間的比較：

contrast.matrix <- makeContrasts("wrky-Ler", levels = design)
fit2 <- contrasts.fit(fit1, contrast.matrix)
fit2 <- eBayes(fit2)

使用topTable函數獲取表達變化倍數和相關統計值：

limma.all <- topTable(fit2, coef = 1, adjust.method = "fdr", number = 3e+05)
limma.fc2 <- topTable(fit2, coef = 1, adjust.method = "fdr", lfc = 1, number = 3e+05)
nrow(limma.all)
## [1] 35697
nrow(limma.fc2)
## [1] 1028
head(limma.fc2)
## logFC AveExpr t P.Value adj.P.Val B
## 13356037 4.540665 10.780301 39.28527 1.070928e-06 0.03822893 3.571807
## 13348341 4.206080 6.004308 31.68465 2.710592e-06 0.04838000 3.402196
## 13348348 3.569355 5.319289 21.22435 1.522879e-05 0.09541435 2.882561
## 13400731 3.029046 6.347611 20.99470 1.595751e-05 0.09541435 2.863891
## 13474474 3.259238 5.847156 20.82997 1.650674e-05 0.09541435 2.850205
## 13454438 2.750208 4.778045 19.55321 2.165705e-05 0.09541435 2.735128

topTable函數中可以直接使用p.value產生進行p值過濾，但由於本數據集生物學重復數太少，上面代碼僅用表達變化倍數進行過濾。number參數用於限定結果的最大行數。

7 探針到基因的轉換
有多種方法，下面介紹兩種。

7.1 方法1：使用芯片注釋文件
每版芯片都會有專門文件對探針進行說明，這些文件多數是公開的。前面通過NCBI FTP下載GEO芯片數據集時已提到matrix、miniml、soft和suppl四個子目錄。miniml和soft目錄都包含有平台相關的信息文件，其中soft目錄下的GSE72154_family.soft.gz文件含有表頭，可以使用該文件：

ff <- "GSE72154/soft/GSE72154_family.soft.gz"
nn <- grep("^[^#!^]", readLines(ff))[1] - 1
pfinfo <- read.table(ff, sep = "\t", quote = "", header = TRUE, skip = nn, fill = TRUE)
colnames(pfinfo)
## [1] "ID" "seqname"
## [3] "strand" "start"
## [5] "stop" "total_probes"
## [7] "gene_assignment" "mrna_assignment"
## [9] "GB_ACC" "ORF"
## [11] "SPOT_ID" "swissprot"
## [13] "GO_biological_process" "GO_cellular_component"
## [15] "GO_molecular_function" "protein_domains"
## [17] "crosshyb_type" "category"

可以看到表格各列所包含的內容。現階段我們只需要ID和ORF列，即探針編號和它們所檢測的mRNA。提取信息后做一些必要的處理

pfinfo <- pfinfo[, c(1, 10)]
pfinfo$ORF <- toupper(pfinfo$ORF)
## 刪除重復信息
pfinfo <- pfinfo[!duplicated(pfinfo), ]
## 刪除非mRNA探針
pfinfo <- pfinfo[pfinfo$ORF != "", ]
rownames(pfinfo) <- pfinfo$ID
nrow(pfinfo)
## [1] 20657

得到以上基本信息后就可以添加到差異表達列表中。首先去除沒有mRNA信息的探針：

result1 <- limma.fc2[rownames(limma.fc2) %in% pfinfo$ID, ]
nrow(result1)
## [1] 214

把AGI信息添加到表達列表，保存表格文件即可：

result1$agi <- pfinfo[rownames(result1), 2]
head(result1)
## logFC AveExpr t P.Value adj.P.Val B
## 13356037 4.540665 10.780301 39.28527 1.070928e-06 0.03822893 3.571807
## 13474474 3.259238 5.847156 20.82997 1.650674e-05 0.09541435 2.850205
## 13454438 2.750208 4.778045 19.55321 2.165705e-05 0.09541435 2.735128
## 13335452 -2.729114 9.397657 -18.21499 2.935056e-05 0.09541435 2.594936
## 13376651 -2.154874 7.127874 -16.55810 4.414524e-05 0.09541435 2.387110
## 13464968 1.667566 8.680761 14.94084 6.845745e-05 0.09541435 2.138016
## agi
## 13356037 AT1G67105.1
## 13474474 AT4G27940.1
## 13454438 AT3G27473.1
## 13335452 AT1G03880.1
## 13376651 AT1G47540.2
## 13464968 AT4G02555.1
write.csv(result1, "results.fc2.1.csv")
7.2 方法2：使用基因組特征文件
芯片廠商提供的探針注釋文件有可能過時。如分析十多年前的試驗芯片，一些基因的注釋可能有很大變化。但不管怎么樣，探針在基因組上的位置信息應該不會有太大的差異，應用基因組進行比對總不會大錯。

Bioconductor提供了專門用於處理基因組特征文件的類和函數，相關的gff格式文件可以從各種途徑免費獲取。本例使用的是TAIR網站下載的擬南芥基因組特征文件。

library("rtracklayer")
ath.gff <- readGFFAsGRanges("TAIR10_GFF3_genes.gff")
class(ath.gff)
## [1] "GRanges"
## attr(,"package")
## [1] "GenomicRanges"

讀入數據存儲為GRanges類，其數據結構可以使用str函數查看，相關的基礎知識已在 R/BioC序列處理之五：Rle和Ranges 一文中介紹過。GRanges的每一條信息都包含了很多內容，這里僅關心外顯子、UTR和miRNA相關的序列：

xtype <- as.character(ath.gff@elementMetadata@listData$type)
unique(xtype)
## [1] "chromosome" "gene"
## [3] "mRNA" "protein"
## [5] "exon" "five_prime_UTR"
## [7] "CDS" "three_prime_UTR"
## [9] "miRNA" "tRNA"
## [11] "ncRNA" "pseudogene"
## [13] "pseudogenic_transcript" "pseudogenic_exon"
## [15] "transposable_element_gene" "mRNA_TE_gene"
## [17] "snoRNA" "snRNA"
## [19] "rRNA"
sels <- grepl("exon|utr|miRNA", xtype, ignore.case = TRUE)
ath.gff <- ath.gff[sels]

接下來獲取探針在基因組上的位置信息。這一步可以用前面的芯片注釋文件數據，也可以使用隨CEL文件讀取時載入的數據。類似數據在前面P/A過濾部分已經用到過，但這里要提取更多的信息：

availProbeInfo(data.raw)
## $pmfeature
## [1] "fid" "fsetid" "atom" "x" "y"
##
## $featureSet
## [1] "fsetid" "strand" "start" "stop"
## [5] "exon_id" "crosshyb_type" "level" "chrom"
## [9] "type"
##
## $core_mps
## [1] "meta_fsetid" "transcript_cluster_id" "fsetid"
pinfo <- getProbeInfo(data.raw, field = strsplit("fid chrom strand start stop exon_id",
" ")[[1]])
head(pinfo)
## fid man_fsetid strand start stop exon chrom
## 1 7 13507007 0 11107381 11107491 884614 chr5
## 2 8 13379105 1 19981386 19981486 778732 chr1
## 3 11 13402401 0 16131514 16131647 798181 chr2
## 4 14 13363588 1 1457201 1457469 765883 chr1
## 5 17 13492405 1 15589898 15590906 872458 chr4
## 6 20 13488138 1 11944077 11944154 868967 chr4

在應用數據前做一些必要的處理：

pinfo <- pinfo[grepl("^chr", pinfo$chrom), ]
pinfo$chrom <- gsub("^c", "C", pinfo$chrom)
pinfo$strand <- sapply(pinfo$strand, FUN = function(x) {
if (x == 0)
"+" else "-"
})

然后構建探針對應的GRanges類對象：

xrange <- GRanges(seqnames = pinfo$chrom, ranges = IRanges(start = pinfo$start,
end = pinfo$stop, names = pinfo$man_fsetid), strand = pinfo$strand)

正式比對並將產生探針-基因注釋數據：

xlap <- findOverlaps(xrange, ath.gff, select = "all", type = "within")
anno <- data.frame(fid = names(xrange)[queryHits(xlap)], agi = ath.gff@elementMetadata$Parent@unlistData[subjectHits(xlap)])
anno <- anno[!duplicated(anno), ]
head(anno)
## fid agi
## 1 13507007 AT5G29044.1
## 2 13379105 AT1G53541.1
## 3 13402401 AT2G38544.1
## 4 13363588 AT1G05070.1
## 5 13492405 AT4G32290.1
## 6 13488138 AT4G22740.1

有部分探針可能和基因組上多個基因對應，有必要進行合並，或者舍棄這部分數據：

anno <- aggregate(agi ~ fid, FUN = c, data = anno, simplify = TRUE)
anno$agi <- sapply(anno$agi, FUN = paste, collapse = ";")
rownames(anno) <- anno$fid

最后合並數據：

result2 <- limma.fc2[rownames(limma.fc2) %in% anno$fid, ]
result2$agi <- sapply(anno[rownames(result2), "agi"], FUN = paste, collapse = ";")
nrow(limma.fc2)
## [1] 1028
nrow(result2)
## [1] 464
write.csv(result2, "results.fc2.2.csv")

為什么最后得到的基因數量比前一種方法要多？自己分析原因吧。

如果需要進行GO/KEGG途徑富集分析，還應導出芯片中的“表達”基因：

xagis <- anno[anno$fid %in% fids, "agi"]
writeLines(xagis, "results.present.genes.txt")
8 SessionInfo
print(sessionInfo(), locale = FALSE)
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 8 (jessie)
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] rtracklayer_1.34.2 GenomicRanges_1.26.4
## [3] GenomeInfoDb_1.10.3 limma_3.30.13
## [5] pd.aragene.1.0.st_3.12.0 DBI_0.6
## [7] RSQLite_1.1-2 oligo_1.38.0
## [9] Biostrings_2.42.1 XVector_0.14.1
## [11] IRanges_2.8.2 S4Vectors_0.12.2
## [13] Biobase_2.34.0 oligoClasses_1.36.0
## [15] BiocGenerics_0.20.0 zblog_0.1.0
## [17] knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.10 BiocInstaller_1.24.0
## [3] formatR_1.4 highr_0.6
## [5] bitops_1.0-6 iterators_1.0.8
## [7] tools_3.3.3 zlibbioc_1.20.0
## [9] digest_0.6.12 bit_1.1-12
## [11] evaluate_0.10 memoise_1.0.0
## [13] preprocessCore_1.36.0 lattice_0.20-35
## [15] ff_2.2-13 Matrix_1.2-8
## [17] foreach_1.4.3 stringr_1.2.0
## [19] affxparser_1.46.0 grid_3.3.3
## [21] BiocParallel_1.8.1 XML_3.98-1.5
## [23] magrittr_1.5 GenomicAlignments_1.10.1
## [25] Rsamtools_1.26.1 codetools_0.2-15
## [27] splines_3.3.3 SummarizedExperiment_1.4.0
## [29] KernSmooth_2.23-15 stringi_1.1.3
## [31] RCurl_1.95-4.8 affyio_1.44.0
作者: ZGUANG@LZU

Created: 2017-03-27 一 18:03
————————————————
版權聲明：本文為CSDN博主「金子哦」的原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/u014801157/article/details/66974577