1

たくさんの読み取りを含む BAM ファイルがあります。scanBamfromで R にロードできますRsamtools。

ただし、読み取りのサブセットのみが必要です。character興味のあるqnamesを持つベクトルがあります。

scanBam数千回の読み取りすべてのデータを含む 13 要素のリストである 1 要素のリストを返します。

qname構造を維持してこのオブジェクトをサブセット化するにはどうすればよいですか? マニュアルやオンラインで何も見つけることができませんでした。

4

2 に答える 2

2

引数に param=ScanBamParam(what="qname") を指定して qname を含めて GenomicAlignments::readGAlignments でデータを入力した方が便利でしょう。次に、%in% でサブセット化できます。以下は、 ExperimentDataパッケージの 1 つを使用した、より完全な例です。

library(GenomicAlignments)
library(RNAseqData.HNRNPC.bam.chr14)

fname <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]    
want <- c("ERR127306.11930623", "ERR127306.24720935",
    "ERR127306.23706011", "ERR127306.22418829", "ERR127306.13372247",
    "ERR127306.20686334", "ERR127306.11412145", "ERR127306.4711647",
    "ERR127306.7479582", "ERR127306.12737243")
aln <- readGAlignments(fname, param=ScanBamParam(what="qname"))
aln[mcols(aln)$qname %in% want]

もちろん、BAM ファイルは大きく、qname はその大部分を占めます。多くの場合、ファイルをチャンクで反復処理することは理にかなっています。これは、（現在の Rsamtools では）yieldReduce で有効化されます。この場合、妥当な（たとえば 1M）読み取り数に設定された yieldSize を持つ BamFile と、データのチャンクを入力して処理するための MAP 関数（たとえば、不要な読み取りのフィルタリング）が提供されます。 )、結果を連結するための (オプションの) REDUCE 関数、および反復がいつ完了したかを示すための (オプションの) DONE 関数。ソリューションは次のようになります (yieldSize は、サンプル データで説明できるように人為的に小さくなっています)。

bfl <- BamFile(fname, yieldSize=100000)  ## larger, e.g., 1M-5M
MAP <- function(bfl, want) {
    ## message("tick")
    aln <- readGAlignments(bfl, param=ScanBamParam(what="qname"))
    if (length(aln) == 0)
        NA                          # semaphore -- DONE
    else
        aln[mcols(aln)$qname %in% want]
}
REDUCE <- c
DONE <- function(x) identical(x, NA)
result <- yieldReduce(bfl, MAP, REDUCE, DONE, want=want)

scanBam を使用して同様のアプローチを採用することもできますが、データ構造 (リストのリスト) を扱うのはより複雑です。

x <- scanBam(fname, param=ScanBamParam(what=c("qname", "pos")))
keep <- lapply(lapply(x, "[[", "qname"), "%in%", want)
result <- Map(function(elts, keep) {
    lapply(elts, "[", keep)
}, x, keep)

これは、yieldReduce でも使用できます。

フィルタリングされた読み取りで新しいbamファイルを作成することに興味がある場合は、

filter_factory <- function(want) {
    list(KeepQname = function(x) x$qname %in% want)
}
filter <- FilterRules(filter_factory(want))
dest <- filterBam(fname, tempfile(), filter=filter,
                  param=ScanBamParam(what="qname"))
readGAlignments(dest)
于 2014-06-13T13:11:43.823 に答える
1

使ってしまいましたsubset(DataFrame(scanBam(bam_file)[[1]]),qname %in% qname_vector)。これはまったく同じ構造 (リストのリスト) を維持するわけではありませんが、すべての情報が維持され、簡単にアクセスできます。

于 2014-06-13T15:53:20.233 に答える