Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Command line tool to produce and print pairwise alignments. #739

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions src/main/scala/com/fulcrumgenomics/alignment/Alignment.scala
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
package com.fulcrumgenomics.alignment

import com.fulcrumgenomics.FgBioDef._
import com.fulcrumgenomics.util.Sequences
import htsjdk.samtools.{CigarElement, CigarOperator, SAMRecord, Cigar => HtsJdkCigar}

import scala.collection.mutable.ArrayBuffer
Expand Down Expand Up @@ -484,4 +485,18 @@ case class Alignment(query: Array[Byte],

copy(queryStart=qStart, targetStart=tStart, cigar=Cigar(elems.result()), score=0)
}

/** Returns a version of the alignment where both sequences have been reverse complemented
* and the alignment has been reversed to match.
*/
def revcomped: Alignment = {
val query = this.query.clone()
val target = this.target.clone()
Sequences.revcomp(query)
Sequences.revcomp(target)
val cigar = this.cigar.reverse
val qs = query.length - (this.queryStart + this.cigar.lengthOnQuery - 1) + 1
val ts = target.length - (this.targetStart + this.cigar.lengthOnTarget - 1) + 1
Alignment(query, target, qs, ts, cigar, this.score)
}
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
/*
* The MIT License
*
* Copyright (c) 2021 Fulcrum Genomics LLC
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
* THE SOFTWARE.
*/

package com.fulcrumgenomics.personal.tfenne

import com.fulcrumgenomics.FgBioDef._
import com.fulcrumgenomics.alignment.{Aligner, Mode}
import com.fulcrumgenomics.cmdline.{ClpGroups, FgBioTool}
import com.fulcrumgenomics.fastq.{FastqSource, FastqWriter}
import com.fulcrumgenomics.sopt.{arg, clp}
import com.fulcrumgenomics.util.{Io, Sequences}

@clp(group=ClpGroups.Personal, description=
"""
|Align one query sequence to one target sequence.
""")
class PairwiseAlign
( @arg(flag='q', doc="Query sequence.") val query: String,
@arg(flag='t', doc="Target sequence.") val target: String,
@arg(flag='m', doc="Alignment mode.") val mode: Mode = Mode.Global,
) extends FgBioTool {

override def execute(): Unit = {
val aligner = Aligner(1, -4, -6, -1, mode)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

scoring params are hard-coded. Nice to have them exposed on the command line.


val left = aligner.align(query, target)
val right = aligner.align(Sequences.revcomp(query), Sequences.revcomp(target)).revcomped

println("Left Aligned:")
left.paddedString().foreach(println)

println()
println("Right Aligned:")
right.paddedString().foreach(println)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice to have the alignment score, cigar, target start/end, query start/end printed too. The latter are useful for local and glocal.

}
}
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,28 @@ class AlignmentTest extends UnitSpec {
alignment.subByTarget(6, 10).cigar.toString() shouldBe "5D"
alignment.subByTarget(6, 11).cigar.toString() shouldBe "5D1="
}

"Alignment.revcomped" should "reverse complement a simple alignment" in {
val in = Alignment("TTTTT", "TTTTT", 1, 1, Cigar("5M"), 5)
val rc = in.revcomped
rc.query.map(_.toChar).mkString shouldBe "AAAAA"
rc.target.map(_.toChar).mkString shouldBe "AAAAA"
rc.queryStart shouldBe 1
rc.targetStart shouldBe 1
rc.cigar.toString shouldBe "5M"
rc.score shouldBe 5
}

it should "correctly reverse a non-global alignment" in {
val in = Alignment("..GTAACCGGTT".filter(_ != '.'), "ACGTAA.CGGTTCCCCC".filter(_ != '.'), 1, 3, Cigar("5M1I4M"), 10)
val rc = in.revcomped
rc.query.map(_.toChar).mkString shouldBe "AACCGGTTAC"
rc.target.map(_.toChar).mkString shouldBe "GGGGGAACCGTTACGT"
rc.queryStart shouldBe 1
rc.targetStart shouldBe 6
rc.cigar.toString shouldBe "4M1I5M"
rc.score shouldBe 10
}
}

class CigarStatsTest extends UnitSpec {
Expand Down