karyoploteR: 画一个测序深度分布图
测序深度的分布图是很常见的一类图，可用于发现染色体片段或整条染色体的缺失或是重复。如果是缺失，测序时捕获的概率就小，最终比对呈现出来的就是有一段区域深度明显偏低；如果是重复，捕获概率就大，测的reads就多，比对得到的深度自然就高。临床上做CNV检测原理也和这差不多。
上一次介绍这个R包的时候，在文末提了一个小问题，非人物种可以用这个R包吗？（见karyoploteR: 将基因组信息在染色体上“画”出来）
答案是可以！比如今天我想画一个植物的测序深度分布图，最终的效果图如下。
下面简单写一下过程
1. 求出每一个碱基上的深度
这一步用到的是samtools的一个模块depth。
samtools depth -b ref.fa.bed xxx.smr.bam > bedbase_depth.txt
这里的bed是含有每一条染色体位置信息的bed文件，左闭右开，从0开始，得到的结果是每一个碱基上的深度，如下。
head -n 4 bedbase_depth.txt
NC_027757.2 32 1
NC_027757.2 33 1
NC_027757.2 34 2
NC_027757.2 35 2
注意：该文件的第二列是0开始还是1开始
怎么看？samtools tview xxx.smr.bam查看一下第一个碱基是在哪个位置即可。结果是第一个碱基就在32位处，故是以1开始的。建议转换为bed坐标，即减1。
以0开始还是以1开始这个地方容易搞错
2. 划定窗口，计算总碱基数和平均深度
接下来的思路是，先分染色体，再对第二列除以50kb结果取整，当作窗口的NO.编号（这时应该想一想每条染色体能分为多少个窗口，窗口起止位置和窗口编号的对应关系，下面的脚本得到的final.all_chr_to_50kb.bed是bed坐标的）。在bedbase_depth.txt文件中，当窗口号相同时，第三列相加，当作该窗口的总碱基数。
划定窗口
$ cat picture.all_chr_to_50kb.pl
#! /ifs1/Software/miniconda3/bin/perl
use warnings;
use strict;
my $parm_name = $ARGV[0];
open my $file_fh, "<", "$parm_name";
open my $file_fh2, ">", "all_chr_to_50kb.bed";
while (<$file_fh>) {
chomp $_;
my @one_line = (split("\t", $_));
my $num = int($one_line[2]/50000);
for(my $i=1; $i <= $num; $i++) {
my $left = 50000*($i-1);
my $right = 50000*($i);
print $file_fh2 "$one_line[0]\t$left\t$right\n";
}
my $left2 = $num*50000;
my $right2 = $one_line[2];
print $file_fh2 "$one_line[0]\t$left2\t$right2\n";
}
close $file_fh;
close $file_fh2;
$ perl picture.all_chr_to_50kb.pl ref.fa.bed
$ head -n 2 all_chr_to_50kb.bed
NC_027757.2 0 50000
NC_027757.2 50000 100000
$ less all_chr_to_50kb.bed | awk '{print $1"\t"$2"\t"$3"\t"$2/50000+1}' > final.all_chr_to_50kb.bed
$ head -n 2 final.all_chr_to_50kb.bed #第四列表示第几个窗口
NC_027757.2 0 50000 1
NC_027757.2 50000 100000 2
计算总碱基数
$ awk '{print $1"_"int(($2-1)/50000)"\t"$3}' bedbase_depth.txt > tmp1 #上面也提到了，建议将位置转成bed坐标，所以减一。
$ head tmp1
NC_027757.2_0 1
NC_027757.2_0 1
NC_027757.2_0 2
NC_027757.2_0 2
接下来用perl脚本结合哈希表可求出窗口和窗口上总的碱基数（除以50000就是平均测序深度了）。
$ cat 1.pl
#! /ifs1/Software/miniconda3/bin/perl
use warnings;
use strict;
my %chuangkou_jianji=();
open my $file_fh, "<", "tmp1";
while (<$file_fh>) {
chomp $_;
my @one_line = (split("\t",$_));
if (exists $chuangkou_jianji{$one_line[0]}) {
$chuangkou_jianji{$one_line[0]} = $chuangkou_jianji{$one_line[0]} + $one_line[1];
} else {
$chuangkou_jianji{$one_line[0]} = $one_line[1];
}
}
close $file_fh;
open my $file_fh2, ">", "tmp2";
foreach my $i ( sort keys %chuangkou_jianji) {
print $file_fh2 "$i\t$chuangkou_jianji{$i}\n";
}
close $file_fh2;
结果如下，每一行表示某条染色体的某个窗口有多少个碱基。
$ head -n 2 tmp2
NC_027757.2_0 592839
NC_027757.2_1 510848
计算平均深度
稍微调整一下，就能得到某条染色体第几个窗口的平均测序深度了。
$ awk -F "_| " '{print $1"_"$2"\t"$3"\t"$4}' tmp2 | awk '{print $1"\t"$2+1"\t"$3/50000}' | sort -k1,1 -k2,2n > tmp3
$ head -n 5 tmp3
NC_027757.2 1 11.8568
NC_027757.2 2 10.217
NC_027757.2 3 12.202
NC_027757.2 4 21.4317
NC_027757.2 5 26.7172
接下来为了画图，需要知道窗口的物理位置。可以用上面tmp3的窗口编号直接计算得到，但这样的话，每一条染色体的最后一个窗口的坐标是错的，尽管对整体画图没有什么影响。为了准确起见，采用哈希表匹配。
#根据上文的final.all_chr_to_50kb.bed得到的哈希表
NC_027757.2-1 0-50000
NC_027757.2-2 50000-100000
#然后再用脚本匹配，在tmp3的基础上得到如下结果
NC_027757.2 0 50000 11.8568
NC_027757.2 50000 100000 10.217
再替换染色体名称就可以了。
由于高度同源序列和重复序列的存在，某些区域的深度会非常高，建议人为调整以便于呈现。如，将平均深度大于100的区域深度全重新定义为100。
awk '{if($4 > 100) print $1"\t"$2"\t"$3"\t100"}{if($4 <= 100) print $0}' depth_distribution.bed > final_depth_distribution.bed
接下来才是本文的主角。
3. 画图
#安装，加载R包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("karyoploteR", version = "3.8")
library(karyoploteR)
Bn.genome <- toGRanges("ref.fa.rename.bed") #染色体名称都改成了chr1A、chr2A......
kp <- plotKaryotype(genome=Bn.genome,labels.plotter = NULL,plot.type=1,main="sample") #labels.plotter = NULL不显示染色体名称，下一行命令是添加染色体名称，这样做的好处是名称的文字大小可以自己设置
kpAddChromosomeNames(kp, cex = 1.2)
sample_bar <- toGRanges("final_depth_distribution.bed") #前面脚本得到的结果文件
kpBars(kp, sample_bar,border="#377EB8",data.panel = 1,r0 = 0,r1 = 1,ymin = 0,ymax = 100)
#下面两行是在染色体的两边加上y轴，左边还是右边由参数side区分，1表示右边，2表示左边，默认是1
#tick.len = 700000, label.margin = -600000这两个值具体多少合适与染色体的长度有关，多试几次就可以了，它们的单位是bp，序列长度的那个bp
kpAxis(kp, data.panel = 1, ymin = 0, ymax=100, tick.pos = c(20,40,60,80), labels = c("20","40","60","80"),tick.len = 700000, label.margin = -600000, side = 2, cex = 0.4, col="#777777")
kpAxis(kp, data.panel = 1, ymin = 0, ymax=100, tick.pos = c(20,40,60,80), labels = c("20","40","60","80"),tick.len = 700000, label.margin = -600000, cex = 0.4,col="#777777")
