Samtools SORT or bedtools bamToFastq giving error when converting reads to fastq.
I am trying to convert mapped reads and unmapped reads into fastq format but I am having trouble with the last bit of my pipeline.
Here is my pipeline:
bwa index refgenome.fa
bwa mem -B 2 -M refgenome.fa cert1.R1.fastq cert1.R2.fastq > cert1_aln.sam
samtools view -Sb cert1_aln.sam > cert1_aln.bam
# I do a sort here but maybe I don't need to? I have just seen many pipelines that sort and index at this point so I did it too.
samtools sort cert1_aln.bam > cert1_aln_sorted.bam | rm -f cert1_aln.bam
samtools index cert1_aln_sorted.bam
# Gather some stats
samtools flagstat cert1_aln_sorted.bam > cert1_stats.txt
# Extract unmapped read whose mate is mapped
samtools view -b -f 4 -F 264 cert1_aln_sorted.bam > cert1_tmp1_unmapped.bam
# Extract mapped read whose mate is unmapped
samtools view -b -f 8 -F 260 cert1_aln_sorted.bam > cert1_tmp2_unmapped.bam
# Extract unmapped read with unmapped mate
samtools view -b -f 12 -F 256 cert1_aln_sorted.bam > cert1_tmp3_unmapped.bam
# Merge the three tmp files into 1
samtools merge cert1_unmapped.bam cert1_tmp*_unmapped.bam
# Extract mapped reads from BAM file
samtools view -b -F 12 cert1_aln_sorted.bam > cert1_mapped.bam
# Sort the BAM files by name
samtools sort -n cert1_unmapped.bam > cert1_unmapped_sorted.bam
samtools sort -n cert1_mapped.bam > cert1_mapped_sorted.bam
# Finally, convert to Fastq
bamToFastq -i cert1_unmapped_sorted.bam -fq cert1_unmapped.R1.fastq -fq2 cert1_unmapped.R2.fastq
bamToFastq -i cert1_mapped_sorted.bam -fq cert1_mapped.R1.fastq -fq2 cert1_mapped.R2.fastq
I am getting an error somewhere in between sorting and converting this is the error I am getting and it crashes my terminal every time:
someseqname is marked as paired, but its mate does not occur next to it in your BAM file. Skipping
The above error is ENDLESSSS just spams my entire terminal and does not stop. My fastq files that I do end up with are incredibly tiny. I tried sorting without the -n flag and with it, same issue each time. I also tried using the header flag with samtools view when extracting the reads (ex: samtools view -bh -F12) but no luck...why is this happening?
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Extracting supplementary reads from bam file
written by g1ang •
-
Samtools mapped and unmapped flags, what is the difference between the flags?
written by DNAngel
-
extract ALL mapped paired-end reads from BAM file
written by romain.derelle •
-
Bamtools stats and samtools flagstat produce different results and confusing my sambamba output.
written by DNAngel
-
How to run SPAdes assembler with bam files?
written by marongiu.luigi
-
Get all unmapped read and mate in bam
written by vinaykusuma •
-
Error when trying to convert BAM file to FASTQ.
written by DNAngel
-
To extract unmapped pair end reads
written by Bioinfonext
-
Extracting mapped paired end FastQ reads from BAM file using Samtools - trying to understand -f 0x0…
written by acxcrew007 •
-
mapped/unmapped mate with samtools
written by Flavia •
Not sure if it will solve the problem but can you the fastq command in samtools?
Trying it now - will see if it works.
I think the main issue here is that
samtools sortoutput tostdout, so you will need to write something like the following in order to save your sort.Then you can check whether your file is properly sorted by using
Another good sanity check before trying to convert fastq would be to verify if the number of left and right reads is the same.
Oh sorry that was a typo on my end! I did write it to a new file. Let me fix it.
using samtools flagstat, R1 and R2 are the same number. But if I use bamtools stats, they are slightly off. Not sure why that is...
You likely have secondary alignments. When you convert BAM's back to fastq be sure to select just primary alignments.
Hey, so I just used samtools fastq cert1_unmapped_sorted.bam > cert1_unmapped.fastq and that seemed to work. Same when I used my mapped reads. However I still have to go through the files and check. Is there any particular reason some protocols suggest separating paired end reads into R1 and R2 fastq files instead of just combining them all into one file? I haven't found a solid explanation for the difference.
Interleaved data (where R1/R2 reads are mixed in a single file) is not usable in most programs and thus there is a general need to separate R1/R2 reads into two separate files.
Hello DNAngel!
Questions similar to yours can already be found at:
We have closed your question to allow us to keep similar content in the same thread.
If you disagree with this please tell us why in a reply below. We'll be happy to talk about it.
Cheers!
Sorry and thanks to the people who commented ans answered, but it is most unproductive @OP to start another wall-of-text-thread with identical content. Please keep it focused in the original one.
I do disagree with it because it is now a different error/question entirely. I did continue the question in my original post but figured based on the posted question header it won't actually be asking what I want to ask now. Two different questions.
Anyways I just added my whole question/pipeline in the original question but I doubt people will realize I have a new question in that post.