有很多種格式來記錄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文件。
