samtools idxstats

有很多种格式来记录mapping的结果,大多数都收录在UCSC的帮助文档中。在mapping之后,如何对mapping的结果进行评估呢? 最简单的,就是通过samtools来评估mapping质量了。

samtools idxstats aln.sorted.bam

注意,这一步之前需要经过sort和index。结果会显示:

 chr1 195471971 6112404 0
 chr10 130694993 3933316 0
 chr11 122082543 6550325 0
 chr12 120129022 3876527 0
 chr13 120421639 5511799 0
 chr14 124902244 3949332 0
 chr15 104043685 3872649 0
 chr16 98207768 6038669 0
 chr17 94987271 13544866 0
 chr18 90702639 4739331 0
 chr19 61431566 2706779 0
 chr2 182113224 8517357 0
 chr3 160039680 5647950 0
 chr4 156508116 4880584 0
 chr5 151834684 6134814 0
 chr6 149736546 7955095 0
 chr7 145441459 5463859 0
 chr8 129401213 5216734 0
 chr9 124595110 7122219 0
 chrM 16299 1091260 0
 chrX 171031299 3248378 0
 chrY 91744698 259078 0
 * 0 0 0

其中第一列是染色体名称,第二列是序列长度,第三列是mapped reads数,第四列是unmapped reads数。 如果是RNAseq,我们可以使用broad institute的RNA-SeQC来得到更加完整的报告。下载到文件之后,也许需要安装BWA来获取更精准的结果,但是如果不安装的话,也可以进行分 析。一般来说,这一步不需要特别精准的结果,所以我很少使用BWA选项。下载的文件如果是.zip结尾的,直接把它改写成.jar就可以运行了。 在它的主页上下载所需要的Example RNA-seq Data。下载结束之后,该解压的解压缩。接下来运行:

 samtools index example/ThousandReads.bam
 samtools faidx example/Homo_sapiens_assembly19.fasta
 java -Xmx2048m -jar RNA-SeQC_v1.1.7.jar -n 1000 -s "TestId|example/ThousandReads.bam|TestDesc" -t example/gencode.v7.annotation_goodContig.gtf -r example/Homo_sapiens_assembly19.fasta -o ./testReport/ -start gc -gc example/gencode.v7.gc.txt

以上的参数只有一个与其说明文档不一样的地方就是使用了-Xmx2048m来指定java虚拟机的内存大小为2G。如果遇到java.lang.OutOfMemoryError,还可以指定得再大些。

当然如果是自己的文件的话,还需要多两步:

1.BAM,reference及GTF文件的基因组名称必须一致。

2.需要使用picard工具包中的CreateSequenceDictionary来构建一个dictionary文件。

原文地址:https://www.cnblogs.com/pennyy/p/4650690.html