title: GDAS004-Bioconductor與基因組級(jí)數(shù)據(jù)分析簡(jiǎn)介
date: 2019-09-03 12:0:00
type: "tags"
tags:
- Bioconductor
categories: - Genomics Data Analysis Series
前言
在生命科學(xué)研究領(lǐng)域的相關(guān)人員街州，經(jīng)常使用電子表格來(lái)記錄和管理實(shí)驗(yàn)數(shù)據(jù)伪阶。針對(duì)這個(gè)問題备埃，我們將會(huì)提出一個(gè)使用R和Bioconductor管理基因組級(jí)實(shí)驗(yàn)數(shù)據(jù)的原則暮屡。使用基于R的軟件來(lái)管理數(shù)據(jù)似乎并不太常見肋层。通常來(lái)說(shuō)，R僅用于編寫程序抚恒，構(gòu)建統(tǒng)計(jì)模型拄衰。但是，針對(duì)復(fù)雜生物實(shí)驗(yàn)數(shù)據(jù)的多方面需求這種情況鬼癣，Bioconductor已經(jīng)非常著手以統(tǒng)一的方式來(lái)詳細(xì)地注釋這些數(shù)據(jù)让腹。Bioconductor的這種做法可以讓我們導(dǎo)入初步分析的結(jié)果（例如多個(gè)微陣列的掃描結(jié)果以及多個(gè)測(cè)序結(jié)果的比對(duì)信息），將這些初步的結(jié)果與有關(guān)實(shí)驗(yàn)證的元數(shù)據(jù)扣溺，以及檢測(cè)的樣本信息結(jié)合起來(lái)，從而將所有的信息都保存在一個(gè)R變量中瓜晤。通過周全的設(shè)計(jì)該變量的數(shù)據(jù)結(jié)構(gòu)锥余，程序可以對(duì)這些數(shù)據(jù)進(jìn)行各種操作，最終獲得更加有效的管理和分析環(huán)境痢掠。
在這一節(jié)中驱犹，我們會(huì)通過一些實(shí)驗(yàn)數(shù)據(jù)來(lái)演示一下基本的操作方法，從而足画，逐步過濾到更加全面的表示方法（
我們將通過對(duì)Bioconductor中所表示的實(shí)驗(yàn)數(shù)據(jù)的一些非承劬裕基本的表示法，逐步走向更先進(jìn)的綜合表示法淹辞，從而理解這種數(shù)據(jù)操作方法所表現(xiàn)出來(lái)的高效率医舆。
案例說(shuō)明 (Youtube上有相應(yīng)視頻)
1.0 Distinct components from an experiment.
library(GSE5859Subset)
library(Biobase)
.oldls = ls()
data(GSE5859Subset)
.newls = ls()
我們導(dǎo)入了一個(gè)GSE5859Subset的實(shí)驗(yàn)數(shù)據(jù)，現(xiàn)在我們使用 setdiff()函數(shù)來(lái)比較一下導(dǎo)入前后文件的區(qū)別象缀，如下所示：
newstuff = setdiff(.newls, .oldls)
newstuff
## [1] "geneAnnotation" "geneExpression" "sampleInfo"
從結(jié)果中我們可以發(fā)現(xiàn)蔬将，當(dāng)導(dǎo)入了 GSE5859Subset實(shí)驗(yàn)數(shù)據(jù)中，這個(gè)數(shù)據(jù)中含有3個(gè)文件央星。它們3個(gè)分別屬于不同的類霞怀，如下所示：
cls = lapply(newstuff, function(x) class(get(x)))
names(cls) = newstuff
cls
## $$geneAnnotation
## [1] "data.frame"
##
## $$geneExpression
## [1] "matrix"
##
## $$sampleInfo
## [1] "data.frame"
我們可以按照種族繪制繪制出它們不同時(shí)間的圖形，如下所示：
boxplot(date~factor(ethnicity), data=sampleInfo, ylab="chip date")
假如我們想查看一下BRCA2在不同種族中的表達(dá)情況莉给，我們可以按以下操作進(jìn)行：
- 我們需要知道
geneExpression這個(gè)矩陣中哪一行代表BRCA2的數(shù)據(jù)毙石，并且我們要知道geneAnnotation哪一列提供了BRCA2的信息（探針名與基因名的對(duì)應(yīng)關(guān)系）廉沮； - 我們需要知道
geneExpression這個(gè)矩陣中的順序與samkpleInfo是一致的，要知道sampleInfo中的信息與geneExpression中的標(biāo)簽的對(duì)應(yīng)關(guān)系徐矩。
因此滞时，我為了將基因表達(dá)矩陣與樣本信息聯(lián)系起來(lái)，我們需要對(duì)三個(gè)數(shù)據(jù)進(jìn)行操作（分別為表達(dá)矩陣丧蘸，樣本信息漂洋，基因注釋信息），從而確保這些信息能夠匹配起來(lái)力喷，如下所示：
all.equal(colnames(geneExpression), sampleInfo$filename)
## [1] TRUE
all.equal(rownames(geneExpression), geneAnnotation$PROBEID)
## [1] TRUE
ind = which(geneAnnotation$SYMBOL=="BRCA2")
boxplot(geneExpression[ind,]~sampleInfo$ethnicity, ylab="BRCA2 by hgfocus")
雖然我們繪制出了圖形刽漂，但是這略微復(fù)雜，為了降低操作難度弟孟，Bioconductor定義了一個(gè)名為ExpressionSet數(shù)據(jù)結(jié)構(gòu)贝咙。
ExpressionSet數(shù)據(jù)結(jié)構(gòu)
我們使用簡(jiǎn)單的命令就能構(gòu)建起ExpressionSet這種數(shù)據(jù)結(jié)構(gòu)，如下所示：
es1 = ExpressionSet(geneExpression)
pData(es1) = sampleInfo
fData(es1) = geneAnnotation
es1
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 8793 features, 24 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: 107 122 ... 159 (24 total)
## varLabels: ethnicity date filename group
## varMetadata: labelDescription
## featureData
## featureNames: 1 30 ... 12919 (8793 total)
## fvarLabels: PROBEID CHR CHRLOC SYMBOL
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
上面的es 這個(gè)單一變量就濃縮了整個(gè)實(shí)驗(yàn)的數(shù)據(jù)拂募，將原來(lái)的三個(gè)數(shù)據(jù)變?yōu)榱艘粋€(gè)⊥バ桑現(xiàn)在我們使用es1就能直接繪圖，如下所示：
boxplot(es1$date~es1$ethnicity)
現(xiàn)在我們已經(jīng)了解了ExpressionSet容器的設(shè)計(jì)陈症，你自己也能構(gòu)建一個(gè)ExpressionSet實(shí)例蔼水。我們將會(huì)在下面的實(shí)驗(yàn)中，使用一個(gè)更加廣泛的數(shù)據(jù)表示方法录肯，如下所示：
library(GSE5859)
data(GSE5859)
e
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 8793 features, 208 samples
## element names: exprs, se.exprs
## protocolData
## sampleNames: GSM25349.CEL.gz GSM25350.CEL.gz ...
## GSM136729.CEL.gz (208 total)
## varLabels: ScanDate
## varMetadata: labelDescription
## phenoData
## sampleNames: 1 2 ... 208 (208 total)
## varLabels: ethnicity date filename
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: hgfocus
元數(shù)據(jù)管理(Metadata management)
及時(shí)地了解關(guān)于實(shí)驗(yàn)意圖和方法的信息非常有用趴腋。GSE5859中使用的一些微陣列數(shù)據(jù)就是類來(lái)源于一篇發(fā)表于2004年（PubmedID為15269782）的論文。在聯(lián)網(wǎng)的情況下论咏，我們可以將論文的信息導(dǎo)入R优炬，如下所示：
library(annotate)
p = pmid2MIAME("15269782")
p
## Experiment data
## Experimenter name: Morley M
## Laboratory: NA
## Contact information:
## Title: Genetic analysis of genome-wide variation in human gene expression.
## URL:
## PMIDs: 15269782
##
## Abstract: A 167 word abstract is available. Use 'abstract' method.
現(xiàn)在查看你R環(huán)境中的e變量，使用experimentData(e) <- p后厅贪，就能夠?qū)ubmedID信息添加到e中蠢护。使用experimentData(e)也能發(fā)現(xiàn)e中的信息改變了，也含有Pubmed ID信息养涮。
關(guān)于ExpressionSet的Why與How
現(xiàn)在我們通過創(chuàng)建以及使用ExpressionSet 實(shí)例已經(jīng)快速了解了這種數(shù)據(jù)結(jié)構(gòu)的特點(diǎn)葵硕。我們現(xiàn)在將要了解一下，為什么要這樣設(shè)計(jì)這種數(shù)據(jù)結(jié)構(gòu)单寂，以及如何使用這些數(shù)據(jù)結(jié)構(gòu)贬芥。
設(shè)計(jì)ExpressionSet 這種數(shù)據(jù)結(jié)構(gòu)的初衷有很多，以下幾點(diǎn)最為關(guān)鍵：
- 降低復(fù)雜分析宣决，注釋以及樣本信息編程的難度蘸劈；
- 支持特征值（其實(shí)就是表達(dá)值）和樣本信息的簡(jiǎn)單的過濾；
- 在將分析數(shù)據(jù)傳遞給分析函數(shù)時(shí)控制內(nèi)存占用尊沸。
以上就是為什么要設(shè)計(jì)ExpressionSet這種數(shù)據(jù)結(jié)構(gòu)威沫。關(guān)于如何使用這種結(jié)構(gòu)贤惯，這些內(nèi)容就要涉及一些設(shè)計(jì)模式問題，在R中棒掠，與面向?qū)ο缶幊滔嚓P(guān)類被稱為S4類孵构。
面向?qū)ο缶幊膛cS4類
面向?qū)ο缶幊?OOP, Object-oriented programming)是計(jì)算機(jī)編程一個(gè)學(xué)科，它有助于降低程序的復(fù)雜性和維護(hù)要求烟很，并能增強(qiáng)程度的互動(dòng)性颈墅。OOP有著各種解決方式，即使在R中雾袱，也可以進(jìn)行OOP恤筛，這就是所謂的S4方法。在R中芹橡，OOP的基本思想就是一個(gè)對(duì)象有些特征：
- 將不同類型的數(shù)據(jù)集中放到一個(gè)統(tǒng)一的結(jié)構(gòu)中毒坛，
- 它被標(biāo)記為某個(gè)正式指定類的成員，并且可以正式地位于相關(guān)對(duì)象的系統(tǒng)中林说，
- 以類特異性方式對(duì)泛型方法進(jìn)行響應(yīng)煎殷。
我們可以使用getClass函數(shù)查看一個(gè)ExpresionSet的類結(jié)構(gòu)，如下所示：
getClass("ExpressionSet")
## Class "ExpressionSet" [package "Biobase"]
##
## Slots:
##
## Name: experimentData assayData phenoData
## Class: MIAME AssayData AnnotatedDataFrame
##
## Name: featureData annotation protocolData
## Class: AnnotatedDataFrame character AnnotatedDataFrame
##
## Name: .__classVersion__
## Class: Versions
##
## Extends:
## Class "eSet", directly
## Class "VersionedBiobase", by class "eSet", distance 2
## Class "Versioned", by class "eSet", distance 3
從上面我們可以看到腿箩，一共有7個(gè)組件定義了該類豪直。只有annoataion這個(gè)組件占據(jù)了一個(gè)單獨(dú)的R類型，也就是character類型珠移。至于其它的組件顶伞，都是S4類的一個(gè)實(shí)例（或者說(shuō)是更多的R數(shù)據(jù)類型）。
舉例說(shuō)明一下剑梳，MIAME類的定義如下所示：
getClass("MIAME")
## Class "MIAME" [package "Biobase"]
##
## Slots:
##
## Name: name lab contact
## Class: character character character
##
## Name: title abstract url
## Class: character character character
##
## Name: pubMedIds samples hybridizations
## Class: character list list
##
## Name: normControls preprocessing other
## Class: list list list
##
## Name: .__classVersion__
## Class: Versions
##
## Extends:
## Class "MIAxE", directly
## Class "characterORMIAME", directly
## Class "Versioned", by class "MIAxE", distance 2
因此，我們可以發(fā)現(xiàn)滑潘，S4可以讓我們定義包含各種信息的新數(shù)據(jù)結(jié)垢乙。我們雖然可以使用R中非常通用的列表結(jié)構(gòu)來(lái)實(shí)現(xiàn)這一目的（就是說(shuō)實(shí)現(xiàn)封裝不同信息），但是S4允許我們能保證有新結(jié)構(gòu)的屬性语卤，以便為這些結(jié)構(gòu)定義的程序不必“檢查”它們正在操作什么追逮。這些檢查過程內(nèi)置在了類系統(tǒng)中。
數(shù)據(jù)框(DataFrame): 表格數(shù)據(jù)的增強(qiáng)型結(jié)構(gòu)
數(shù)據(jù)框是R中非炒舛妫基礎(chǔ)的數(shù)據(jù)結(jié)構(gòu)钮孵。Bioconductor中引入了DataFrame類，它可以實(shí)現(xiàn)其它一些額外的功能眼滤，后面我們會(huì)提到巴席。
library(S4Vectors)
sa2 = DataFrame(sampleInfo)
sa2
## DataFrame with 24 rows and 4 columns
## ethnicity date filename group
## <factor> <Date> <character> <numeric>
## 1 ASN 2005-06-23 GSM136508.CEL.gz 1
## 2 ASN 2005-06-27 GSM136530.CEL.gz 1
## 3 ASN 2005-06-27 GSM136517.CEL.gz 1
## 4 ASN 2005-10-28 GSM136576.CEL.gz 1
## 5 ASN 2005-10-07 GSM136566.CEL.gz 1
## ... ... ... ... ...
## 20 ASN 2005-06-27 GSM136524.CEL.gz 0
## 21 ASN 2005-06-23 GSM136514.CEL.gz 0
## 22 ASN 2005-10-07 GSM136563.CEL.gz 0
## 23 ASN 2005-10-07 GSM136564.CEL.gz 0
## 24 ASN 2005-10-07 GSM136572.CEL.gz 0
show方法描述頂部和底部的行數(shù)據(jù)，提供有關(guān)列數(shù)據(jù)類型的信息诅需，并且能夠比data.frame數(shù)據(jù)結(jié)構(gòu)更有效地管理列中的復(fù)雜數(shù)據(jù)類所漾唉，下面是類定義：
getClass("DataFrame")
## Class "DataFrame" [package "S4Vectors"]
##
## Slots:
##
## Name: rownames nrows listData
## Class: character_OR_NULL integer list
##
## Name: elementType elementMetadata metadata
## Class: character DataTable_OR_NULL list
##
## Extends:
## Class "DataTable", directly
## Class "SimpleList", directly
## Class "DataTable_OR_NULL", by class "DataTable", distance 2
## Class "List", by class "SimpleList", distance 2
## Class "Vector", by class "SimpleList", distance 3
## Class "Annotated", by class "SimpleList", distance 4
源于多樣本NGS數(shù)據(jù)的成熟數(shù)據(jù)表示方法
短讀長(zhǎng)測(cè)序產(chǎn)生了遠(yuǎn)比微陣列多得多的數(shù)據(jù)荧库。測(cè)序?qū)嶒?yàn)也可以通過更加多樣化的手段進(jìn)行分析，因此測(cè)序?qū)嶒?yàn)?zāi)軌蛟诓煌６燃?jí)別(at various levels of granularity)提供有效的輸入赵刑。
處理RNA-seq數(shù)據(jù)的典型流程包括：
- 展示每個(gè)周期中堿基所調(diào)用的色譜（2008年后就不再提供了）分衫；
- 堿基的讀取質(zhì)量值，它們通常以FASTQ格式進(jìn)行儲(chǔ)存般此；
- 通過基因組或轉(zhuǎn)錄本比對(duì)后蚪战，對(duì)各個(gè)基因統(tǒng)計(jì)的結(jié)果；
- 基因組極的多個(gè)樣本的數(shù)據(jù)儲(chǔ)存铐懊，用于推斷感興趣的生物學(xué)問題邀桑。通常一個(gè)樣本對(duì)應(yīng)兩個(gè)FASTQ文件。
Bioconductor中的BAM文件
RNAseqData.HNRNPC.bam.chr14 中含有一些BAM文件居扒，這些BAM文件是8個(gè)Hela細(xì)胞系樣本的雙端測(cè)序數(shù)據(jù)概漱，其中4個(gè)樣本通過RNAi的方式敲低了異質(zhì)核核糖核蛋白C1和C2（hnRNPC,heterogeneous nuclear ribonucleoproteins C1 and C2)。這個(gè)家庭的核蛋白參與pre-mRNA的加工∠参梗現(xiàn)在我們演示一下瓤摧，如何從BAM文件中來(lái)統(tǒng)計(jì)與HNRNPC有關(guān)的讀長(zhǎng)數(shù)目(reads count)。
BAM文件儲(chǔ)存在RNAseqData.HNRNPC.bam.chr14 包中玉吁，如下所示：
library(RNAseqData.HNRNPC.bam.chr14)
bfp = RNAseqData.HNRNPC.bam.chr14_BAMFILES
length(bfp)
## [1] 8
bfp[1]
## ERR127306
## "/Library/Frameworks/R.framework/Versions/3.4/Resources/library/RNAseqData.HNRNPC.bam.chr14/extdata/ERR127306_chr14.bam"
我們使用Rsamtools包以一個(gè)變量的形式管理這些信息照弥，如下所示：
library(Rsamtools)
bfl = BamFileList(file=bfp)
bfl
## BamFileList of length 8
## names(8): ERR127306 ERR127307 ERR127308 ... ERR127303 ERR127304 ERR127305
BAM文件包括關(guān)于讀長(zhǎng)(reads)與基因組比對(duì)后的有限數(shù)量的元數(shù)據(jù)。這里我們以BamFileList實(shí)例來(lái)說(shuō)明一下：
seqinfo(bfl)
## Seqinfo object with 93 sequences from an unspecified genome:
## seqnames seqlengths isCircular genome
## chr1 249250621 <NA> <NA>
## chr10 135534747 <NA> <NA>
## chr11 135006516 <NA> <NA>
## chr11_gl000202_random 40103 <NA> <NA>
## chr12 133851895 <NA> <NA>
## ... ... ... ...
## chrUn_gl000247 36422 <NA> <NA>
## chrUn_gl000248 39786 <NA> <NA>
## chrUn_gl000249 38502 <NA> <NA>
## chrX 155270560 <NA> <NA>
## chrY 59373566 <NA> <NA>
我們看到基因”沒有指定“进副，但是我們碰巧知道參考基因組hg19用于進(jìn)入基因序列的比對(duì)这揣。使用[genome(bfl) <- “hg19”可以糾正這點(diǎn)，但是還存在著一個(gè)bug影斑。這個(gè)bug應(yīng)該修正一下给赞，以便可以標(biāo)記使用不同的參考基因組來(lái)進(jìn)行比對(duì)或糾正錯(cuò)誤，注：這段不太理解矫户，原文如下：
We see that the genome is “unspecified”, but we happen to know that the reference build ‘hg19’ was used to define the genomic sequence for alignment. [genome(bfl) <- “hg19” should rectify this but at present there is a bug. This should be fixed so that attempts to compare alignments using different
references can be flagged or throw errors.]
為了驗(yàn)證是否成功敲除片迅，我們將檢查編碼HNRNPC的基因組附近的序列的比對(duì)情況。稍后皆辽，我們將展示如何使用Bioconductor注釋工具獲取基因地址柑蛇；不過這個(gè)地址僅用于hg19基因組，如下所示：
hnrnpcLoc = GRanges("chr14", IRanges(21677296, 21737638))
hnrnpcLoc
## GRanges object with 1 range and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr14 [21677296, 21737638] *
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
這標(biāo)識(shí)了HNRNPC編碼基因?yàn)閔g19注釋的染色體14的區(qū)域驱闷。評(píng)估表達(dá)式量化的一種簡(jiǎn)單方法是使用GenomicAlignmentspackage中的summary izeOverlaps耻台。此計(jì)數(shù)方法將隱式使用多核計(jì)算。
上面的代碼標(biāo)識(shí)了14號(hào)染色體中編碼HNRNPC附近的序列空另。評(píng)估基因表達(dá)量化的一種簡(jiǎn)單方法就是使用GenomicAlignmentspackage包中的summarizeOverlaps 函數(shù)盆耽，這個(gè)計(jì)算方法隱式地使用了多核計(jì)算，如下所示：
library(GenomicAlignments)
library(BiocParallel)
register(SerialParam())
hnse = summarizeOverlaps(hnrnpcLoc,bfl)
練習(xí)
- 檢查
hnse中的count組件內(nèi)容。與HNRNPC敲除相對(duì)應(yīng)的有關(guān)量化的指標(biāo)是什么征字？ - 向
hnse的colData組件中添加一個(gè)變量condition都弹，用于指示敲除株與野生株。 - 具有后綴306和307的樣本是單個(gè)對(duì)照的技術(shù)重復(fù)匙姜；308和309同樣如此畅厢。因此使用了6個(gè)不同的生物學(xué)樣本。現(xiàn)在向
colData中添加一個(gè)變量氮昧，用于區(qū)別這6個(gè)樣本框杜。
ExpressionSet設(shè)計(jì)的一些復(fù)雜性(進(jìn)階)
生物信息學(xué)家是在一個(gè)不斷變化的環(huán)境中工作。自Bioconductor建立以來(lái)袖肥，計(jì)算機(jī)的容量和通量已經(jīng)大幅度增加咪辱。R語(yǔ)言也性現(xiàn)了很大的變量。在早期的時(shí)候椎组，大家都比較注意內(nèi)存的使用油狂，因?yàn)楫?dāng)時(shí)在R中，向一個(gè)函數(shù)傳遞參數(shù)時(shí)寸癌，會(huì)復(fù)制內(nèi)存中的對(duì)象专筷， 這種操作比較低效。而“環(huán)境”類中的R對(duì)象則不會(huì)發(fā)生復(fù)制蒸苇，這就導(dǎo)致了使用環(huán)境來(lái)儲(chǔ)存大量的數(shù)據(jù)用于分析磷蛹。
class(slot(es1, "assayData"))
## [1] "environment"
ae = slot(es1, "assayData")
ae
## <environment: 0x7f813b4dc698>
“環(huán)境”是R中一種非常特殊的對(duì)象類型，可以將它視為對(duì)表達(dá)值矩陣引用的容器溪烤。
ls(ae)
## [1] "exprs"
dim(ae$exprs)
## [1] 8793 24
dim(exprs(es1))
## [1] 8793 24
object.size(ae)
## 56 bytes
object.size(ae$exprs)
## 2253824 bytes
練習(xí)（可選）
請(qǐng)解釋以下運(yùn)行結(jié)果：
e1 = new.env()
l1 = list()
e1$x = 5
l1$x = 5
all.equal(e1$x, l1$x)
## [1] TRUE
f = function(z) {z$x = 4; z$x}
all.equal(f(e1), f(l1))
## [1] TRUE
all.equal(e1$x, l1$x) # different from before
## [1] "Mean relative difference: 0.25"
運(yùn)行結(jié)果如下所示：
> e1 = new.env()
> l1 = list()
> e1$x = 5
> l1$x = 5
> all.equal(e1$x, l1$x)
[1] TRUE
> ## [1] TRUE
> f = function(z) {z$x = 4; z$x}
> all.equal(f(e1), f(l1))
[1] TRUE
> ## [1] TRUE
> all.equal(e1$x, l1$x) # different from before
[1] "Mean relative difference: 0.25"
> ## [1] "Mean relative difference: 0.25"
