statistic GC content by interval. BED format file include arget interval information. BEDTool statistic read number and extract sequence, awk statistic GC content.

bedtools map -a interval.bed -b sample.bam -c 10,10 -o count,concat | awk -v OFS="t" '{n=length($5); gc=gsub("[gcGC]", "", $5); print $1,$2,$3,$4,gc/n}'

思路:利用bedtools的map工具,首先找到map到interval.bed中的每个interval的reads的序列,然后统计这些序列中有多少GC。

缺点是,该map找到的是覆盖到interval的reads的序列,有可能read有一半与interval重叠,但bedtools会把该reads的所有序列都提出来。

PS:Bedtools是基因组特征genome feature比对工具,用起来非常非常的爽。

#####################################################################

#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任

#Author: Jason Zhu

#####################################################################