使用bamdst完成公司外显子测序数据分析报告的重要环节
目前主流的ngs科研服务，包括WES和RNA-seq价格都是透明的，反正建库四五百块钱，测序都是60元一个G左右，也就是说10G数据量的wes也就是一千块钱，同理RNA-seq也是如此，而且有趣的是标准的生物信息分析报告通常是附赠的，因为也的确是数据出来后必备的，只不过跟你的科研结论还差十万八千里。但是看懂这个报告又是必须滴，比如每个样本的下面这些指标
1）Sample：样本编号；
2）Raw_reads：原始reads数；
3）Raw_base(G)：原始碱基数；
4）Clean_reads：经过滤后的有效reads数；
5）Clean_base(G)：经过滤后的有效碱基数；
6）Effective rate（%）：过滤得到的clean data碱基数占raw data碱基数的比例；
7）Q20(%)：测序质量值大于20的碱基占总碱基的比例；
8）Q30(%)：测序质量值大于30的碱基占总碱基的比例；
9）Error rate(%)：所有碱基的平均错误率；
10）GC ratio(%)：GC碱基的总含量。
比较好理解， 只要是稍微在我们生信技能树上面学习过，都明白，如果你还没有，那么应该是跟着我们的十万学员一起学习：
3.一万人陪你学习GEO数据库挖掘知识（公益视频听课笔4.记分享）
8.B站的11套生物信息学公益视频配套讲义，练习题及思维导图第一弹
并不是所有的指标都好理解，或者说好计算。
多样化的QC指标
比如：
（1） SampleID：样本编号；
（2） Total_Reads：clean data中双端总reads数目（read1和read2的总数）；
（3） Duplicates：重复的read数目（比例：重复的read数/clean read数）；
（4） Mapped_Reads：比对到参考基因组上的read数（比例）；
（5） Properly_mapped：成对双端read比对的方向和位置都正确的read数量（比例）；
（6） PE_mapped：成对的双端read均比对到参考基因组的read数量（比例）；
（7） SE_mapped：成对的双端read中只有一端比对到参考基因组的read数量（比例）；
（8） With_Mate_Mapped_to_Diff_Chr：成对的双端read分别比对到参考基因组的不同的染色体的read占总read数的比例；
（9） With_Mate_Mapped_to_Diff_Chr (mapQ>=5)：成对的双端read分别比对到参考基因组的不同的染色体的read，且比对质量大于等于5的read总数（比例）；
（10） Mean_Target_Depht: 在目标区域，平均测序深度；
（11） Target_Coverage：目标区域平均覆盖度，即被覆盖基因组区域区域占目标区域总大小的比例；
（12） Target_Coverage _above_4 X：目标区域中测序深度大于4X的碱基所占比例；
（13） Target_Coverage _above_10X：目标区域中测序深度大于10X的碱基所占比例；
（14） Target_Coverage _above_20X：目标区域中测序深度大于20X的碱基所占比例；
（15） Target_Coverage _above_50X：目标区域中测序深度大于50X的碱基所占比例；
（16） Target_Coverage _above_100X：目标区域中测序深度大于100X的碱基所占比例；
（17）Reads_on_target：捕获到外显子区域癿reads占Mapped reads数目及比例，即捕获率；
（18）Reads_Near_Target(500bp_flanking_region)：捕获到外显子区域的侧翼500bp范围的reads占Mapped reads数目及比例；
（19）Reads_On_and_Near_Target：捕获到外显子区域以及外显子侧翼500bp范围的总reads数占Mapped reads数目及比例。
这里面就涉及到了探针靶向区域的一些统计，不仅仅是区域本身，还可以统计其侧翼区域，或者是自己写脚本来做，我在4年前的WES教程就提到过, 教程目录如下：
-
WES（一）测序质量控制
-
WES（二）snp-calling
-
WES（三）snp-filter
-
WES（四）不同个体的比较
-
WES（五）不同软件比较
-
WES（六）用annovar注释
-
WES（七）看de novo变异情况
但是自己写的脚本不好用，而且局限性很大，因为当时毕竟不是专门搞这个开发。
最近在某美女的推荐下试用了 https://github.com/shiquan/bamdst 工具，效果非常棒，所以强烈安利一下。
安装bamdst
听过我的软件安装课程的朋友应该是对软件安装没啥疑问了，很简单的make即可，很标准的安装套路，如下：
cd ~/biosoft
git clone https://github.com/shiquan/bamdst.git
cd bamdst
make
~/biosoft/bamdst/bamdst --help
使用bamdst对bam文件进行批量统计
首先需要一个bed格式的外显子芯片捕获区域记录文件, 如果公司没有提供，你也可以使用CCDS记录的：
cat CCDS.20160908.txt |perl -alne '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort -i >exon_probe.hg38.gene.bed
批量运行代码是：
ls ../alignment/*.bam|while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $sample
mkdir $sample
~/biosoft/bamdst/bamdst -p exon_probe.GRCh38.gene.bed -o $sample $id
done
对每个样本统计出来的信息非常丰富！
样本汇总
当然，大部分情况我们做WES和RNA-seq 都是大样本量，所以每个样本的统计信息进行汇总是非常有必要的，而且还可以进行一些统计可视化。
这个时候大名鼎鼎的 multiqc 就可以派上用场了， 虽然它目前还没有 bamdst 的插件，但是我相信不远的将来，会有的。
或者可以自行写脚本进行汇总，当然，我也会发给作者看看，他是否还愿意继续开发维护，写multiqc的插件。
网友评论