根据bed文件计算CG含量

看到一个临床小panel的测序数据的质控结果，CG含量 37%，（一般的结果都是45-55%的范围）那是实验和测序造成的，还是本来就是这个数值呢？
思路：根据bed文件提取基因组的的序列，然后再计算序列的CG含量。写脚本有点麻烦，下面这个方法稍微简单一点：

提取序列：

for region in $( awk '{print $1":"$2"-"$3}'  my_two_gene.bed ) ;do samtools  faidx  ucsc.hg19.fasta  ${region }; done > my_two_gene.fasta

samtools的1.9 版本http://www.htslib.org/doc/samtools.html支持批量提取 ，也可以写成：

awk  '{print $1":"$2"-"$3}'  my_two_gene.bed  > my_two_gene_region.txt
samtools  faidx  ucsc.hg19.fasta  -r my_two_gene_region.txt > my_two_gene.fasta

计算CG含量：

all=$(grep -v "^>"  my_two_gene.fasta | grep -o -i '[ATCG]'  | wc -l ) 
CG=$( grep -v "^>"  my_two_gene.fasta | grep -o -i '[CG]'  | wc -l )
echo "scale=3;${CG}/${all}" | bc

得到的结果，为0.359 , 和下机数据基本相符，看来GC含量本来就是这么低。