title: GDAS006-管理與處理大量BED文件
date: 2019-09-05 12:0:00
type: "tags"
tags:

• Bioconductor
  categories:
• Genomics Data Analysis Series

前言

我們已經(jīng)通過ExpressionSet類和BamFileList文件（含有BAM比對后的文件）來了解了一些關(guān)于數(shù)據(jù)管理的基本知識枫匾。

• S4類有助于控制某個數(shù)據(jù)結(jié)構(gòu)中不同組件之間的復(fù)雜程度和一致性。
• R包能夠進(jìn)一步用于協(xié)調(diào)和發(fā)布數(shù)據(jù)精置、軟件程度和文檔睁本。在這個課程的后面部分里我們會學(xué)習(xí)如何構(gòu)建R包主卫。

在這個一部分的結(jié)尾處局骤，我們來解決一個明顯沒有最佳解決方案的問題。我們考慮打包GRanges的序列化后的實例以及已經(jīng)建立了索引的壓縮文本文件讲竿，這些文件都可以使用GenomicFiles進(jìn)行并行計算處理泥兰。

為什么這是一個比較難解決的問題呢？因為在某些情況下题禀，含有大量功能的GRanges實例將會非常大鞋诗，加載它們可能會花費大量的時間，并且會占用大量的內(nèi)存迈嘹。這可能導(dǎo)致在RAM有限時削彬，在群集節(jié)點上使用這些表示方法出現(xiàn)困難。如果數(shù)據(jù)以允許隨機訪問的索引格式保存在磁盤上秀仲，則根據(jù)需要將文件內(nèi)容的目標(biāo)提取融痛，并轉(zhuǎn)換為GRanges，就可能為某些應(yīng)用程度帶來更好的性能神僵。創(chuàng)建erma包的一個原因就是獲取關(guān)于這些不同方法之間平衡后的數(shù)據(jù)雁刷。

erma：表觀遺傳路線圖簡介

為了促進(jìn)對DNA突變體的理解，我們對細(xì)胞類型之間的表觀基因組變異有所興趣保礼。NIH表觀基因組路線圖收集了很多關(guān)于ChIP-seq實驗分析的數(shù)據(jù)沛励，在這個公眾數(shù)據(jù)庫里記錄了數(shù)百個不同細(xì)胞類型和標(biāo)志特異性文件供其他研究者使用。

管理這類數(shù)據(jù)的有效方法是什么呢炮障？其中erma 包就提供了這樣一種方法目派。

BED文件集合

如同Ernst和Kellis使用大量細(xì)胞系的數(shù)據(jù)所報道的那樣，他們已經(jīng)進(jìn)行了染色質(zhì)狀態(tài)的統(tǒng)計“估算”胁赢，如下所示：

library(erma)
ef = dir(system.file("bed_tabix", package="erma"), patt="bed.gz$")
length(ef)
## [1] 31
head(ef)
## [1] "E002_25_imputed12marks_mnemonics.bed.gz"
## [2] "E003_25_imputed12marks_mnemonics.bed.gz"
## [3] "E021_25_imputed12marks_mnemonics.bed.gz"
## [4] "E032_25_imputed12marks_mnemonics.bed.gz"
## [5] "E033_25_imputed12marks_mnemonics.bed.gz"
## [6] "E034_25_imputed12marks_mnemonics.bed.gz"

使用GenomicFiles進(jìn)行管理

我們認(rèn)為文件名是不變的企蹭，但是儲存BED文件路徑的向量則可以被命名得更詳細(xì)。通過makeErmaSet函數(shù)可以獲取元數(shù)據(jù)，如下所示：

mm = makeErmaSet()
## NOTE: input data had non-ASCII characters replaced by ' '.
mm
## ErmaSet object with 0 ranges and 31 files: 
## files: E002_25_imputed12marks_mnemonics.bed.gz, E003_25_imputed12marks_mnemonics.bed.gz, ..., E088_25_imputed12marks_mnemonics.bed.gz, E096_25_imputed12marks_mnemonics.bed.gz 
## detail: use files(), rowRanges(), colData(), ...
head(colData(mm))
## (showing narrow slice of  6 x 95  DataFrame)   narrDF with 6 rows and 6 columns
##      Epigenome.ID..EID.          GROUP Epigenome.Mnemonic     ANATOMY
##             <character>    <character>        <character> <character>
## E002               E002            ESC            ESC.WA7         ESC
## E003               E003            ESC             ESC.H1         ESC
## E021               E021           iPSC        IPSC.DF.6.9        IPSC
## E032               E032   HSC & B-cell       BLD.CD19.PPC       BLOOD
## E033               E033 Blood & T-cell        BLD.CD3.CPC       BLOOD
## E034               E034 Blood & T-cell        BLD.CD3.PPC       BLOOD
##                TYPE         SEX
##         <character> <character>
## E002 PrimaryCulture      Female
## E003 PrimaryCulture        Male
## E021 PrimaryCulture        Male
## E032    PrimaryCell        Male
## E033    PrimaryCell     Unknown
## E034    PrimaryCell        Male
## use colnames() for full set of metadata attributes.
names(files(mm)) = colData(mm)$Epigenome.Mnemonic
head(files(mm))
##                                                                                                                 ESC.WA7 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E002_25_imputed12marks_mnemonics.bed.gz" 
##                                                                                                                  ESC.H1 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E003_25_imputed12marks_mnemonics.bed.gz" 
##                                                                                                             IPSC.DF.6.9 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E021_25_imputed12marks_mnemonics.bed.gz" 
##                                                                                                            BLD.CD19.PPC 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E032_25_imputed12marks_mnemonics.bed.gz" 
##                                                                                                             BLD.CD3.CPC 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E033_25_imputed12marks_mnemonics.bed.gz" 
##                                                                                                             BLD.CD3.PPC 
## "/Library/Frameworks/R.framework/Versions/3.3/Resources/library/erma/bed_tabix/E034_25_imputed12marks_mnemonics.bed.gz"

