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")