I have a set of bam files that include reads with quality scores in a mix of incompatible formats. This includes Solexa, Illumina 1.5, and Sanger encoding. The reads were aligned with bwa, so as far as I understand the misspecified quality scores didn't affect the alignment and were just copied through to the final bam file.
Now I'd like to call variants with these datasets using samtools mpileup, but I am stuck because the quality scores in all input files need to have the same encoding. That is, I am aware of the -6 flag to mpileup that would work if all the samples were in Illumina 1.5 format, but it's not applicable since I have a mix.
So my question is: can I do anything better than the brute force approaches of (1) fixing the original fastq files and realigning or (2) mucking through the bam files and changing quality scores myself with e.g. pysam?
2 answers
You wouldn't need to go through pysam. Just convert the BAM to SAM then rescale the 11th column and convert back to BAM. You won't need to even sort again probably not even index, though the latter is usually not that time consuming.
I'm not positive, but you can probably use samtools merge. It'd be easy to tell if it works after looking at IGV or Tablet.
But then again, is it is that hard to write a script that parses the quality scores column and adds or subtracts 33? Here's my untested/pseudo code:
samtools view file.bam | \
perl -lane 'if('/^#/){print; next;}
@qual=split //, $F[$qualIndex]; $_=chr(ord($_)+33) for(@qual); $F[$qualIndex]=join("",@qual); print join("\t",@F)' | \
samtools -bS > file.corrected.bam
the principle is correct but it will be more than just adding 33 - they will first need to detect from the read id which format the data is in then convert this encoded value to a numerical quality value and then re-encode that into Sanger encoding.
Thanks, that's correct. I targeted a specific re-scaling to each bam. samtools merge won't change the quality scores so it doesn't help. Just a quick clarification to your example, which is basically what I did, is to output the header so that the conversion back to bam works:
samtools view -h in.bam | rescale_quals.py | samtools view -bS - > out.bam
Could you share this rescale_quals.py script? I am in a similar position. Thank!
Hi, where could I find the rescale_quals.py script? I am in the same position. Thanks a lot.
Here it is (sadly I'm too lazy to post to github or anything better...):
#!/usr/bin/python
import sys
def process(fname):
fin = open(fname)
for line in fin:
if line[0] == '@':
print line,
continue
line = line.strip().split("\t")
if len(line[10]) != len(line[9]):
print >>sys.stderr, "number of quality scores does not match number of bases: [%s] vs. [%s]" % (line[9], line[10])
line[10] = 'O' * len(line[9])
line[10] = "".join([chr(ord(c)-31) for c in line[10]])
print "\t".join(line)
fin.close()
if __name__ == "__main__":
if len(sys.argv) <= 1:
print >>sys.stderr, "usage (converts fastq 1.3/1.5 to sanger): samtools view -h in.bam | bam_rescale_quals.py - | samtools view -bS - > out.bam"
else:
if sys.argv[1] == "-":
sys.argv[1] = "/dev/stdin"
process(sys.argv[1])
Cool. Thanks a lot.
Just to clarify, this script converts the Illumina1.5+ (PHRED64+) to Sanger encoding (Phred33+), right?
Log in to answer this question.
-
Mpileup Command Line
written by Giffredo •
-
How To Sort Sam File
written by Chen Sun
-
Sam To Bam - Loss Of Data Or Just Great Compression?
written by Clare
-
Samtools View Of A Sam File Entire Chromosome
written by dfernan
-
Samtools Merge/Reheader Of Bam Files With Different @Sq Lines
written by Nathan S. Watson-Haigh
-
Length Of Contigs In Bam File Off By 1 Versus Reference Fasta
written by Tancata •
-
Denovogear Support Marked Dup And Sorted Bam Files?
written by kulkarni.aishwarya87 •
-
How To Convert Bigwig To Bam?
written by ashis.csedu •
-
Problem Running Pindel With Bam Files
written by xenophiliuslovegood •
-
Is It Necessary To Do Base Quality Score Recalibration Before Variants Detection In Gatk
written by Lds •
Maybe merging the three data sets isn't a good idea, I think that you cannot improve your variant calling because each technology has particular bias in errors and coverages. I suggest you to analyze each set independently and report the combination of high quality calls in each data set.
Point well taken, and I'll be on the lookout for that. In this instance, the mixing is more a function of data source (SRA with Sanger quality scores versus in-house sequencing with Illumina 1.5 quality scores) than technology (but certainly there are center-specific effects that might be at play).