使用表格描述染色質(zhì)狀態(tài)信息

定義染色體狀態(tài)與標(biāo)記概述

下表顯示了不同類型的染色質(zhì)標(biāo)記是如何組合并用于判斷染色體質(zhì)片段的狀態(tài)的谅摄。那些縮寫的意義為：活躍TSS(TssA)徒河，兩種類型的上下流啟動子(PromU, PromD1, PromD2)，轉(zhuǎn)錄為5'或3‘優(yōu)先或略強或弱(Tx5’, Tx3’, Tx, TxWk)螟凭，聯(lián)合轉(zhuǎn)錄或調(diào)控(TxReg)虚青，優(yōu)先轉(zhuǎn)錄，增強子活性和可能弱(TxEnh5‘螺男，TxEnh3‘棒厘，TxEnhW)，兩種類型的活躍增強子下隧，可能是側(cè)翼或弱的(EnhA1奢人，EnhA2，EnhAF淆院，EnhW1何乎，EnhW2)或由H3K27ac(EnhAc)，初級DNase(DNase)土辩，ZNF基因和重復(fù)序列(ZNF/rpts)支救，異染色質(zhì)(Het)，二價或平衡型啟動子(Propp拷淘，PromBiv)各墨，抑制多梳(ReprPC)或靜默(Quies)。

library(png)
im = readPNG(system.file("pngs/emparms.png", package="erma"))
grid.raster(im)

plot of chunk lkxy

erma包中含有的文件包括31個不同細(xì)胞系和基因組的信息（原文里面是tilings启涯，這個我的理解就是表示了一些信息）贬堵，并根據(jù)上面顯示了狀態(tài)圖將這些信息分配給染色體狀態(tài)。

下面展示了兩個基因BRCA2和Eome啟動子區(qū)域的狀態(tài)標(biāo)記结洼，如下所示：

stateProfile(mm, "BRCA2")
## 'select()' returned 1:many mapping between keys and columns
## Warning: executing %dopar% sequentially: no parallel backend registered
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

[圖片上傳失敗...(image-7b3885-1573126148074)]

stateProfile(mm, "EOMES")
## 'select()' returned 1:many mapping between keys and columns
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

plot of chunk lkbr2

有兩個基因在啟動子區(qū)域的狀態(tài)與常規(guī)的狀態(tài)不同黎做，以及與細(xì)胞類型的一致性方面也不同。

并行目標(biāo)查詢

雖然可視化很直觀松忍，但是我們也希望能夠進(jìn)行文本編輯蒸殿，從而提取一些信息供下游使用。在這里我們簡單地將每個細(xì)胞類型的狀態(tài)制成表格鸣峭。使用readuceByFile函數(shù)就能將文件分配給可用的內(nèi)核或主機進(jìn)行并行處理伟桅，如下所示：

gm = promoters(range(genemodel("BRCA2"))) # 2000 upstream, 200 down by default
library(BiocParallel)
register(MulticoreParam(workers=2))
library(GenomicFiles)
ans = reduceByFile( gm, files(mm), MAP=function(range,file) {
   table( import(file, genome="hg19", which=range)$name ) } )
ans= unlist(ans, recursive=FALSE)
names(ans) = colData(mm)$Epigenome.Mnemonic
ans[1:4]
## $$ESC.WA7
## 
##   1_TssA 16_EnhW1  2_PromU 25_Quies 
##        1        1        1        1 
## 
## $$ESC.H1
## 
##   1_TssA 16_EnhW1 19_DNase  2_PromU 25_Quies 3_PromD1 
##        1        1        1        1        1        1 
## 
## $$IPSC.DF.6.9
## 
##   1_TssA 16_EnhW1 19_DNase  2_PromU 25_Quies 
##        1        1        1        1        1 
## 
## $$BLD.CD19.PPC
## 
##   1_TssA 16_EnhW1 19_DNase  2_PromU 25_Quies 3_PromD1 
##        1        1        1        1        1        1

因此，我們看到叽掘，在某個特定的基因的啟動子區(qū)域的基因序列方面，不同的細(xì)胞類型存在著數(shù)量和狀態(tài)類型的變化玖雁。在后面的練習(xí)中更扁，我們將會考慮如何以這種方式開發(fā)其它潛在感興趣的信息。

總結(jié)

在這一部分中我們看到：

• 在磁盤中，使用了大量添加了索引的文本來表示基因組注釋或?qū)嶒灲Y(jié)果（原文是：“l(fā)eave on disk” large numbers of indexed textual representations of genomic annotation or experimental result）浓镜；
• 使用R包的方式來簡化標(biāo)記這些數(shù)據(jù)溃列；
• 使用GenomicFiles來管理有關(guān)文件的元數(shù)據(jù)；
• 使用reduceByFile來并行處理文件內(nèi)容膛薛。

這些技術(shù)避免了創(chuàng)建听隐，加載大量的R對象，從而提供了某種程度的便捷性哄啄。在內(nèi)存資源有限的情況下雅任，對于那些需要大量集群(clusters)的應(yīng)用來說，這種策略就顯得性價比較高咨跌。

參考資料

1. Managing and processing large numbers of BED files