title: GDAS006-深入理解SummarizedExperiment類(lèi)
date: 2019-09-05 12:0:00
type: "tags"
tags:
- Bioconductor
categories: - Genomics Data Analysis Series
前言
在數(shù)據(jù)管理部分，我們應(yīng)用summarizeOverlaps 來(lái)管理RNAseqData.HNRNPC.bam.chr14中的BAM文件≈匣冢現(xiàn)在我們?cè)俅蝸?lái)操作一下瓢对。
使用SummarizedExperiment來(lái)管理BAM文件
library(RNAseqData.HNRNPC.bam.chr14)
bfp = RNAseqData.HNRNPC.bam.chr14_BAMFILES
library(Rsamtools)
bfl = BamFileList(file=bfp)
hnrnpcLoc = GRanges("chr14", IRanges(21677296, 21737638))
library(GenomicAlignments)
library(BiocParallel)
register(SerialParam())
hnse = summarizeOverlaps(hnrnpcLoc,bfl)
hnse
## class: RangedSummarizedExperiment
## dim: 1 8
## metadata(0):
## assays(1): counts
## rownames: NULL
## rowData names(0):
## colnames(8): ERR127306 ERR127307 ... ERR127304 ERR127305
## colData names(0):
hnse 是RangedSummarizedExperiment 類(lèi)的一實(shí)例照捡。這個(gè)類(lèi)就類(lèi)似于 ExpressionSet ，不過(guò)有著更多的內(nèi)容來(lái)管理元數(shù)據(jù)叠萍，它的流程如下所示：
有效地使用SummarizedExperiment實(shí)例涉及學(xué)習(xí)它的一些方法负敏。為了獲取HNRNPC基因的讀長(zhǎng)/區(qū)域(read/region)霞掺，我們可以使用assay方法砸紊，如下所示：
assay(hnse)
## ERR127306 ERR127307 ERR127308 ERR127309 ERR127302 ERR127303 ERR127304
## [1,] 5422 6320 5896 5558 172 196 316
## ERR127305
## [1,] 282
以上就是最基本的結(jié)果表示方法传于。列名則是樣本標(biāo)識(shí)符，但是有關(guān)區(qū)域檢查的一些信息則已經(jīng)丟失批糟。
SummarizedExperiment中的元數(shù)據(jù)
hnse 對(duì)象還有一些其它的信息格了，如下所示：
rowRanges(hnse)
## 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
seqinfo(hnse)
## Seqinfo object with 1 sequence from an unspecified genome; no seqlengths:
## seqnames seqlengths isCircular genome
## chr14 NA NA <NA>
metadata(hnse)
## list()
我們還可以進(jìn)一步分析差異信息，從而輸出更多徽鼎，更廣泛的信息盛末。
通過(guò)添加輸入信息來(lái)更有效地生成SummarizedExperiment
使用元數(shù)據(jù)定義感興趣的區(qū)域
We have seen that it is sufficient to define a single GRanges to drive summarizeOverlaps over a set of BAM files. We’d like to preserve more metadata about the regions examined. We’ll use the TxDb infrastructure, to be described in more detail later, to get a structure defining gene regions on chr14. We’ll also use the Homo.sapiens annotation package to add gene symbols.
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
gr14 = genes(txdb, vals=list(tx_chrom="chr14"))
## Warning: 'elementLengths' is deprecated.
## Use 'elementNROWS' instead.
## See help("Deprecated")
## Warning: 'elementLengths' is deprecated.
## Use 'elementNROWS' instead.
## See help("Deprecated")
gr14$symbol = mapIds(Homo.sapiens, keys=gr14$gene_id, keytype="ENTREZID",
column="SYMBOL")
## 'select()' returned 1:1 mapping between keys and columns
gr14
## GRanges object with 781 ranges and 2 metadata columns:
## seqnames ranges strand | gene_id
## <Rle> <IRanges> <Rle> | <character>
## 10001 chr14 [ 71050957, 71067384] - | 10001
## 100113389 chr14 [ 45580078, 45580176] + | 100113389
## 100113391 chr14 [ 20794600, 20794698] - | 100113391
## 100124539 chr14 [ 91592770, 91592896] + | 100124539
## 100126297 chr14 [101507700, 101507781] + | 100126297
## ... ... ... ... . ...
## 9870 chr14 [ 75127955, 75179807] - | 9870
## 9878 chr14 [ 21945335, 21967319] + | 9878
## 9895 chr14 [102829300, 102968818] + | 9895
## 9950 chr14 [ 93260650, 93306304] + | 9950
## 9985 chr14 [ 24641234, 24649463] + | 9985
## symbol
## <character>
## 10001 MED6
## 100113389 SNORD127
## 100113391 SNORD126
## 100124539 SNORA11B
## 100126297 MIR300
## ... ...
## 9870 AREL1
## 9878 TOX4
## 9895 TECPR2
## 9950 GOLGA5
## 9985 REC8
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genome
定義BAM文件的樣本信息
現(xiàn)在我們有三個(gè)樣本，一個(gè)用于控制否淤，兩個(gè)用于敲低悄但。我們使用GenomicFiles來(lái)綁定樣本信息的元數(shù)據(jù)，如下所示：
char = rep(c("hela_wt", "hela_hkd"), each=4)
bff = GenomicFiles(files=path(bfl))
colData(bff)$condition = char
sid = c(1,1,1,1,2,2,3,3)
bff$sample = sid
bff
## GenomicFiles object with 0 ranges and 8 files:
## files: ERR127306_chr14.bam, ERR127307_chr14.bam, ..., ERR127304_chr14.bam, ERR127305_chr14.bam
## detail: use files(), rowRanges(), colData(), ...
比較讀長(zhǎng)覆蓋區(qū)石抡，保留元數(shù)據(jù)
我們來(lái)查看5個(gè)基因檐嚣，其中就包括HNRNPC。當(dāng)計(jì)算后啰扛，我們會(huì)將樣本信息再綁定回結(jié)果嚎京，如下所示：
hnse = summarizeOverlaps(gr14[c(1:4,305)],files(bff))
colData(hnse) = cbind(colData(hnse), colData(bff))
hnse
## class: RangedSummarizedExperiment
## dim: 5 8
## metadata(0):
## assays(1): counts
## rownames(5): 10001 100113389 100113391 100124539 3183
## rowData names(2): gene_id symbol
## colnames(8): ERR127306 ERR127307 ... ERR127304 ERR127305
## colData names(2): condition sample
assay(hnse)
## ERR127306 ERR127307 ERR127308 ERR127309 ERR127302 ERR127303
## 10001 114 156 129 144 175 213
## 100113389 0 0 0 0 0 0
## 100113391 0 0 0 0 0 0
## 100124539 0 0 0 0 0 0
## 3183 2711 3160 2948 2779 86 98
## ERR127304 ERR127305
## 10001 210 165
## 100113389 0 0
## 100113391 0 0
## 100124539 0 0
## 3183 158 141
從上面我們可以看出，行標(biāo)識(shí)符是隨計(jì)數(shù)矩陣一起出現(xiàn)的隐解，現(xiàn)在我們繪制出這幾個(gè)基因的箱線圖鞍帝，如下所示：
par(mfrow=c(2,2))
for (i in 2:5) {
boxplot(assay(hnse)[i,]~hnse$condition, ylab=rowRanges(hnse)$symbol[i])
}
從ExpressionSet提取數(shù)據(jù)
從ExpressionSet中提取數(shù)據(jù)很容易，但是還需要基于基因組范圍查詢(xún)的陣列探針的子集煞茫，如下所示：
library(ALL)
data(ALL)
allse = makeSummarizedExperimentFromExpressionSet(ALL)
allse
## class: RangedSummarizedExperiment
## dim: 12625 128
## metadata(3): experimentData annotation protocolData
## assays(1): exprs
## rownames(12625): 1000_at 1001_at ... AFFX-YEL021w/URA3_at
## AFFX-YEL024w/RIP1_at
## rowData names(0):
## colnames(128): 01005 01010 ... 83001 LAL4
## colData names(21): cod diagnosis ... f.u date.last.seen
rowRanges(allse)
## Warning: 'elementLengths' is deprecated.
## Use 'elementNROWS' instead.
## See help("Deprecated")
## GRangesList object of length 12625:
## $$1000_at
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
##
## $$1001_at
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
##
## $$1002_f_at
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
##
## ...
## <12622 more elements>
## -------
## seqinfo: no sequences
總結(jié)
RangedSummarizedExperiment類(lèi)實(shí)例化了一些Bioconductor數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)的一些關(guān)鍵原則：
- 基于樣本特征的數(shù)據(jù)分析和元數(shù)據(jù)分析帕涌，并將它們以某種方式綁定在一起。
- 類(lèi)似于矩陣的子集直接用于分析和樣本數(shù)據(jù)分析续徽。
- 基于范圍的子集設(shè)置適用于可通過(guò)基因組坐標(biāo)尋址的分析蚓曼。
- 可以在
mcol中提供關(guān)鍵分析特征的任意元數(shù)據(jù)(rowRanges(se))。 - 可以通過(guò)
metadata(se)<-來(lái)添加任何元數(shù)據(jù)钦扭。
在后面的內(nèi)容里我們還將會(huì)了解更多的關(guān)于SummarizedExperiment改造的一些數(shù)據(jù)結(jié)構(gòu)纫版，從而用于專(zhuān)門(mén)用于RNA-seq的多個(gè)階段處理和數(shù)據(jù)分析。
