如果您不太希望精確讀取5000個數字，則可以使用一個samtools命令執行此操作：

  samtools視圖-bosubset.bam -s 123.4對齊方式。 bam chr1 chr2  

這將選擇40％的讀物（ .4 部分）（ 123 是種子，即便於重現）。方便的部分是，如果您具有配對末端讀段，它將使伴侶保持配對。對於每條染色體5000個讀數，只需將 .4 部分更改為足夠小的數字即可。

通常，您實際上不需要對標頭進行子集化。如果這樣做，某些工具的性能會更好一些，但是無論如何，您通常都會得到相同的結果。

您可以使用 SAMsift：

  samsift \ -i文件。bam\ -0'c = {“ chr1”：5000，“ chr2”：5000 }'\ -f'c [RNAME] >0'\ -c'c [RNAME]-= 1'\ -m不間斷刪除 

說明：

• -i file.bam –輸入文件
• -0'c = {“ chr1”：5000，“ chr2”：5000}' –初始化（為感興趣的染色體創建一個遞減計數字典）
• -f'c [RNAME] >0' –過濾條件（是當前染色體的計數器還是> 0？）
• -c'c [RNAME]-= 1' –代碼遞減當前染色體的計數（5000-> 4999-> ...- > 1-> 0-> -1-> ...）
• -m nonstop-remove –刪除導致Python錯誤且不停止的行（在這種情況下，訪問另一個染色體的不存在計數器（例如chr3）可能會導致錯誤。

有關更多信息，請參見 SAMsift自述文件。

我通常會推薦Devon Ryan的答案。但是，如果您確實希望每個Chr具有相同的數字，則可以使用以下python / pysam代碼（將從每個Chr輸出大約5000）：

  from pysam import AlignmentFilefrom隨機導入randomnreads = 5000infile = AlignmentFile（” mybam.bam“）outfile = AlignmentFile（” outbam.bam“，” wb“，template = infile）for inhr中的chr。reference_names：reads_in_chr = infile.count（chr）frac = float（nreads）/ reads_in_chr count = 0＃-將該循環替換為恰好5000次的第一次讀取-＃讀取infile.fetch（chr，multiple_iterators = True）：如果random（） < frac：count + = 1 outfile.write（read）print（“輸出的％i從％s讀取”％（count，chr））outfile.close（） 

染色體。如果只想使用，請說chr1和chr2用 [“ chr1”，“ chr2”] 替換 infile.reference_names 。

如果需要從每個chr精確地讀取5000次，但是不管它們是否是第一個，那麼您可以將內部for循環替換為：

 對於infile.fetch（chr，multiple_iterators = True）讀入：如果count > = nreads，則計數+ = 1：如果需要這些文件的伙伴，則打破outfile.write（read） 

同時讀取，您可以在 outfile.write（read）之後添加以下內容：

  mate_read = infile.mate（ read）outfile.write（mate_read） 

請注意，這很慢。

使用 samjdk ： http://lindenb.github.io/jvarkit/SamJdk.html

  $ java -jar dist /samjdk.jar --body -e \'Map<String，Integer> c = new HashMap<>（）;公共對象apply（SAMRecord r）{int n = c.getOrDefault（r.getContig（），0）; if（n> = 5000）返回false; c.put（r.getContig（），n + 1）; return r;}'\ input.bam  

如果您只想創建帶有標頭的截斷的BAM文件，則可以顯著簡化原始代碼：

 （samtools視圖-H input.bam; samtools視圖input.bam |頭-5000）| samtools -bo output.bam  

此命令避免了中間文件的調用，並減少了一些無用的中間命令的調用，而以子外殼為代價（由（…）調用）

如上所述，但格式如下：

 （samtools視圖-H input.bam samtools視圖input.bam | head -5000＃（*））\ | samtools -bo output.bam  

…當然，可以通過將上面標有（*）的行替換為一個或多個，將其擴展為通過多條染色體進行過濾按染色體名稱子集劃分的行（ x ，如果已為原始BAM文件建立索引，則不需要 grep ！）