bam格式转换为Fastq/Fasta格式

Samtools Fastq
GATK SamToFastq
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

samtools fastq -@ $(nproc) -1 D1_1.fq.gz -2 D1_2.fq.gz -s /dev/null -0 /dev/null D1.rn.sort.bam -c 9

-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

zcat D1_1.fq.gz | awk 'NR%4==1 && $1!~/^@/ {print; print NR; exit}'

b. 压缩文件是否完整

gunzip -t D1_1.fq.gz

GATK SamToFastq

gatk SamToFastq --INPUT D1.rn.sort.bam --FASTQ D1_1.fq.gz --SECOND_END_FASTQ D1_2.fq.gz --TMP_DIR . --INCLUDE_NON_PRIMARY_ALIGNMENTS true --COMPRESSION_LEVEL 9

注意，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端的大小也可能不一样。