bam格式轉換為Fastq/Fasta格式

1. Samtools Fastq
2. GATK SamToFastq
3. Bedtools bamtofastq

舉例說明，比如說我們現在有一個轉錄組比對檔案D1_D1.sort.bam：

samtools view -H D1_D1.sort.bam | tail
@SQ	SN:19	LN:58617616
@SQ	SN:20	LN:64444167
@SQ	SN:21	LN:46709983
@SQ	SN:22	LN:50818468
@SQ	SN:X	LN:156040895
@SQ	SN:Y	LN:57227415
@SQ	SN:MT	LN:16569
@RG	ID:D1_D1_D1	PL:illumina	PU:D1	LB:D1	SM:D1	CN:BGI
@PG	ID:bwa	PN:bwa	VN:0.7.10-r789	CL:/home/lixf/RNA_module_V1.0/Alignment/Alignment_BWA/bin/bwa mem -t 4 -a -R @RG\tID:D1_D1_D1\tPL:illumina\tPU:D1\tLB:D1\tSM:D1\tCN:BGI /data_center_11/Project/wenpp/RNA/01_RY20190826D002_RNAref/prepare/Homo_sapiens.GRCh38.dna.toplevel.fa /data_center_11/Project/wenpp/RNA/01_RY20190826D002_RNAref/01.deal_reads/process/Filter_SOAPnuke/D1/D1/D1_1.fq.gz /data_center_11/Project/wenpp/RNA/01_RY20190826D002_RNAref/01.deal_reads/process/Filter_SOAPnuke/D1/D1/D1_2.fq.gz
@PG	ID:samtools	PN:samtools	PP:bwa	VN:1.12	CL:samtools view -H D1_D1.sort.bam

           

從@PG可以看出來比對軟體是BWA，這是一個與hg38 toplevel基因組序列比對後産生的，根據染色體上位置進行排序的BAM檔案。那麼如果我們要與不同的基因組序列進行比對，比如更換基因組序列的版本為hg19，就需要将BAM檔案還原成原來的fastq檔案，而從上述資料可以看出，原始測序應該是一個雙端測序檔案，那麼我們應該轉換後得到的是兩個fastq檔案。

Samtools Fastq

1. 重新根據read name進行排序

原本D1_D1.sort.bam檔案的序列為：

samtools view D1_D1.sort.bam | head | cut -f 1,2        
ST-E00600:147:HNTKYALXX:1:1144:13675:27524	337
ST-E00600:147:HNTKYALXX:1:2144:15700:30405	369
ST-E00600:147:HNTKYALXX:1:2124:6876:4977	321
ST-E00600:147:HNTKYALXX:1:2124:6876:4977	401
ST-E00600:147:HNTKYALXX:1:2361:24813:5415	337
ST-E00600:147:HNTKYALXX:1:2361:23267:4805	337
ST-E00600:147:HNTKYALXX:1:2240:12274:12759	97
ST-E00600:147:HNTKYALXX:1:2344:3260:11835	353
ST-E00600:147:HNTKYALXX:1:2375:15266:19476	417
ST-E00600:147:HNTKYALXX:1:1374:16694:4476	385
           

可以看到read name是沒有規律的，因為這時BAM檔案是按照reads在染色體上比對的位置排序的，這樣就不友善轉換，因為雙端測序的reads應該是成對排列的。

是以，要先按照read name進行排序，根據定義，我們可以知道read name各部分的含義：

@           : 代表read name開頭
ST-E00600   : 測序平台的編号
147         : Run的編号
HNTKYALXX   : flow cell的編号
1           : lane的編号
1144        : Tile的編号 
13675       : tile中所處cluster的x軸坐标
27524       : tile中所處cluster的x軸坐标
## 比對軟體也可能在read name最後添加:1/:2代表來自哪一端的測序。
           

比如說bam中的這兩條read：

ST-E00600:147:HNTKYALXX:1:2124:6876:4977	321
ST-E00600:147:HNTKYALXX:1:2124:6876:4977	401
           

第二列的數字代表bitwise flags，标記該read的比對情況。可以用

samtools flags

進行換算。

samtools flags 321
0x141	321	PAIRED,READ1,SECONDARY
samtools flags 401
0x191	401	PAIRED,REVERSE,READ2,SECONDARY
           

這是什麼意思呢？

代表第1條read是雙端測序，屬于R1端，而且不是primary alignment；

第2條read也是雙端測序，Reverse表示比對到的是參考序列的負鍊，屬于R2端，不是primary alignment。

回顧了以上的内容，我們就可以對bam檔案進行重新排序：

samtools sort -n -@ $(nproc) -o D1.rn.sort.bam D1_D1.sort.bam
           

$(nproc)代表伺服器可用的線程數，多線程處理會更快。

轉換後的bam檔案：

samtools view D1.rn.sort.bam| head | cut -f 1,2
ST-E00600:147:HNTKYALXX:1:1101:1036:18035	99
ST-E00600:147:HNTKYALXX:1:1101:1036:18035	147
ST-E00600:147:HNTKYALXX:1:1101:1036:18035	2195
ST-E00600:147:HNTKYALXX:1:1101:1036:18035	403
ST-E00600:147:HNTKYALXX:1:1101:1036:20948	83
ST-E00600:147:HNTKYALXX:1:1101:1036:20948	163
ST-E00600:147:HNTKYALXX:1:1101:1036:21042	83
ST-E00600:147:HNTKYALXX:1:1101:1036:21042	163
ST-E00600:147:HNTKYALXX:1:1101:1036:21637	83
ST-E00600:147:HNTKYALXX:1:1101:1036:21637	163
           

2. 轉換為Fastq

-1 代表輸出的R1端Fastq.gz檔案名

-2 代表輸出的R2端Fastq.gz檔案名

-s singleton代表将無法找到mate pair的reads輸出到檔案，這裡不需要，是以輸出到/dev/null

-0 将bitwise flags中同時有READ1/READ2标記或缺乏任何READ1/READ2标記的reads輸出到檔案，這裡不需要，是以輸出到/dev/null

-c 如果輸出的是fastq.gz，必須使用-c 選項，選擇壓縮級别，一般選擇9，最高壓縮

*比對軟體某些情況下會将一端比對上，另一端無法比對上的reads保留到bam檔案，這樣的reads稱為singleton。

3. 檢查Fastq檔案

最好檢查轉換來的Fastq檔案是否格式正确，并且沒有受到損壞

a. 是否有非@開頭的identifier lines

b. 壓縮檔案是否完整

gunzip -t D1_1.fq.gz
           

GATK SamToFastq

注意，GATK方法僅适用于無重複read name的bam檔案，也就是說read name的字尾中必須要标明:1/:2，不然GATK會報錯。也就是說D1.rn.sort.bam其實是不能用這個方法的。

bedtools bamtofastq

bedtools bamtofastq -i D1.rn.sort.bam -fq D1.1.fq -fq2 D1.2.fq 2>>bamtofastq.log
           

注意，bedtools方法在singleton的reads中會報錯

*****WARNING: Query ST-E00600:147:HNTKYALXX:1:1272:31566:33771 is marked as paired, but its mate does not occur next to it in your BAM file.  Skipping. 
           

也就是說它隻會保留paired的reads，是以生成的fastq檔案可能會比應有的大小要小，而且R1端與R2端的大小也可能不一樣。