I am carrying out some metagenomic comparison using whole genome sequencing approach. I would like to filter out reads from my files in a sequential manner, thus I would like to align to a first reference genome, then use the generated SAM/BAM file to extract the mapped reads with samtols view -F4, re-generate fastq files and align those to another reference genome. This decreases the size of the input and speeds the procedure. Using BWA, though, I get an error of file synchronization when doing the second alignment:
[mem_sam_pe] paired reads have different names:
"HISEQ:90:C3UNJACXX:1:1107:16388:61141",
"HISEQ:90:C3UNJACXX:3:1303:2008:58095"
[mem_sam_pe] [mem_sam_pe] paired reads have different names:
"HWI-ST1437:64:C3UM1ACXX:4:1214:3833:62731",
"HWI-ST1437:64:C3UM1ACXX:4:1202:18519:45906"
This is weird since the fastq files were generated by the filtered BAM file, samtools should have known how to keep the synchronization of the files. My question is: is it possible to realign a BAM file? maybe not with BWA, at least with another tool?
Thank you.
The procedure I carried out is:
# align
bwa mem -R <read_group1> <ref.fa> <file1_1.fq> <file1_2.fq> -o <file_aln.sam>
# convert
samtools view -Sb <file_aln.sam> > <file_aln.bam> # sort samtools sort <file_aln.bam> -o <file_alnSRT.bam>
# select mapped
samtools view -h -F 4 <file_alnSRT.bam> -o <file_alnMAP.bam>
# convert to fastq
samtools fastq -1 <filemap_1.fq.gz> -2 <filemap_1.fq.gz> <file_alnMAP.bam>
# re-align
bwa mem -R <read_group> <ref2.fa> <filemap_1.fq.gz> <filemap_1.fq.gz> | samtools sort -o <file_realign.bam>