使用picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE（文件）讀取列表包含將被包含在輸出SAM或BAM文件中或從中排除的讀取的文件。默認值：null。

（1）不會很快，但是您可以為grep提供QNAMES文件。

  samtools視圖文件。 grep -f'qnames.txt >子集.sam  

其中qnames.txt具有

  EXAMPLE：QNAME1EXAMPLE：QNAME2EXAMPLE：QNAME3EXAMPLE：QNAME4EXAMPLE：QNAME5 代碼> 

（2）這會稍微複雜一點，但是您可以舉一個示例，說明grep可能具有正確的QNAME嗎？

（3）要保留BAM標頭我會使用 samtools -H file.bam > header.txt 來獲取標頭，然後使用grep文件

下面是一小段python代碼，顯示了一種幼稚但可管理的方式（注意，我希望 grep 更快，儘管讓其輸出標頭會很煩人）：

  import pysamqnames = set（...）＃讀取名稱到此處bam = pysam.AlignmentFile（“ some input file.bam”）obam = pysam。 alignmentFile（“ some output file.bam”，“ w”，template = bam）用於bam.fetch（until_eof = True）中的b：如果qname中的b.query_name：obam.write（b）bam.close（）obam。 close（） 

現在可以了，但是如果您要查看的名稱很多，那麼結果將變得非常慢（ grep -f也會如此）。 .. ）。如果需要按名稱選擇較大的讀取列表，則更有效的策略是：

1. 對BAM文件進行名稱排序（或picard） ，這可以用多個線程來完成。
2. 執行與上述相同的迭代，但僅與步驟2中列表中的最高讀取名稱進行比較（匹配後從堆棧中刪除該元素）。
3. li> ol>

對 $ n $ span>對齊方式執行grep到 $ m $ span> qnames會給你 $ O（mn）$ span>操作。如果您決定先對BAM和qname進行排序，則可以通過彈出讀取和qnames將其減少為 $ O（\ min \ {m，n \}）$ span>他們的堆棧。您也不會以這種方式一次讀取整個BAM，從而限制了內存消耗。我很確定Picard也會執行類似的操作。

這是我為此目的編寫的Python示例：

  import pysamdef filter_qname（ bamfile，idfile，outfile）：'''篩選器按查詢名稱讀取SAM / BAM文件參數---------- bamfile：str排序的BAM文件路徑idfile：str包含要保留的qname的文本文件路徑outfile：str要寫入'''的輸出文件＃讀入按名稱排序的BAM文件bam = pysam.AlignmentFile（bamfile，“ rb”）如果outfile.endswith（'。sam'）：output = pysam.AlignmentFile（outputfile， “ w”，template = bam）elif outfile.endswith（'。bam'）：output = pysam.AlignmentFile（outputfile，“ wb”，template = bam）否則：提高ValueError（“ outfile`的未知輸出文件格式： {}“。format（outfile））＃讀取要刪除的ID（僅保留唯一ID的list（set（...）））ids = list（set（set（[[l.rstrip（）for l in open（idfile ，'r'）。readlines（）]））＃排序ID以進行匹配BAM以進行有效處理（破壞性功能）ids.sort（key = str.lower，reverse = True）n_ids = len（ids）＃用於確保BAM按查詢名稱排序的變量last_q = None讀取= bam.fetch（until_eof = True ）read = next（reads），但為True：＃如果read名稱大於當前堆棧頂部（如果read.query_name > ids [-1]）：＃如果這是最後一個ID，如果n_ids == 1：＃將剩餘的讀取寫入新文件並退出while循環output.write（read）
對於bam中的r：output.write（read）break＃否則彈出該ID，再嘗試其他操作：ids.pop（）n_ids-= 1＃跳過與堆棧頂部匹配的讀取，如果是read.query_name == ids [- 1]：read = next（reads）＃如果讀取名稱小於堆棧頂部，則將其寫入並繼續進行其他操作：output.write（read）read = next（reads）output.close（） 

如果您想在實際中查看它，我已經在自己的​​其他生物信息學工具包此處中實現了它

您可以使用 SAMsift：

  samsift -i file.bam -0'q = open（“ qnames .txt“）。read（）。splitlines（）'-f'q中的QNAME' 

-i file.bam 指定輸入的BAM文件。 -0'q = open（“ qnames.txt”）。read（）。splitlines（）'會加載您的QNAME列表。 -f'q中的QNAME'指定用於對齊的過濾器。

使用SAMsift，可以使用任何Python表達式進行過濾。缺點是它比其他工具慢。

到目前為止，發布的大多數解決方案要么很慢，要么在某些情況下會產生虛假結果。 （我沒有測試Picard，我認為它會按預期工作。但是，就我個人而言，當我必須設置JVM時，我傾向於退縮！）

讓我們從另一個角度解決問題：

好的哈希函數的查找複雜度約為O（1）。如果我們可以創建QNAME的哈希表，則搜索BAM文件的記錄數應為O（n）。 Python開發人員為cpython中的哈希函數感到非常自豪，因此讓我們使用它：

 ＃！/ usr / bin / env python3import syswith open（sys.argv [1]，'r' ）作為索引文件：ids = set（l.rstrip（'\ r \ n'）對於索引文件中的l）對於sys.stdin中的行：qname，_ = line.split（'\ t'，1）如果qnames在ids中：sys.stdout.write（line） 

，測試集由一個〜3GB BAM文件和一個帶有〜65K條目的 qnames.txt 組成，運行 samtools視圖input.bam | ./idfilter.py qnames.txt 在我的測試服務器上大約需要 6.5秒。相比之下，將SAM傳遞到 fgrep -f qnames.txt （使用GNU grep）大約需要8.5秒（但可能會產生虛假結果）。

（為什麼 grep 會產生假結果嗎？如果您的 qnames.txt 中包含QNAME foo ，但是您的BAM文件還包含 fooX 和 fooY ，它們都將由 fgrep 匹配。解決方案是將您的每種搜索模式都固定在行首和製表符之間，例如 ^ foo \ t ，但是您必須使用標準的 grep （這將必須為每個模式構造一個NFA），而不是使用實現Aho-Corasick的 fgrep ，並且你會慢幾個數量級。）

我已經使用 std :: collections :: HashMap 的哈希函數重寫了我在Rust中的小Python腳本，並使用rust_htslib板條箱直接從BAM進行讀取和寫入，這是更多的方法即使讀寫BAM，速度也比Python快一倍。但是，我的Rust並不真正適合公眾消費……