0
$\begingroup$

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>
$\endgroup$
1
$\begingroup$

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.

$\endgroup$

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .