使用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
的插件。