这个完全是项目实战经验分享咯,有大样本量NGS多组学数据处理经验的朋友应该能很容易理解,动辄几个T的数据,上百个样本很难一个个的检查是否出现问题,需要一个简单方便快捷质控方案。而我认为qualimap+multiqc完美解决多组学比对结果的质控,当然也欢迎大家在我们生信技能树平台推荐自己的实战经验!
qualimap的英文文档本来就很清楚,但是需要一些时间来理解: http://qualimap.bioinfo.cipf.es/doc_html/analysis.html
qualimap的测试数据也给的很全面: http://qualimap.bioinfo.cipf.es/doc_html/samples.html#counts-example-output
需要自己制作 mm10.exon.chr.bed
文件,我在生信菜鸟团博客多次讲解过它的制作方式,如下:
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS//archive/21/CCDS.20161208.txt
cat CCDS.20161208.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 >mm10.exon
awk '{print $0"\t0\t+"}' mm10.exon >mm10.exon.bed
然后就可以运行 qualimap,下面的shell脚本我也在生信菜鸟团博客多次讲解过,相信你们肯定能看懂了。
## $1 for the config file: bam_path.txt
## $2 and $3 for submit jobs.
exon_bed='/home/jianmingzeng/annotation/CCDS/mouse/mm10.exon.chr.bed'
qualimap='/home/jianmingzeng/biosoft/Qualimap/qualimap_v2.2.1/qualimap'
cat $1 |while read id
do
echo $id
if((i%$2==$3))
then
$qualimap bamqc --java-mem-size=20G -gff $exon_bed -bam $id
fi
i=$((i+1))
done
可以看外显子的测序情况。
这里其实应该是首推RSeQC这个软件,可惜那是个python的,而且运行超慢,还具耗费内存。所以不得已转为
示例报告: http://kokonech.github.io/qualimap/kidney_rnaseqqc/qualimapReport.html
## $1 for the config file: bam_path.txt
## $2 and $3 for submit jobs.
gtf='/home/jianmingzeng/reference/gtf/gencode/gencode.v25.annotation.gtf'
qualimap='/home/jianmingzeng/biosoft/Qualimap/qualimap_v2.2.1/qualimap'
cat $1 |while read id
do
file=$(basename $id )
sample=${file%%.*}
echo $sample
if((i%$2==$3))
then
$qualimap rnaseq --java-mem-size=20G -gtf $gtf -bam $id -pe -oc $sample
fi
i=$((i+1))
done
属于转录组数据质控的一部分,比如:6 samples in 2 conditions 的报告,这个时候的input数据是表达矩阵了:
示例报告:
qualimap会给每一个样本单独进行质控,得到质控的html报告,不过那个报告本身非常丑陋,而且单独的html报告依然是不方便浏览,需要归纳汇总,这个时候multiqc就能大展身手。multiqc已经发展成为了一个质控平台,大家可以在其平台上面开发各种质控软件的可视化归纳汇总插件,而qualimap就已经被开发了。
对WES数据汇总,下面我截图其中一个例子:
对RNA-seq数据汇总,我截图两个例子:
因为我这里展示的公共数据的质控结果,所以非常优秀,但实际在处理自己的真实数据,总是会发现各种各样的问题。