TARS: Transcript Alignment for RNA Splicing (AUS-418)

.githooks/pre-commit

@@ -0,0 +1,65 @@
+#!/usr/bin/env bash
+#
+# Reformats staged .java files with IntelliJ's command-line formatter using the
+# repo's hmf_code_scheme.xml, so commits match the in-IDE style.
+#
+# One-time setup per clone:   git config core.hooksPath .githooks
+# Skip a single commit:       git commit --no-verify
+# Override the formatter:      export FORMAT_SH=/path/to/format.sh
+#
+set -euo pipefail
+
+repo_root="$(git rev-parse --show-toplevel)"
+scheme="$repo_root/hmf_code_scheme.xml"
+
+# Staged Java files (added/copied/modified), NUL-safe (bash 3.2 compatible).
+java_files=()
+while IFS= read -r -d '' f; do
+  java_files+=("$f")
+done < <(git diff --cached --name-only --diff-filter=ACM -z -- '*.java')
+[[ ${#java_files[@]} -eq 0 ]] && exit 0
+
+# Resolve IntelliJ's format.sh.
+format_sh="${FORMAT_SH:-}"
+if [[ -z "$format_sh" ]]; then
+  for candidate in \
+    "/Applications/IntelliJ IDEA.app/Contents/bin/format.sh" \
+    "/Applications/IntelliJ IDEA CE.app/Contents/bin/format.sh" \
+    "$HOME/.local/share/JetBrains/Toolbox/apps/"*/*/bin/format.sh \
+    /opt/intellij*/bin/format.sh; do
+    [[ -x "$candidate" ]] && { format_sh="$candidate"; break; }
+  done
+fi
+
+if [[ -z "$format_sh" || ! -x "$format_sh" ]]; then
+  echo "pre-commit: IntelliJ format.sh not found; skipping formatting." >&2
+  echo "            Set FORMAT_SH=/path/to/format.sh to enable." >&2
+  exit 0
+fi
+
+# Absolute paths for the formatter.
+abs_files=()
+for f in "${java_files[@]}"; do
+  abs_files+=("$repo_root/$f")
+done
+
+# Run with an isolated config/system dir so the formatter doesn't collide with
+# an already-open IntelliJ ("Only one instance of IDEA can be run at a time").
+# Cached under .git for warm-start speed.
+hook_idea_dir="$(git rev-parse --git-path hmf-format-idea)"
+mkdir -p "$hook_idea_dir"
+hook_idea_dir="$(cd "$hook_idea_dir" && pwd)"
+props="$hook_idea_dir/idea.properties"
+printf 'idea.config.path=%s/config\nidea.system.path=%s/system\n' "$hook_idea_dir" "$hook_idea_dir" > "$props"
+export IDEA_PROPERTIES="$props"
+
+echo "pre-commit: formatting ${#abs_files[@]} Java file(s) with IntelliJ scheme..."
+if ! "$format_sh" -s "$scheme" "${abs_files[@]}"; then
+  echo "pre-commit: format.sh failed; commit with --no-verify to bypass." >&2
+  exit 1
+fi
+
+# Re-stage anything the formatter touched.
+for f in "${java_files[@]}"; do
+  git add -- "$repo_root/$f"
+done

hmf-common/src/main/java/com/hartwig/hmftools/common/bam/CigarUtils.java

@@ -7,10 +7,12 @@
 import static com.hartwig.hmftools.common.bam.SamRecordUtils.NO_CIGAR;
 import static com.hartwig.hmftools.common.bam.SamRecordUtils.NO_POSITION;
 
+import static htsjdk.samtools.CigarOperator.EQ;
 import static htsjdk.samtools.CigarOperator.H;
 import static htsjdk.samtools.CigarOperator.I;
 import static htsjdk.samtools.CigarOperator.M;
 import static htsjdk.samtools.CigarOperator.S;
+import static htsjdk.samtools.CigarOperator.X;
 
 import java.util.ArrayList;
 import java.util.Collections;
@@ -77,6 +79,34 @@ public static int cigarAlignedLength(final Cigar cigar)
         return cigar.getCigarElements().stream().filter(x -> x.getOperator().consumesReferenceBases()).mapToInt(x -> x.getLength()).sum();
     }
 
+    public static int cigarBaseLength(final List<CigarElement> elements)
+    {
+        return elements.stream().filter(x -> x.getOperator().consumesReadBases()).mapToInt(CigarElement::getLength).sum();
+    }
+
+    public static int cigarAlignedLength(final List<CigarElement> elements)
+    {
+        return elements.stream().filter(x -> x.getOperator().consumesReferenceBases()).mapToInt(CigarElement::getLength).sum();
+    }
+
+    // total M/=/X bases (aligned bases excluding indels and skips); matches STAR's overhang definition
+    public static int matchedBases(final List<CigarElement> elements)
+    {
+        int total = 0;
+        for(CigarElement element : elements)
+        {
+            CigarOperator operator = element.getOperator();
+            if(operator == M || operator == EQ || operator == X)
+                total += element.getLength();
+        }
+        return total;
+    }
+
+    public static boolean hasHardClip(final List<CigarElement> elements)
+    {
+        return elements.stream().anyMatch(x -> x.getOperator() == H);
+    }
+
     public static int calcCigarAlignedLength(final String cigarStr)
     {
         // a string-parsing version of the method above
@@ -142,6 +172,12 @@ public static int rightSoftClipLength(final Cigar cigar)
         return (lastElement != null && lastElement.getOperator() == S) ? lastElement.getLength() : 0;
     }
 
+    public static int rightSoftClipLength(final List<CigarElement> elements)
+    {
+        CigarElement lastElement = elements.get(elements.size() - 1);
+        return lastElement.getOperator() == S ? lastElement.getLength() : 0;
+    }
+
     public static int leftHardClipLength(final Cigar cigar)
     {
         CigarElement firstElement = cigar.getFirstCigarElement();

hmf-common/src/main/java/com/hartwig/hmftools/common/bam/SamRecordUtils.java

@@ -40,6 +40,7 @@ public final class SamRecordUtils
     public static final String READ_GROUP_ATTRIBUTE = SAMTag.RG.name();
 
     public static final String XS_ATTRIBUTE = "XS";
+    public static final String XA_ATTRIBUTE = "XA";
 
     // Redux tags
     public static final String CONSENSUS_READ_ATTRIBUTE = "CR";

hmf-common/src/main/java/com/hartwig/hmftools/common/logging/HmfLogger.java

@@ -0,0 +1,43 @@
+package com.hartwig.hmftools.common.logging;
+
+import org.apache.logging.log4j.Level;
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.spi.ExtendedLogger;
+import org.apache.logging.log4j.spi.ExtendedLoggerWrapper;
+
+// Logger with a custom DEBUG_2 level: finer than DEBUG (intLevel 500), coarser than TRACE (600), so it can be enabled to show
+// extra per-read detail without pulling in third-party TRACE noise (htsjdk and friends log at TRACE). The wrapper mirrors what
+// log4j2's own custom-level generator produces, exposing TARS_LOGGER.debug2(...) alongside the standard level methods.
+public final class HmfLogger extends ExtendedLoggerWrapper
+{
+    public static final Level DEBUG_2 = Level.forName("DEBUG_2", 550);
+
+    private static final String FQCN = HmfLogger.class.getName();
+
+    private HmfLogger(final ExtendedLogger logger)
+    {
+        super(logger, logger.getName(), logger.getMessageFactory());
+    }
+
+    public static HmfLogger getLogger(final Class<?> clazz)
+    {
+        return new HmfLogger((ExtendedLogger) LogManager.getLogger(clazz));
+    }
+
+    public static HmfLogger getLogger(final String name)
+    {
+        return new HmfLogger((ExtendedLogger) LogManager.getLogger(name));
+    }
+
+    // Resolves standard and hmftools custom level names. Referencing this class loads DEBUG_2 (via the static field above) so a
+    // later Level.valueOf("DEBUG_2") does not throw.
+    public static Level parseLevel(final String name)
+    {
+        return Level.valueOf(name);
+    }
+
+    public void debug2(final String message, final Object... params)
+    {
+        logIfEnabled(FQCN, DEBUG_2, null, message, params);
+    }
+}

hmf-common/src/main/java/com/hartwig/hmftools/common/utils/config/ConfigUtils.java

@@ -5,6 +5,8 @@
 import static com.hartwig.hmftools.common.utils.file.FileDelimiters.CSV_DELIM;
 import static com.hartwig.hmftools.common.utils.file.FileReaderUtils.createFieldsIndexMap;
 
+import com.hartwig.hmftools.common.logging.HmfLogger;
+
 import java.io.File;
 import java.io.IOException;
 import java.nio.file.Files;
@@ -26,7 +28,7 @@ public class ConfigUtils
     public static final String LOG_DEBUG_DESC = "Log at DEBUG level";
 
     public static final String LOG_LEVEL = "log_level";
-    public static final String LOG_LEVEL_DESC = "Specify log level: ERROR, WARN, INFO, DEBUG or TRACE";
+    public static final String LOG_LEVEL_DESC = "Specify log level: ERROR, WARN, INFO, DEBUG, DEBUG_2 or TRACE";
 
     public static final String SAMPLE_ID_FILE = "sample_id_file";
     public static final String SAMPLE_ID_FILE_DESC = "Sample ID CSV file with 'SampleId' column";
@@ -54,7 +56,7 @@ public static void setLogLevel(final ConfigBuilder configBuilder)
         }
         else if(configBuilder.isRegistered(LOG_LEVEL) && configBuilder.hasValue(LOG_LEVEL))
         {
-            Configurator.setRootLevel(Level.valueOf(configBuilder.getValue(LOG_LEVEL)));
+            Configurator.setRootLevel(HmfLogger.parseLevel(configBuilder.getValue(LOG_LEVEL)));
         }
     }
 

hmf-common/src/main/java/com/hartwig/hmftools/common/utils/config/LoggingOptions.java

@@ -2,6 +2,8 @@
 import org.apache.logging.log4j.core.config.Configurator;
 import org.apache.logging.log4j.Level;
 
+import com.hartwig.hmftools.common.logging.HmfLogger;
+
 import com.beust.jcommander.IStringConverter;
 import com.beust.jcommander.Parameter;
 
@@ -19,15 +21,15 @@ static class Log4jLevelTypeConverter implements IStringConverter<Level>
         @Override
         public Level convert(String value)
         {
-            return Level.valueOf(value);
+            return HmfLogger.parseLevel(value);
         }
     }
 
     @Parameter(names = "-" + ConfigUtils.LOG_DEBUG, description = "Log verbose")
     public boolean LogDebug;
 
     @Parameter(names = "-" + ConfigUtils.LOG_LEVEL,
-               description = "Specify log level (values: [OFF, FATAL, ERROR, WARN, INFO, DEBUG, TRACE, ALL])",
+               description = "Specify log level (values: [OFF, FATAL, ERROR, WARN, INFO, DEBUG, DEBUG_2, TRACE, ALL])",
                converter = Log4jLevelTypeConverter.class)
     public org.apache.logging.log4j.Level LogLevel = null;
 

