-
Notifications
You must be signed in to change notification settings - Fork 308
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Added support for coverage calculation from alignment files
- Loading branch information
1 parent
e7e1adf
commit f15afe5
Showing
15 changed files
with
687 additions
and
60 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
85 changes: 85 additions & 0 deletions
85
adam-cli/src/main/scala/org/bdgenomics/adam/cli/Reads2Coverage.scala
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,85 @@ | ||
| /** | ||
| * Licensed to Big Data Genomics (BDG) under one | ||
| * or more contributor license agreements. See the NOTICE file | ||
| * distributed with this work for additional information | ||
| * regarding copyright ownership. The BDG licenses this file | ||
| * to you under the Apache License, Version 2.0 (the | ||
| * "License"); you may not use this file except in compliance | ||
| * with the License. You may obtain a copy of the License at | ||
| * | ||
| * http://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| */ | ||
| package org.bdgenomics.adam.cli | ||
|
|
||
| import org.apache.spark.SparkContext | ||
| import org.bdgenomics.adam.projections.AlignmentRecordField._ | ||
| import org.bdgenomics.adam.projections.Projection | ||
| import org.bdgenomics.adam.rdd.ADAMContext._ | ||
| import org.bdgenomics.adam.rdd.features.CoverageRDD | ||
| import org.bdgenomics.adam.rdd.read.{ AlignedReadRDD, AlignmentRecordRDD } | ||
| import org.bdgenomics.utils.cli._ | ||
| import org.kohsuke.args4j.{ Argument, Option => Args4jOption } | ||
|
|
||
| /** | ||
| * Reads2Coverage (accessible as the command 'reads2coverage' through the CLI) takes two arguments, | ||
| * an INPUT and OUTPUT, and calculates the number of reads from INPUT at every location in | ||
| * the file. Optional arguments are only_negative_strands, only_positive_strands and collapse. | ||
| * only_negative_strands and only_positive_strands save coverage computed from only negative and positive strands, | ||
| * respectively. Collapse specifies whether saved coverage should merge neighboring coverage with the same counts | ||
| * to one record. | ||
| */ | ||
| object Reads2Coverage extends BDGCommandCompanion { | ||
| val commandName: String = "reads2coverage" | ||
| val commandDescription: String = "Calculate the coverage from a given ADAM file" | ||
|
|
||
| def apply(cmdLine: Array[String]): BDGCommand = { | ||
| new Reads2Coverage(Args4j[Reads2CoverageArgs](cmdLine)) | ||
| } | ||
| } | ||
|
|
||
| class Reads2CoverageArgs extends Args4jBase with ParquetArgs { | ||
| @Argument(required = true, metaVar = "INPUT", usage = "The reads file to use to calculate depths", index = 0) | ||
| var inputPath: String = null | ||
| @Argument(required = true, metaVar = "OUTPUT", usage = "Location to write the coverage data in ADAM/Parquet format", index = 1) | ||
| var outputPath: String = null | ||
| @Args4jOption(required = false, name = "-collapse", usage = "Collapses neighboring coverage records " + | ||
| "of equal counts into the same record") | ||
| var collapse: Boolean = true | ||
| @Args4jOption(required = false, name = "-only_negative_strands", usage = "Compute coverage for negative strands") | ||
| var onlyNegativeStrands: Boolean = false | ||
| @Args4jOption(required = false, name = "-only_positive_strands", usage = "Compute coverage for positive strands") | ||
| var onlyPositiveStrands: Boolean = false | ||
| } | ||
|
|
||
| class Reads2Coverage(protected val args: Reads2CoverageArgs) extends BDGSparkCommand[Reads2CoverageArgs] { | ||
| val companion: BDGCommandCompanion = Reads2Coverage | ||
|
|
||
| def run(sc: SparkContext): Unit = { | ||
|
|
||
| val proj = Projection(contigName, start, end, cigar) | ||
|
|
||
| // If saving strand specific coverage, require that only one direction is specified | ||
| require(!(args.onlyNegativeStrands && args.onlyPositiveStrands), | ||
| "Cannot compute coverage for both negative and positive strands separately") | ||
|
|
||
| // load reads | ||
| val readsRdd: AlignmentRecordRDD = sc.loadAlignments(args.inputPath, projection = Some(proj)) | ||
|
|
||
| val finalReads = if (args.onlyNegativeStrands && !args.onlyPositiveStrands) { | ||
| readsRdd.transform(rdd => rdd.filter(_.getReadNegativeStrand)) | ||
| } else if (!args.onlyNegativeStrands && args.onlyPositiveStrands) { | ||
| readsRdd.transform(rdd => rdd.filter(!_.getReadNegativeStrand)) | ||
| } else { | ||
| readsRdd | ||
| } | ||
|
|
||
| finalReads.toCoverage(args.collapse) | ||
| .save(args.outputPath) | ||
| } | ||
| } |
37 changes: 37 additions & 0 deletions
37
adam-cli/src/test/scala/org/bdgenomics/adam/cli/Reads2CoverageSuite.scala
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,37 @@ | ||
| /** | ||
| * Licensed to Big Data Genomics (BDG) under one | ||
| * or more contributor license agreements. See the NOTICE file | ||
| * distributed with this work for additional information | ||
| * regarding copyright ownership. The BDG licenses this file | ||
| * to you under the Apache License, Version 2.0 (the | ||
| * "License"); you may not use this file except in compliance | ||
| * with the License. You may obtain a copy of the License at | ||
| * | ||
| * http://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| */ | ||
| package org.bdgenomics.adam.cli | ||
|
|
||
| import org.bdgenomics.adam.rdd.ADAMContext._ | ||
| import org.bdgenomics.adam.util.ADAMFunSuite | ||
| import org.bdgenomics.utils.cli.Args4j | ||
|
|
||
| class Reads2CoverageSuite extends ADAMFunSuite { | ||
|
|
||
| sparkTest("correctly calculates coverage from small sam file") { | ||
| val inputPath = copyResource("artificial.sam") | ||
| val outputPath = tmpFile("coverage.adam") | ||
|
|
||
| val args: Array[String] = Array(inputPath, outputPath) | ||
| new Reads2Coverage(Args4j[Reads2CoverageArgs](args)).run(sc) | ||
| val coverage = sc.loadCoverage(outputPath) | ||
|
|
||
| val pointCoverage = coverage.flatten.rdd.filter(_.start == 30).first | ||
| assert(pointCoverage.count == 5) | ||
| } | ||
| } |
85 changes: 85 additions & 0 deletions
85
adam-core/src/main/scala/org/bdgenomics/adam/models/Coverage.scala
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,85 @@ | ||
| /** | ||
| * Licensed to Big Data Genomics (BDG) under one | ||
| * or more contributor license agreements. See the NOTICE file | ||
| * distributed with this work for additional information | ||
| * regarding copyright ownership. The BDG licenses this file | ||
| * to you under the Apache License, Version 2.0 (the | ||
| * "License"); you may not use this file except in compliance | ||
| * with the License. You may obtain a copy of the License at | ||
| * | ||
| * http://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| */ | ||
| package org.bdgenomics.adam.models | ||
|
|
||
| import org.apache.spark.rdd.RDD | ||
| import org.bdgenomics.formats.avro.Feature | ||
|
|
||
| /** | ||
| * Converts from avro Feature to Coverage. | ||
| */ | ||
| object Coverage { | ||
|
|
||
| /** | ||
| * Creates Coverage from ReferenceRegion and coverage count in that ReferenceRegion. | ||
| * | ||
| * @param region ReferenceRegion in which Coverage spans | ||
| * @param count Coverage count for each base pair in region | ||
| * @return Coverage spanning the specified ReferenceRegion | ||
| */ | ||
| private[adam] def apply(region: ReferenceRegion, count: Double): Coverage = { | ||
| Coverage(region.referenceName, region.start, region.end, count) | ||
| } | ||
|
|
||
| /** | ||
| * Creates Coverage from Feature, extracting region information and feature score for coverage. | ||
| * | ||
| * @param feature Feature to create coverage from | ||
| * @return Coverage spanning the specified feature | ||
| */ | ||
| private[adam] def apply(feature: Feature): Coverage = { | ||
| Coverage(feature.getContigName, feature.getStart, feature.getEnd, feature.getScore) | ||
| } | ||
|
|
||
| /** | ||
| * Creates an RDD of Coverage from RDD of Features. | ||
| * | ||
| * @param rdd RDD of Features to extract Coverage from | ||
| * @return RDD of Coverage spanning all features in rdd | ||
| */ | ||
| private[adam] def apply(rdd: RDD[Feature]): RDD[Coverage] = { | ||
| rdd.map(f => Coverage(f)) | ||
| } | ||
| } | ||
|
|
||
| /** | ||
| * Coverage record for CoverageRDD. | ||
| * Contains Region indexed by contig name, start and end, as well as count of coverage at | ||
| * each base pair in that region. | ||
| * | ||
| * @param contigName Specifies chromosomal location of coverage | ||
| * @param start Specifies start position of coverage | ||
| * @param end Specifies end position of coverage | ||
| * @param count Specifies count of coverage at location | ||
| */ | ||
| case class Coverage(contigName: String, start: Long, end: Long, count: Double) { | ||
|
|
||
| /** | ||
| * Converts Coverage to Feature, setting Coverage count in the score attribute. | ||
| * | ||
| * @return Feature built from Coverage | ||
| */ | ||
| def toFeature: Feature = { | ||
| val fb = Feature.newBuilder() | ||
| fb.setContigName(contigName) | ||
| fb.setStart(start) | ||
| fb.setEnd(end) | ||
| fb.setScore(count) | ||
| fb.build() | ||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.