當測序得到的fastq文件map到基因組之后,我們通常會得到一個sam或者bam為擴展名的文件。SAM的全稱是sequence alignment/map format。而BAM就是SAM的二進制文件(B取自binary)。 那么SAM文件的格式是什么樣子的呢?如果你想真實地了解SAM文件,可以查看它的說明文檔。SAM由頭文件和map結果組成。頭文件由一行行以@起始的注釋構成。而map結果是類似下面的東西:
HWI-ST1001:137:C12FPACXX:7:1115:14131:66670 0 chr1 12805 1 42M4I5M * 0 0 TTGGATGCCCCTCCACACCCTCTTGATCTTCCCTGTGATGTCACCAATATG CCCFFFFFHHGHHJJJJJHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJ AS:i:-28 XN:i:0 XM:i:2 XO:i:1XG:i:4 NM:i:6 MD:Z:2C41C2 YT:Z:UU NH:i:3 CC:Z:chr15 CP:i:102518319 XS:A:+ HI:i:0 HWI-ST1001:137:C12FPACXX:7:2313:17391:30032 272 chr1 13494 1 51M * 0 0 ACTGCCTGGCGCTGTGCCCTTCCTTTGCTCTGCCCGCTGGAGACAGTGTTT CFFFFHHJJJJIJJJJIJJJJJJJJJJJJJJJJJJJJJHHHHFA+FFFC@B AS:i:-3 XN:i:0 XM:i:1 XO:i:0 XG:i:0NM:i:1 MD:Z:44G6 YT:Z:UU XS:A:+ NH:i:3 CC:Z:chr15 CP:i:102517626 HI:i:0 HWI-ST1001:137:C12FPACXX:7:1109:17518:53305 16 chr1 13528 1 51M * 0 0 CGCTGGAGCCGGTGTTTGTCATGGGCCTGGGCTGCAGGGATCCTGCTACAA #############AB=?:*B?;A?<2+233++;A+A2+<7==@7,A<A<=> AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0NM:i:2 MD:Z:8A21T20 YT:Z:UU XS:A:+ NH:i:4 CC:Z:chr15 CP:i:102517592 HI:i:0
看上去很類似fastq文件,它也有read名稱,序列,質量等信息,但是又不完全一樣。首先,每個read只占一行,只是它被tab分成了很多列,一共有12列,分別記錄了:
1. read名稱
2. SAM標記
3. chromosome
4. 5′端起始位置
5. MAPQ(mapping quality,描述比對的質量,數字越大,特異性越高)
6. CIGAR字串,記錄插入,刪除,錯配以及splice junctions(后剪切拼接的接頭)
7. mate名稱,記錄mate pair信息
8. mate的位置
9. 模板的長度
10. read序列
11. read質量
12. 程序用標記
顯然,其中chromosome至CIGAR的信息都是非常重要的。但是這些對我們不重要,我們只需要了解SAM/BAM文件是什么,就可以了。重要的是如果進行下游的操作。 要操作SAM/BAM文件,首先需要安裝samtools。它的安裝過程和所有的linux/unix程序一樣,都是經過make之后生成可執行程序,然后把它的路徑告知系統,或者放在系統可以找到的位置就可以了。 比如:
tar zxvf samtools-0.1.18.tar.bz2
cd samtools-0.1.18/
make
samtoolpath=`pwd`
PATH=PATH:$samtoolpath
然后就可以按照samtools主頁上介紹的工具進行各種操作了。我們最常見的幾步操作比如 0. SAM,BAM轉換
samtools view -h file.bam > file.sam
samtools view -b -S file.sam > file.bam
1. sorting BAM文件。大多數下游程序都要求BAM文件是被排過序的。
samtools sort file.bam outputPrefix
2. 創建BAM index。這也是被大多數下游程序所要求。
samtools index sorted.bam
3. index模板基因組。這也是被大多數下游程序所要求。
samtools faidx Homo_sapiens_assembly19.fasta
在很多時候,我們還會看到一種擴展名為BED的mapping文件。其具體格式也是幾經變化,但是現在以UCSC的描述為准。從BAM文件轉換成BED文件,我們需要安裝BEDtools。下載安裝就不多說了。示例一個如何從BAM文件轉換成BED文件的命令:
bamToBed -i reads.bam > reads.bed
更多的具體內容可以參見其說明文檔。 當然,還有很多種格式來記錄mapping的結果,大多數都收錄在UCSC的幫助文檔中。比如上次有人問及的.bw是什么文件(bigWig文件)之類的,都可以在那里找到答案。 上次談及fastq文件時,有講過其質量評估的問題,那么在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文件。
原文來自:http://pgfe.umassmed.edu/ou/archives/3050