Issues with Variant Calling from old bam files

Question

[williamsbrian5064]

Hi,

I am in the process of realigning some samples to an updated genome. For the majority of the samples, I had the original fastq files. For a couple of my samples, I was only given a bam file. So what I did was use Picards "SamToFastq" function to convert the bam file back to a fastq file. Since I didn't know how the data was trimmed before I received it, I ran it through trimmomatic (ended up being a waste of time), aligned with bwa, and sorted/marked duplicated with Picard.

I then used samstool for variant calling. All of the samples that were originally bam files, weren't get any of their variants called when I looked at the vcf file. So I don't know if I messed something up after I converted the bam files back to fastq files. When I compare I compare the vcf to a newly generate bam file, you can see that there are SNPs and what not in the bam file but they aren't being called by samstool. If it was happening to other samples then I would think that I am having an issue with samtools, but this is only happening to the samples that were originally bam files that were converted back to fastq files. Am I doing something wrong here? I'm just trying to figure out what I did wrong here

Answer 1

[Devon Ryan]

Have a look at the flags of the reads, it's likely that mates are flipped the wrong way or something.

[williamsbrian5064]

hmm.. how do you do that?

[Pierre Lindenbaum]

show us the output of samtools flagstat input.bam for example.

[williamsbrian5064]

551741879 + 0 in total (QC-passed reads + QC-failed reads)
2724021 + 0 secondary
0 + 0 supplementary
59155766 + 0 duplicates
549580515 + 0 mapped (99.61% : N/A)
549017858 + 0 paired in sequencing
274508929 + 0 read1
274508929 + 0 read2
0 + 0 properly paired (0.00% : N/A)
546856494 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
21832898 + 0 with mate mapped to a different chr
5313 + 0 with mate mapped to a different chr (mapQ>=5)

This is what I got back.. looks like they aren't properly paired? Am I interpreting that correctly?

[Pierre Lindenbaum]

looks like they aren't properly paired? Am I interpreting that correctly?

yes Devon is right. Check you haven't switched reads R1 and reads R2.

[williamsbrian5064]

Wooooooooooow. The problem came up when I did that pointless trimming. I made an error and trimmed R1 twice instead of R1 and R2. Thanks for the help everyone. I really really appreciate your help.

[WouterDeCoster]

I have moved the original comment from Devon which hinted you to the correct solution to an answer, as such you can mark it as accepted and mark this question as solved.