需求将fastq按bed文件进行拆分。
首先要按bed文件的区域进行拆分，必然是要引入位置信息的。所以，先将fastq文件与参考基因组进行比对。比对后得到bam文件，这时，使用bedtools intersect就可以很方便的将bam文件根据bed文件拆开。
bedtools intersect -abam sample.bam -b target.bed > sample.target.bam
如果加入-v参数，则是不要target区域。
bedtools intersect -abam sample.bam -b target.bed -v > sample.nontarget.bam
在得到需求区域的bam文件后，再还原为fastq文件。
bedtools和samtools都有现成的将bam还原为fastq的方法：
首先是bedtools的bamtofastq
bedtools bamtofastq -i sample.target.bam \
-fq sample.target_1.fq.gz \
-fq2 sample.target_2.fq.gz
如果用samtools
samtools fastq -1 sample.target_1.fq -2 sample.target_2.fq \
-0 /dev/null -s /dev/null \
-n -F sample.target.bam
上面两种方法，都存在一定的问题，就是read name会改变，因为有时fastq第一行的”@”后面跟的信息，如果存在空格，在比对后bam文件中只会剩下空格最前的部分作为read name。而bedtools和samtools根据这个bam文件来处理，就会丢失了这部分信息。同时，bedtools在还原时不停的报找不到pair的Warning，大概是因为双端的reads不对到了不同位置，而只有一条reads在target区域（存疑）。因此还原出来的fastq文件特别小。samtools没有报Warning，但是还原出来的文件也是一样小，我觉得可能与bedtools原理相同。
事实上，只要一端的reads能比对到target区域，这个双端reads都要保留。那么，其实可以先把read name抠下来，再从原始的fastq里根据read name把fastq过滤出来。
抠read name非常简单
samtools view sample.target.bam | awk '{print $1}' | sort | uniq > name.list
然后使用seqtk。
seqtk subseq sample_1.fq.gz name.list > sample.target_1.fq
seqtk subseq sample_2.fq.gz name.list > sample.target_2.fq
gzip sample.target_1.fq
gzip sample.target_2.fq
即可。当然这其实就不是bam文件转为fastq了，但是结果表现更好，read name也是完整的。