【原】对参考基因组按照200k分区间统计测序深度及GC含量

以前是自己写脚本： 【直播】我的基因组47:测序深度和GC含量的关系 可能是太复杂，大多数读者表示看不懂，所以我重新使用已有的轮子来做这件事。

下载hg38参考基因组

直接谷歌搜索即可：

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

http://hgdownload.cse./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也可以做到：

bedtools   makewindows -g  sizes.genome  -w  200000 > 200k.bed


而且还完美的解决了我自己的perl脚本的无关痛痒的小bug

使用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

要想顺利完成这个小任务，其实需要掌握好bedtools：bedtools 用法大全(一文就够吧)

而这个教程仍然是属于单细胞数据处理的一部分，大约十分之一的工作量，后面才有可能完成单细胞拷贝数全景图如下：

后续，我们的线下单细胞课程会细致入微的讲解这篇文章的图表复现，敬请期待！