有很多種格式來記錄mapping的結(jié)果，大多數(shù)都收錄在UCSC的幫助文檔中。在mapping之后，如何對(duì)mapping的結(jié)果進(jìn)行評(píng)估呢？ 最簡(jiǎn)單的，就是通過samtools來評(píng)估m(xù)apping質(zhì)量了。

samtools idxstats aln.sorted.bam

注意，這一步之前需要經(jīng)過sort和index。結(jié)果會(huì)顯示：

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

其中第一列是染色體名稱，第二列是序列長(zhǎng)度，第三列是mapped reads數(shù)，第四列是unmapped reads數(shù)。 如果是RNAseq，我們可以使用broad institute的RNA-SeQC來得到更加完整的報(bào)告。下載到文件之后，也許需要安裝BWA來獲取更精準(zhǔn)的結(jié)果，但是如果不安裝的話，也可以進(jìn)行分 析。一般來說，這一步不需要特別精準(zhǔn)的結(jié)果，所以我很少使用BWA選項(xiàng)。下載的文件如果是.zip結(jié)尾的，直接把它改寫成.jar就可以運(yùn)行了。 在它的主頁上下載所需要的Example RNA-seq Data。下載結(jié)束之后，該解壓的解壓縮。接下來運(yùn)行：

samtools index example/ThousandReads.bamsamtools faidx example/Homo_sapiens_assembly19.fastajava -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

以上的參數(shù)只有一個(gè)與其說明文檔不一樣的地方就是使用了-Xmx2048m來指定java虛擬機(jī)的內(nèi)存大小為2G。如果遇到j(luò)ava.lang.OutOfMemoryError，還可以指定得再大些。

當(dāng)然如果是自己的文件的話，還需要多兩步：

1.BAM，reference及GTF文件的基因組名稱必須一致。

2.需要使用picard工具包中的CreateSequenceDictionary來構(gòu)建一個(gè)dictionary文件。

轉(zhuǎn)載于:https://www.cnblogs.com/pennyy/p/4650690.html

總結(jié)

以上是生活随笔為你收集整理的samtools idxstats的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。