Skip to content
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
2 changes: 1 addition & 1 deletion src/main/scala/com/fulcrumgenomics/bam/ClipBam.scala
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ class ClipBam
@arg(flag='m', doc="Optional output of clipping metrics.") val metrics: Option[FilePath] = None,
@arg(flag='r', doc="Reference sequence fasta file.") val ref: PathToFasta,
@arg(flag='c', doc="The type of clipping to perform.") val clippingMode: ClippingMode = ClippingMode.Hard,
@arg(flag='a', doc="Automatically clip extended attributes that are the same length as bases.") val autoClipAttributes: Boolean = false,
@arg(flag='a', doc="Automatically clip extended attributes that are the same length as bases, except listed tags that are not per-base data for the read, such as RG.") val autoClipAttributes: Boolean = false,
@arg(flag='H', doc="Upgrade all existing clipping in the input to the given clipping mode prior to applying any other clipping criteria.") val upgradeClipping: Boolean = false,
@arg( doc="Require at least this number of bases to be clipped on the 5' end of R1") val readOneFivePrime: Int = 0,
@arg( doc="Require at least this number of bases to be clipped on the 3' end of R1") val readOneThreePrime: Int = 0,
Expand Down
15 changes: 14 additions & 1 deletion src/main/scala/com/fulcrumgenomics/bam/SamRecordClipper.scala
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,17 @@ object SamRecordClipper {
/** The set of tags that should be invalidated if a read undergoes clipping. */
val TagsToInvalidate: Seq[String] = Bams.AlignmentTags

/** Tags whose values are not this read's per-base data, so they are not auto-clipped even when their length equals
* the read's. Best-effort: an unlisted tag whose length matches the read's is still clipped. */
val TagsNeverAutoClipped: Set[String] = Set(
"RG", "LB", "PU", "PG", "CO", "MI",
"BC", "QT", "RX", "QX", "OX", "BZ", "CB", "CR", "CY", "UB", "UR", "UY", "BX",
"MC", "MD", "SA", "OA", "OC", "CG", "XA", "cs", "jM", "jI",
"CC", "CT", "FS", "PT", "GX", "GN",
"R2", "Q2",
"FZ", "MM", "ML", "mv", "pi", "st", "fn"
)
Comment thread
Copilot marked this conversation as resolved.

private val NoCallBase = 'N'.toByte
private val NoCallQual = 2.toByte
}
Expand All @@ -76,7 +87,8 @@ object SamRecordClipper {
* @param autoClipAttributes if true attributes that are the same length as the bases
* and qualities be automatically clipped in the same way as
* the bases and qualities, otherwise attributes are not
* touched.
* touched. Tags in [[SamRecordClipper.TagsNeverAutoClipped]]
* are never clipped.
*/
class SamRecordClipper(val mode: ClippingMode, val autoClipAttributes: Boolean) {
import SamRecordClipper._
Expand Down Expand Up @@ -491,6 +503,7 @@ class SamRecordClipper(val mode: ClippingMode, val autoClipAttributes: Boolean)
val newLength = rec.length
val oldLength = newLength + remove
rec.attributes.foreach {
case (tag, _) if TagsNeverAutoClipped.contains(tag) => ()
case (tag, s: String) if s.length == oldLength => rec(tag) = if (fromStart) s.drop(remove) else s.take(newLength)
case (tag, a: Array[_]) if a.length == oldLength => rec(tag) = if (fromStart) a.drop(remove) else a.take(newLength)
case _ => ()
Expand Down
2 changes: 1 addition & 1 deletion src/main/scala/com/fulcrumgenomics/bam/TrimPrimers.scala
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ class TrimPrimers
@arg(flag='S', doc="Match to primer locations +/- this many bases.") val slop: Int = 5,
@arg(flag='s', doc="Sort order of output BAM file (defaults to input sort order).") val sortOrder: Option[SamOrder] = None,
@arg(flag='r', doc="Optional reference fasta for recalculating NM, MD and UQ tags.") val ref: Option[PathToFasta] = None,
@arg(flag='a', doc="Automatically trim extended attributes that are the same length as bases.") val autoTrimAttributes: Boolean = false,
@arg(flag='a', doc="Automatically trim extended attributes that are the same length as bases, except listed tags that are not per-base data for the read, such as RG.") val autoTrimAttributes: Boolean = false,
@arg(doc="Trim only first of pair reads (R1s) or fragment reads, otherwise both ends of a pair.") val firstOfPair: Boolean = false

)extends FgBioTool with LazyLogging {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -535,6 +535,31 @@ class SamRecordClipperTest extends UnitSpec with OptionValues {
withAuto[String]("az") shouldBe "67890123456789012345678901234567890"
}

it should "not auto-clip tags that are not this read's per-base data, even when their length equals the read's" in {
val expected = Set(
"RG", "LB", "PU", "PG", "CO", "MI",
"BC", "QT", "RX", "QX", "OX", "BZ", "CB", "CR", "CY", "UB", "UR", "UY", "BX",
"MC", "MD", "SA", "OA", "OC", "CG", "XA", "cs", "jM", "jI",
"CC", "CT", "FS", "PT", "GX", "GN",
"R2", "Q2",
"FZ", "MM", "ML", "mv", "pi", "st", "fn"
)
SamRecordClipper.TagsNeverAutoClipped shouldBe expected

val rec = r(10, "5S10M5S")
val value = "ACGTACGTACGTACGTACGT"
val depth = Array.tabulate[Short](20)(_.toShort)
val perBase = Seq("OQ", "E2", "XB")
expected.foreach(tag => rec(tag) = value)
perBase.foreach(tag => rec(tag) = value)
rec("cd") = depth

clipper(Hard, autoClip=true).upgradeAllClipping(rec) shouldBe (5, 5)
expected.foreach(tag => rec[String](tag) shouldBe value)
perBase.foreach(tag => rec[String](tag) shouldBe value.slice(5, 15))
rec[Array[Short]]("cd") shouldBe depth.slice(5, 15)
}

it should "not convert reads that have no soft-clipping" in {
val noSoft = r(10, "55M", attrs=Map("az" -> "12345678901234567890123456789012345678901234567890"))
val hard = r(10, "5H55M10H", attrs=Map("az" -> "12345678901234567890123456789012345678901234567890"))
Expand Down
Loading