使用bamdst完成公司外显子测序数据分析报告的重要环节

Posted on by

使用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碱基的总含量。

比较好理解， 只要是稍微在我们生信技能树上面学习过，都明白，如果你还没有，那么应该是跟着我们的十万学员一起学习：

1.3个学生的linux视频学习笔记

2.生信人应该这样学R语言系列视频学习心得笔记分享

3.一万人陪你学习GEO数据库挖掘知识（公益视频听课笔4.记分享）

5.公共数据库挖掘视频学习心得体会

6.生信小技巧系列第一季完结版视频教程学习笔记分享

7.人类全外显子测序数据分析视频教程学习笔记

8.B站的11套生物信息学公益视频配套讲义，练习题及思维导图第一弹

9.转录组测序数据分析公益视频学习笔记分享

并不是所有的指标都好理解，或者说好计算。

多样化的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的插件。

​

Comments are closed.