【直播】我的基因组47:测序深度和GC含量的关系 以前自己写脚本，太复杂，大家看不懂

下载hg38参考基因组

直接谷歌搜索即可：

拿到下载链接，近1G大小的文件，里面存放分类约3G的参考基因组。

http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

拿到hg38各个染色体长度

使用python里面的pyfaidx模块的faidx命令，代码如下

pip install pyfaidx
faidx hg38.fasta -i chromsizes > sizes.genome

结果如下：

chr1 248956422
chr10 133797422
chr11 135086622
chr12 133275309
chr13 114364328
chr14 107043718
chr15 101991189
chr16 90338345
chr17 83257441
chr18 80373285
chr19 58617616
chr2 242193529
chr20 64444167
chr21 46709983
chr22 50818468
chr3 198295559
chr4 190214555
chr5 181538259
chr6 170805979
chr7 159345973
chr8 145138636
chr9 138394717
chrM 16569
chrX 156040895
chrY 57227415

生成200Kb的区间bed文件

这里需要写脚本，我使用自己擅长的perl语言，代码如下：

cat sizes.genome |perl -alne '{print "$F[0]\t",200000*$_,"\t",200000*($_+1) foreach 0..$F[1]/200000-1}' > 200k.bed

输出文件开头如下：

chr1 0 200000
chr1 200000 400000
chr1 400000 600000
chr1 600000 800000
chr1 800000 1000000
chr1 1000000 1200000

值得注意的是线粒体长度是小于200K的，所以后面会出现警告：

Feature (chrM:0-200000) beyond the length of chrM size (16569 bp). Skipping.

不过问题不大。

使用bedtools统计200Kb的区间的基因组GC含量

因为使用的是bedtools这样成熟的轮子， 所以就是一行代码而已：

bedtools nuc -fi hg38.fa -bed 200k.bed | cut -f 1-3,5 > 200k_gc.bed
# 4_pct_at 5_pct_gc 6_num_A 7_num_C 8_num_G 9_num_T 10_num_N

文件如下：

$head 200k_gc.bed
#1_usercol 2_usercol 3_usercol 5_pct_gc
chr1 0 200000 0.420110
chr1 200000 400000 0.220065
chr1 400000 600000 0.315425
chr1 600000 800000 0.427140
chr1 800000 1000000 0.534730
chr1 1000000 1200000 0.608690
chr1 1200000 1400000 0.616775
chr1 1400000 1600000 0.567760
chr1 1600000 1800000 0.550655

使用bedtools统计200Kb的区间的测序情况

用的bedtools的multicov 这个小命令 ，代码如下：

bedtools multicov -bams alignment/SRR3081110.bam -bed 200k.bed > 200K_counts.bed

输出文件如下：

$head 200K_counts.bed
chr1 0 200000 53
chr1 200000 400000 44
chr1 400000 600000 58
chr1 600000 800000 146
chr1 800000 1000000 39
chr1 1000000 1200000 16
chr1 1200000 1400000 16
chr1 1400000 1600000 33
chr1 1600000 1800000 43
chr1 1800000 2000000 63