hmf-common/src/test/java/com/hartwig/hmftools/common/bam/CigarUtilsTest.java

@@ -33,6 +33,20 @@
 
 public class CigarUtilsTest
 {
+    @Test
+    public void testElementListHelpers()
+    {
+        final List<CigarElement> elements = CigarUtils.cigarElementsFromStr("10S94M79N57M3S");
+
+        assertEquals(164, CigarUtils.cigarBaseLength(elements));     // read-consuming: 10S + 94M + 57M + 3S
+        assertEquals(230, CigarUtils.cigarAlignedLength(elements));  // ref-consuming: 94M + 79N + 57M
+        assertEquals(151, CigarUtils.matchedBases(elements));        // M/=/X only: 94M + 57M
+        assertEquals(10, CigarUtils.leftSoftClipLength(elements));
+        assertEquals(3, CigarUtils.rightSoftClipLength(elements));
+        assertFalse(CigarUtils.hasHardClip(elements));
+        assertTrue(CigarUtils.hasHardClip(CigarUtils.cigarElementsFromStr("5H94M")));
+    }
+
     @Test
     public void testCigarFromStr()
     {

hmf-common/src/test/java/com/hartwig/hmftools/common/logging/HmfLoggerTest.java

@@ -0,0 +1,40 @@
+package com.hartwig.hmftools.common.logging;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+
+import org.apache.logging.log4j.Level;
+import org.junit.Test;
+
+public class HmfLoggerTest
+{
+    @Test
+    public void debug2RegisteredAndParseable()
+    {
+        assertEquals(HmfLogger.DEBUG_2, HmfLogger.parseLevel("DEBUG_2"));
+        assertEquals(550, HmfLogger.DEBUG_2.intLevel());
+    }
+
+    @Test
+    public void parseLevelResolvesStandardLevels()
+    {
+        assertEquals(Level.INFO, HmfLogger.parseLevel("INFO"));
+        assertEquals(Level.TRACE, HmfLogger.parseLevel("TRACE"));
+    }
+
+    @Test
+    public void debug2SitsBetweenDebugAndTrace()
+    {
+        // finer than DEBUG (shows everything DEBUG does), but coarser than TRACE (excludes third-party TRACE noise)
+        assertTrue(HmfLogger.DEBUG_2.isLessSpecificThan(Level.DEBUG));
+        assertTrue(HmfLogger.DEBUG_2.isMoreSpecificThan(Level.TRACE));
+    }
+
+    @Test
+    public void debug2MethodLogsWithoutError()
+    {
+        // smoke: the extended-logger method resolves and runs with parameterized placeholders
+        HmfLogger logger = HmfLogger.getLogger(HmfLoggerTest.class);
+        logger.debug2("debug2 smoke {} {}", "a", 1);
+    }
+}

pom.xml

@@ -42,6 +42,7 @@
         <module>sage</module>
         <module>sigs</module>
         <module>sv-tools</module>
+        <module>tars</module>
         <module>teal</module>
         <module>virus-interpreter</module>
         <module>v-chord</module>
@@ -95,6 +96,7 @@
         <sage.version>5.1</sage.version>
         <sigs.version>1.2.1</sigs.version>
         <sv-tools.version>1.0</sv-tools.version>
+        <tars.version>1.0</tars.version>
         <teal.version>1.4</teal.version>
         <virus-interpreter.version>1.7.1</virus-interpreter.version>
         <v-chord.version>1.0.1</v-chord.version>

redux/src/main/java/com/hartwig/hmftools/redux/duplicate/ReadCache.java

@@ -55,6 +55,8 @@ public class ReadCache
     public static final int DEFAULT_POP_DISTANCE_CHECK = 100; // how often in base terms to check for popping read groups
     public static final int DEFAULT_LOG_READ_COUNT_THRESHOLD = 100000; // based on observed cache sizes for deep panels
     public static final int DEFAULT_DYNAMIC_READ_COUNT_THRESHOLD = 100_000; // level at which steps are taken to reduce the cache size
+    // log a single fragment-coord group reaching this many reads (extreme-expression locus)
+    public static final int LARGE_FRAG_GROUP_LOG_SIZE = 1000;
 
     public ReadCache(
             int groupSize, int maxSoftClipLength, boolean useFragmentOrientation, final DuplicatesConfig duplicatesConfig,
@@ -465,6 +467,12 @@ public void addRead(final FragmentCoords fragCoords, final SAMRecord read)
             }
 
             reads.add(read);
+
+            if(reads.size() == LARGE_FRAG_GROUP_LOG_SIZE)
+            {
+                RD_LOGGER.debug("large frag-coord group({} reads) at {}:{} cigar({})",
+                        reads.size(), Chromosome, read.getAlignmentStart(), read.getCigarString());
+            }
         }
 
         public int readCount() { return FragCoordsMap.values().stream().mapToInt(x -> x.size()).sum(); }

redux/src/main/java/com/hartwig/hmftools/redux/write/BamWriterNoSync.java

@@ -105,9 +105,10 @@ public void close()
         {
             mSortedBamWriter.flush();
 
-            RD_LOGGER.debug("sorted-writer records written({} writes={} avg={} max={}) maxCache({}) to BAM({})",
+            RD_LOGGER.debug("sorted-writer records written({} writes={} avg={} max={}) maxCache({} at {}) to BAM({})",
                     mSortedBamWriter.written(), mSortedBamWriter.writeCount(), mSortedBamWriter.avgWriteCount(),
-                    mSortedBamWriter.maxWrite(), mSortedBamWriter.maxCache(), filenamePart(mFilename));
+                    mSortedBamWriter.maxWrite(), mSortedBamWriter.maxCache(), mSortedBamWriter.maxCacheLocation(),
+                    filenamePart(mFilename));
         }
         else if(mBamWriter != null)
         {