Bugfix for BQSR: Offset into qualityScore list was wrong

[arahuja]

The wrong offset was being used to retrieve the quality score for a read. The original code filtered the BaseCovar iterator to only process high quality bases and all of the "offset" value was to the first high quality base, but then was used to look up the qualityScore in the full quality score sequence. This change

1. Renames and uses proper indices
2. Tests for ReadCovariates
3. Removed RDD arguments to covar classes as they were unsused

[AmplabJenkins]

All automated tests passed.
Refer to this link for build results: https://amplab.cs.berkeley.edu/jenkins/job/ADAM-prb/48/

[massie]

Curious about the min quality value here. Is it arbitrary or a commonly used filter value?

[arahuja]

I actually don't know the significance of using 2. This was another reason I was getting different read counts vs GATK as their default min base quality is 6. At least this way its parameterizable, but if there is a more justified default we should probably change it to that. Or add a command line argument for it.

On that note, the argument options for Transform are growing large - is the goal to keep all of the operations in there (and all their various configurations) or to start splitting them out?

[fnothaft]

Phred 2 is P=0.5.

[massie]

I think Phred 2 is P=0.63

https://www.google.com/search?q=10+%5E+(-2%2F10)&oq=10+%5E+(-2%2F10)&aqs=chrome..69i57j0j69i64l2.7131j0j7&sourceid=chrome&espv=210&es_sm=119&ie=UTF-8)

[fnothaft]

Ah—you're right—Phred 3 is P=0.5... My apologies for the mistake!

[chartl]

Q2 is a legacy value; it's effectively Illumina's garbage value; they tend
to be bad bases even after recalibration. Q2 is what the machine produces
when it's read through a nucleotide repeat and has no real idea what the
bases are. See below: all the faint bases are Q2.
[image: Inline image 1]

On Mon, Jan 27, 2014 at 1:54 PM, Matt Massie notifications@git.luolix.topwrote:

In
adam-core/src/main/scala/edu/berkeley/cs/amplab/adam/rdd/recalibration/ReadCovariates.scala:

@@ -28,27 +28,35 @@ object ReadCovariates {
}

class ReadCovariates(val read: RichADAMRecord, qualByRG: QualByRG, covars: List[StandardCovariate],

•                 val dbSNP: SnpTable) extends Iterator[BaseCovariates] with Serializable {
  

•                 val dbSNP: SnpTable, val minQuality:Int = 2) extends Iterator[BaseCovariates] with Serializable {
  

I think Phred 2 is P=0.63

https://www.google.com/search?q=10+%5E+(-2%2F10)&oq=10+%5E+(-2%2F10)&aqs=chrome..69i57j0j69i64l2.7131j0j7&sourceid=chrome&espv=210&es_sm=119&ie=UTF-8
)

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/60/files#r9200596
.

[chartl]

These are also one of the primary targets for aligners that clip
low-quality tails.

On Mon, Jan 27, 2014 at 2:01 PM, christopherlhartl@gmail.com <
christopherlhartl@gmail.com> wrote:

Q2 is a legacy value; it's effectively Illumina's garbage value; they tend
to be bad bases even after recalibration. Q2 is what the machine produces
when it's read through a nucleotide repeat and has no real idea what the
bases are. See below: all the faint bases are Q2.
[image: Inline image 1]

On Mon, Jan 27, 2014 at 1:54 PM, Matt Massie notifications@git.luolix.topwrote:

In
adam-core/src/main/scala/edu/berkeley/cs/amplab/adam/rdd/recalibration/ReadCovariates.scala:

@@ -28,27 +28,35 @@ object ReadCovariates {
}

class ReadCovariates(val read: RichADAMRecord, qualByRG: QualByRG, covars: List[StandardCovariate],

•                 val dbSNP: SnpTable) extends Iterator[BaseCovariates] with Serializable {
  

•                 val dbSNP: SnpTable, val minQuality:Int = 2) extends Iterator[BaseCovariates] with Serializable {
  

I think Phred 2 is P=0.63

https://www.google.com/search?q=10+%5E+(-2%2F10)&oq=10+%5E+(-2%2F10)&aqs=chrome..69i57j0j69i64l2.7131j0j7&sourceid=chrome&espv=210&es_sm=119&ie=UTF-8
)

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/60/files#r9200596
.

[massie]

The 'transform' commandline growing large is ok for now. We want all the transformations to be part of a single job (to minimize i/o, improve pipelining, etc). In the future, we'll probably need to use a configuration file instead which defined the exact pipeline of operations and options the user wants.

[massie]

This looks good to me. I'll wait to hear from @jey before I merge it. Jey and I just had a code walkthrough last Friday and he found this very issue. You beat us to the punch fixing it, Arun. Kudos!

[massie]

Thanks, Arun!