如何从BAM文件中提取fastq

虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件，但偶尔我们也需要提取BAM文件中特定区域中fastq。最开始我认为这是一个非常简单的操作，因为samtools其实已经提供了相应的工具samtools fastq.

以biostar handbook的Ebola病毒数据为例，首先获取比对得到的BAM文件。

# 建立文件夹
mkdir -p refs
# 根据Accession下载
ACC=AF086833
REF=refs/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF
bwa index $REF
SRR=SRR1553500
fastq-dump -X 100000 --split-files $SRR
BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

让我们尝试提取下前1kb比对的fastq文件

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools fastq -1 read_1.fq -2 read_2.fq -
#[M::bam2fq_mainloop] discarded 0 singletons
#[M::bam2fq_mainloop] processed 8702 reads

看起来似乎我们成功了，那让我们把这些序列比对回去吧

bwa mem refs/AF086833.fa read_1.fq read_2.fq
[mem_sam_pe] paired reads have different names: "SRR1553500.11772", "SRR1553500.43775"

你会发现这里产生了错误。这个错误说的是配对的read名字不同。导致问题的原因是：samtools fastq会按照fastq在bam里原本的顺序进行提取，而按照位置排序的bam的两个配对read并不是前后紧挨着的。所以我们需要进行一波排序。

cat read_1.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_1.fq.gz
cat read_2.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_2.fq.gz

这里将5个shell命令用管道连起来完成这个任务，咋看起来比较复杂，但是理解每一步的目的，每个参数的作用，以后就能活用起来了。

首先用cat逐行读取read，传入到paste中。paste的作用是将多个文件合并成一个文件，做法比较简单粗暴，举个例子

cat A.txt
A
B
C
cat B.txt
D
E
F
paste A.txt B.txt
A   D
B   E
C   F

cat read_1.fq | paste - - - -表示读取4行合并成一行。相当于fastq里的read都会用一行而不是四行表示。接着用sort按照第一行(-k1,1)排序，内存的缓冲数据为3G(-S 3G). 然后将制表符替换成换行符，让一行又变回到4行，后续是压缩成gz，降低磁盘占用。

这样做其实还没又解决问题，因为这两个文件里的reads数其实不相同，我们需要进一步从里面剔除掉只有一方才有的read。这里使用的seqkit，一个非常强大的序列处理轮子。

seqkit grep -f <(seqkit seq -n -i read_1.fq.gz ) read_2.fq.gz  > read_flt_2.fq
seqkit grep -f <(seqkit seq -n -i read_2.fq.gz ) read_1.fq.gz  > read_flt_1.fq
wc -l read_flt_1.fq
   15028 read_flt_1.fq
wc -l read_flt_2.fq
   15028 read_flt_2.fq

终于，你得到了配对存放的两个fastq文件。

通过这个案例我想表达的观点是：虽然有些时候没有现成的工具可以解决你的问题，但是你可以通过组合几个现有的工具来完成。也就是将一个大问题拆解成几个独立的小问题，然后逐个击破。

当然，你首先得熟悉一些已有的工具，比如说常用unix命令:sort, paste, tr等，和一些非常强大生信命令行工具：samtools, bedtools, seqkit等。

一个更佳方便的方法

前面只是问题的一种方法，可能还不是最好的，但是我当初最先想到的一个比较粗暴的方法。当然更佳巧妙的方法是，提取序列之后可以按照read name排序，然后提取

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools sort -n | samtools fastq -1 read_1.fq -2 read_2.fq -s singleton.fq -

是不是更好理解了。