测序完成得到的reads我们会比对到参考基因组得到bam文件，bam文件一般很大，很多时候我们只需要提取部分内容。
根据参考基因组位置提取
根据指定基因组区域的提取bam，可以使用以下命令。
samtools
samtools view -hb chr:start-end wgs.sort.bam > target.region.bam
# 根据bed文件来提取
samtools view -hb -L target.bed wgs.sort.bam > target.region.bam
bedtools
bedtools intersect -a wgs.sort.bam -b target.bed > target.region.bam
bedtools 是神器，参考 https://bedtools.readthedocs.io/en/latest/
sambamba
sambamba view -hb chr:start-end wgs.sort.bam > target.region.bam
# 根据bed文件来提取可以用 `sambamba slice `
sambamba slice -L target.bed wgs.sort.bam > target.region.bam
# sambamba view -L 报错，没有找到原因
sambamba 的很多用法跟samtools类似，可以参考 https://github.com/biod/sambamba
测评
选择5个基因的bed文件，对3.8G文件大小的WES的bam提取
target.bed
5 1253262 1295184 TERT
13 32889611 32973805 BRCA2
10 89622870 89731687 PTEN
17 41196312 41277500 BRCA1
17 7565097 7590856 TP53
samtools
bam=deduped.bam
if [ -f samtools_view_target.sam ];then
rm samtools_view_target.sam
fi
while read chr start end gene;do
samtools view $bam $chr:$start-$end >> samtools_view_target.sam
done < target.bed
资源消耗：0.50 user 0.10 system 0:04.74 elapsed 12%CPU
samtools view -L
samtools view -hb -L target.bed deduped.bam > samtools_view_L_target.bam
资源消耗：67.23 user 1.83 system 1:22.61 elapsed 83%CPU
sambamba
bam=deduped.bam
if [ -f sambamba_view_target.sam ];then
rm sambamba_view_target.sam
fi
while read chr start end gene;do
sambamba view $bam $chr:$start-$end >> sambamba_view_target.sam
done < target.bed
0.36 user 0.18 system 0:05.11 elapsed 10%CPU
sambamba slice -L
sambamba slice -L target.bed deduped.bam > sambamba_view_L_target.bam
资源消耗：0.08 user 0.03 system 0:06.85 elapsed 1%CPU
bedtools
bedtools intersect -a deduped.bam -b target.bed > bedtools_intersect_target.bam
282.51 user 2.11 system 4:47.22 elapsed 99%CPU
几个软件得到了相同结果，但是资源消耗各异。具体的资源统计如下表所示，利用sambamba的CPU使用是最低的，利用samtools 逐个提取速度是最快的，其次是sambamba，bedtools最慢。
但是我们必须要注意此次测评仅用了5个区域，如果用更多的区域，情况会大不一样，我猜测sambamba slice -L 会是速度最快资源消耗最少的。
883
