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: 
[mem_sam_pe] [mem_sam_pe] paired reads have different names: 

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>

Looks like I solved it by using PICARD: extracting the aligned reads but without sorting, I used

java -jar picard.jar SamToFastq I=<file_alnMAP.bam> FASTQ=<filemap_1.fq> SECOND_END_FASTQ=<filemap_2.fq>
bwa mem -R <read_group> <ref2.fa> <filemap_1.fq> <filemap_1.fq>

From samtools flagstat the output looks fine.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.