Samtools filtering question

[Phantom123]

I used Samtools and bcftools to generate bcf files and then vcf files from the original BAM file containing RNA-seq data. Some of the publications i have read which make use of Samtools talk of filtering based on read quality, phred score, and depth. What are the typical parameters you use for these? The commands are as follows:

samtools view -bh -q(min_MapQ) in.bam > out.bam

samtools mpileup -uf in.ref.fa –Q(min_phred_Score) in.bam > out.bam

bcftools view var.raw.bcf | vcfutils.pl varFilter –D(max_depth) > var.flt.vcf

I did the alignment once using Tophat and once using BWA so that i have two sets of BAM files for the same RNA-seq data. Do the Samtools filtering parameters depend on the alignment algorithm used?

[Heisman]

There are no hard and fast rules when it comes to filtering, typically. In terms of depth, the idea is that all else being equal, the more reads you have supporting a call, the more likely the call is to be correct. If you have a homozygous position, all reads should support it if there are no errors anywhere. The reason more read depth can be really useful is because positions can be heterozygous, for example, if you have 6 reads covering a heterozygous site and all 6 match the wildtype, you may call it homozygous wildtype. With 12 reads, it's much less likely to have all 12 match the wildtype for a heterozygous site.

Read quality and phred score of the call are somewhat intertwined. You can have a fairly high read depth but have many of the reads align with a poor mapping quality score; this implies a lot of those reads may not in fact come from the region they were aligned too (or perhaps they do). In that case, your read depth is artificially high, in a sense.

I have not analyzed RNA-Seq data so I do not want to pretend I know the normal metrics used for filtering. I do regularly call variants from sequencing data. For that, I use the varFilter as you showed. Other than that, the only thing I consider regarding filtering is coverage. I feel 10x is generally an acceptable cutoff to use, with higher or lower being fine depending on the situation.

[Richard Finney]

I did the [rna seq] alignment once using Tophat and once using BWA

What was your reference genome that you aligned to using BWA?

[Phantom123]

I did the [rna seq] alignment once using Tophat and once using BWA

What was your reference genome that you aligned to using BWA?

UCSC hg 19.

[Phantom123]

There are no hard and fast rules when it comes to filtering, typically. In terms of depth, the idea is that all else being equal, the more reads you have supporting a call, the more likely the call is to be correct. If you have a homozygous position, all reads should support it if there are no errors anywhere. The reason more read depth can be really useful is because positions can be heterozygous, for example, if you have 6 reads covering a heterozygous site and all 6 match the wildtype, you may call it homozygous wildtype. With 12 reads, it's much less likely to have all 12 match the wildtype for a heterozygous site.

Read quality and phred score of the call are somewhat intertwined. You can have a fairly high read depth but have many of the reads align with a poor mapping quality score; this implies a lot of those reads may not in fact come from the region they were aligned too (or perhaps they do). In that case, your read depth is artificially high, in a sense.

I have not analyzed RNA-Seq data so I do not want to pretend I know the normal metrics used for filtering. I do regularly call variants from sequencing data. For that, I use the varFilter as you showed. Other than that, the only thing I consider regarding filtering is coverage. I feel 10x is generally an acceptable cutoff to use, with higher or lower being fine depending on the situation.

Thank you.