对参考基因组按照200k分区间统计测序深度及GC含量
以前是自己写脚本： 【直播】我的基因组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也可以做到：
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 用法大全(一文就够吧)
而这个教程仍然是属于单细胞数据处理的一部分，大约十分之一的工作量，后面才有可能完成单细胞拷贝数全景图如下：
后续，我们的线下单细胞课程会细致入微的讲解这篇文章的图表复现，敬请期待！