ISOFOX: integrate tars support (AUS-418)

compar/src/main/java/com/hartwig/hmftools/compar/isofox/IsofoxSummaryComparer.java

@@ -6,7 +6,6 @@
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_ALT_FRAG_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_CHIMERIC_FRAG_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_DUPLICATE_FRAGS;
-import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_ENRICHED_GENE_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FORWARD_STRAND_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FRAG_LENGTH_50TH;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FRAG_LENGTH_5TH;
@@ -71,7 +70,6 @@ public void registerThresholds(final DiffThresholds thresholds)
         thresholds.addFieldThreshold(FLD_FRAG_LENGTH_5TH, -1, 0.05);
         thresholds.addFieldThreshold(FLD_FRAG_LENGTH_50TH, -1, 0.05);
         thresholds.addFieldThreshold(FLD_FRAG_LENGTH_95TH, -1, 0.05);
-        thresholds.addFieldThreshold(FLD_ENRICHED_GENE_PERC, 0.01, -1);
         thresholds.addFieldThreshold(FLD_MEDIAN_GC_RATIO, 0.01, -1);
         thresholds.addFieldThreshold(FLD_FORWARD_STRAND_PERC, 0.01, -1);
     }
@@ -81,8 +79,7 @@ public List<String> comparedFieldNames()
     {
         return List.of(
                 FLD_QC_STATUS, FLD_TOTAL_FRAGS, FLD_DUPLICATE_FRAGS, FLD_SPLICED_FRAG_PERC, FLD_UNSPLICED_FRAG_PERC, FLD_ALT_FRAG_PERC,
-                FLD_CHIMERIC_FRAG_PERC, FLD_READ_LENGTH, FLD_FRAG_LENGTH_5TH, FLD_FRAG_LENGTH_50TH, FLD_FRAG_LENGTH_95TH,
-                FLD_ENRICHED_GENE_PERC
+                FLD_CHIMERIC_FRAG_PERC, FLD_READ_LENGTH, FLD_FRAG_LENGTH_5TH, FLD_FRAG_LENGTH_50TH, FLD_FRAG_LENGTH_95TH
         );
     }
 

compar/src/main/java/com/hartwig/hmftools/compar/isofox/IsofoxSummaryData.java

@@ -32,7 +32,6 @@ public record IsofoxSummaryData(RnaStatistics RnaStatistics) implements Comparab
     static final String FLD_FRAG_LENGTH_5TH = Column.FragLength5th.toString();
     static final String FLD_FRAG_LENGTH_50TH = Column.FragLength50th.toString();
     static final String FLD_FRAG_LENGTH_95TH = Column.FragLength95th.toString();
-    static final String FLD_ENRICHED_GENE_PERC = Column.EnrichedGenePercent.toString();
     static final String FLD_MEDIAN_GC_RATIO = Column.MedianGCRatio.toString();
     static final String FLD_FORWARD_STRAND_PERC = Column.ForwardStrandPercent.toString();
 
@@ -63,7 +62,6 @@ public List<String> displayValues()
         values.add(format("%.1f", RnaStatistics.fragmentLength5thPercent()));
         values.add(format("%.1f", RnaStatistics.fragmentLength50thPercent()));
         values.add(format("%.1f", RnaStatistics.fragmentLength95thPercent()));
-        values.add(format("%.2f", RnaStatistics.enrichedGenePercent()));
         return values;
     }
 
@@ -101,7 +99,6 @@ public Mismatch findMismatch(
         checkDiff(diffs, FLD_FRAG_LENGTH_5TH, ref.fragmentLength5thPercent(), otherData.fragmentLength5thPercent(), thresholds);
         checkDiff(diffs, FLD_FRAG_LENGTH_50TH, ref.fragmentLength50thPercent(), otherData.fragmentLength50thPercent(), thresholds);
         checkDiff(diffs, FLD_FRAG_LENGTH_95TH, ref.fragmentLength95thPercent(), otherData.fragmentLength95thPercent(), thresholds);
-        checkDiff(diffs, FLD_ENRICHED_GENE_PERC, ref.enrichedGenePercent(), otherData.enrichedGenePercent(), thresholds);
         checkDiff(diffs, FLD_MEDIAN_GC_RATIO, ref.medianGCRatio(), otherData.medianGCRatio(), thresholds);
         checkDiff(diffs, FLD_FORWARD_STRAND_PERC, ref.forwardStrandPercent(), otherData.forwardStrandPercent(), thresholds);
 

compar/src/test/java/com/hartwig/hmftools/compar/isofox/IsofoxSummaryDataTest.java

@@ -3,7 +3,6 @@
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_ALT_FRAG_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_CHIMERIC_FRAG_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_DUPLICATE_FRAGS;
-import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_ENRICHED_GENE_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FORWARD_STRAND_PERC;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FRAG_LENGTH_50TH;
 import static com.hartwig.hmftools.compar.isofox.IsofoxSummaryData.FLD_FRAG_LENGTH_5TH;
@@ -45,7 +44,6 @@ public void setUp()
                 Map.entry(FLD_FRAG_LENGTH_5TH, b -> b.fragmentLength5thPercent = alternateValueSource.RnaStatistics().fragmentLength5thPercent()),
                 Map.entry(FLD_FRAG_LENGTH_50TH, b -> b.fragmentLength50thPercent = alternateValueSource.RnaStatistics().fragmentLength50thPercent()),
                 Map.entry(FLD_FRAG_LENGTH_95TH, b -> b.fragmentLength95thPercent = alternateValueSource.RnaStatistics().fragmentLength95thPercent()),
-                Map.entry(FLD_ENRICHED_GENE_PERC, b -> b.enrichedGenePercent = alternateValueSource.RnaStatistics().enrichedGenePercent()),
                 Map.entry(FLD_MEDIAN_GC_RATIO, b -> b.medianGCRatio = alternateValueSource.RnaStatistics().medianGCRatio()),
                 Map.entry(FLD_FORWARD_STRAND_PERC, b -> b.forwardStrandPercent = alternateValueSource.RnaStatistics().forwardStrandPercent())
         );

compar/src/test/java/com/hartwig/hmftools/compar/isofox/TestIsofoxSummaryDataBuilder.java

@@ -21,7 +21,6 @@ public class TestIsofoxSummaryDataBuilder
     public double fragmentLength5thPercent = 31;
     public double fragmentLength50thPercent = 100;
     public double fragmentLength95thPercent = 150;
-    public double enrichedGenePercent = 0.2;
     public double medianGCRatio = 0.5;
     public double forwardStrandPercent = 0.9;
 
@@ -38,7 +37,6 @@ public class TestIsofoxSummaryDataBuilder
         b.fragmentLength5thPercent = 25;
         b.fragmentLength50thPercent = 80;
         b.fragmentLength95thPercent = 120;
-        b.enrichedGenePercent = 0.4;
         b.medianGCRatio = 0.6;
         b.forwardStrandPercent = 0.5;
     };
@@ -61,7 +59,6 @@ private IsofoxSummaryData build()
                 .fragmentLength5thPercent(fragmentLength5thPercent)
                 .fragmentLength50thPercent(fragmentLength50thPercent)
                 .fragmentLength95thPercent(fragmentLength95thPercent)
-                .enrichedGenePercent(enrichedGenePercent)
                 .medianGCRatio(medianGCRatio)
                 .forwardStrandPercent(forwardStrandPercent)
                 .build();

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/rna/RnaStatisticFile.java

@@ -50,7 +50,6 @@ public enum Column
         FragLength5th,
         FragLength50th,
         FragLength95th,
-        EnrichedGenePercent,
         MedianGCRatio,
         ForwardStrandPercent;
     }
@@ -84,7 +83,6 @@ public static String writeLine(final String sampleId, final RnaStatistics statis
         sj.add(String.format("%.0f", statistics.fragmentLength5thPercent()));
         sj.add(String.format("%.0f", statistics.fragmentLength50thPercent()));
         sj.add(String.format("%.0f", statistics.fragmentLength95thPercent()));
-        sj.add(String.format("%.3f", statistics.enrichedGenePercent()));
         sj.add(String.format("%.3f", statistics.medianGCRatio()));
         sj.add(String.format("%.3f", statistics.forwardStrandPercent()));
         return sj.toString();
@@ -128,7 +126,6 @@ public static RnaStatistics fromLines(final List<String> lines)
                 .fragmentLength5thPercent(getDoubleValue(fieldsIndexMap, Column.FragLength5th.toString(), values))
                 .fragmentLength50thPercent(getDoubleValue(fieldsIndexMap, Column.FragLength50th.toString(), values))
                 .fragmentLength95thPercent(getDoubleValue(fieldsIndexMap, Column.FragLength95th.toString(), values))
-                .enrichedGenePercent(getDoubleValue(fieldsIndexMap, Column.EnrichedGenePercent.toString(), values))
                 .medianGCRatio(getDoubleValue(fieldsIndexMap, Column.MedianGCRatio.toString(), values))
                 .forwardStrandPercent(getDoubleValue(fieldsIndexMap, Column.ForwardStrandPercent.toString(), values))
                 .build();

hmf-common/src/main/java/com/hartwig/hmftools/common/rna/RnaStatistics.java

@@ -22,10 +22,6 @@ public abstract class RnaStatistics
     public abstract double fragmentLength50thPercent();
     public abstract double fragmentLength95thPercent();
 
-    // proportion of fragments in 7 highly expressed genes
-    public abstract double enrichedGenePercent();
-
-    // Median GC (excluding 7 highly expressed genes)
     public abstract double medianGCRatio();
 
     public abstract double forwardStrandPercent();

isofox/README.md

@@ -11,7 +11,7 @@ For transcript abundance, Isofox uses a similar methodology to several previous
 * Avoids overfitting of 'retained intron' transcripts which may simply be intronic reads 
 * Individual or combinations of splice junctions which are unique to a transcript will be weighed strongly. Does not overfit variability of coverage within exons
 
-The input for Isofox is mapped paired end reads (we use STAR for our aligner).
+The input for Isofox is mapped paired end reads. We align with bwa-mem2 against a transcriptome-augmented reference and lift the alignments back to genomic coordinates with tars, then mark duplicates with redux; Isofox takes the resulting post-tars, post-redux BAM.
 
 
 ### A note on duplicates, highly expressed genes, raw and adjusted TPM
@@ -22,19 +22,9 @@ We find that 6 genes in particular (RN7SL2, RN7SL1, RN7SL3, RN7SL4P, RN7SL5P & R
 In addition, any junction which maps in the Poly-G region of LINC00486 is filtered from all analyses (v38: chr2:32,916,190-32,916,630; v37: 2:33,141,260-33,141,700) as they are likely the result of Poly-G sequencer artefacts.
 
 ### A note on alignment and multi-mapping
-We use STAR as our aligner. ISOFOX expects BAM output with chimeric reads in the BAM itself, so it is essential when using STAR to set the outSAMtype to 'BAM Unsorted' and the chimOutType to 'WithinBAM'
+Reads are aligned with bwa-mem2 against a transcriptome-augmented reference and lifted back to genomic coordinates by tars, then duplicate-marked by redux. Chimeric and supplementary alignments are retained in the BAM.
 
-The full list of non default parameters we use internally is:
-
-```
---outSAMtype BAM Unsorted --outSAMunmapped Within --outBAMcompression 0 --outSAMattributes All --outFilterMultimapNmax 10 
---outFilterMismatchNmax 3 limitOutSJcollapsed 3000000 -chimSegmentMin 10 --chimOutType WithinBAM SoftClip 
---chimJunctionOverhangMin 10 --chimSegmentReadGapMax 3 --chimScoreMin 1 --chimScoreDropMax 30 --chimScoreJunctionNonGTAG 0 
---chimScoreSeparation 1 --outFilterScoreMinOverLread 0.33 --outFilterMatchNminOverLread 0.33 --outFilterMatchNmin 35 
---alignSplicedMateMapLminOverLmate 0.33 --alignSplicedMateMapLmin 35 --alignSJstitchMismatchNmax 5 -1 5 5
-```
-
-STAR allows setting of the `—outFilterMultimapNmax` parameter to specify the maximum number of multimaps to allow for an alignment, and we use the default value (10).   STAR will mark one of the multi-mappable reads as primary and the remainder as secondary reads. By default, STAR sets MAPQ for mappable reads to 255, whereas multi-mappable reads will have qual scores of 3 or less. For transcript abundance only, Isofox counts both primary and secondary reads at all locations, but reduces the weight of the reads to reflect the multi-mapping  (MAPQ 3 = 50%,  MAPQ 2 = 33%, MAPQ 1 = 20%, MAPQ 0 => 10%).    Reads with MAPQ of <10 are excluded from novel splice junction and chimeric analysis
+Isofox supports both bwa-tars and STAR alignments, selected by `-aligner` (`bwa-tars` is the default, or `star`); the flag only affects how multi-mapped fragments are handled, which is the one place the two aligners differ. Under `bwa-tars` a multi-mapped read is a single primary alignment carrying its alternate loci in the bwa `XA` tag (no secondary records; map qualities 0 to 60, with a confident single-locus read at 60), and the fragment is counted once at its primary locus and flagged multi-mapped. Under `star` the alternate mappings are separate secondary records and ambiguity is encoded in the map quality (255 unique, 3 or lower multi-mapped); a multi-mapped fragment is down-weighted by map-quality tier so its mass is shared across the loci it maps to, reproducing pre-tars behaviour. Under either aligner, multi-mapped reads are excluded from novel splice junction and chimeric analysis. The optional `MULTI_MAP_LOCI` write type emits a tsv of each multi-mapped read's primary and XA alternate loci per gene collection for auditing (bwa-tars only).
 
 ## Configuration
 The functions of Isofox are controlled by the 'functions' argument:
@@ -74,7 +64,7 @@ excluded_regions | Drop reads in regions of high multi-mappability
 excluded_gene_id_file | Exclude genes in file, format EnsemblGeneId,GeneName
 enriched_gene_ids | List of EnsemblGeneIds separated by ';', see Enriched Genes information below
 drop_dups | Default is false. By default duplicate fragments will be counted towards transcript expression.
-write_type | Debug only: EXON, SPLICE_JUNC, FRAG_LENGTH, FRAG_LENGTH_BY_GENE, READ, CHIMERIC_READ, CHIMERIC_POSITION_DATA, SPLICE_SITE, TRANS_COMBO;GC_RATIO, separated by ';'
+write_type | Debug only, separated by ';': EXON, SPLICE_JUNC, FRAG_LENGTH, FRAG_LENGTH_BY_GENE, READ, CHIMERIC_READ, CHIMERIC_POSITION_DATA, SPLICE_SITE, TRANS_COMBO, GC_RATIO, MULTI_MAP_LOCI
 
 ### Reference Files
 
@@ -98,7 +88,8 @@ frag_length_min_count | Minimum number of fragments to observe for length distri
 exp_rate_frag_lengths | Discrete buckets for fragment lengths, either with frequency specified or left as zero if to be calculated (ie with -apply_calc_frag_lengths). eg '50-0;75-0;100-0;125-0;150-0;200-0;250-0;300-0;400-0;550-0' 
 read_length | Expected RNA read length (eg 76 or 151), will be computed if not provided
 long_frag_limit | Default 550 bases, fragments longer than this without a splice junction are not considered to support a gene for the purposes of expression
-single_map_qual | Default 255, discard reads with map quality below this unless using the config 'apply_map_qual_adjust'
+single_map_qual | Default 60, or 255 when 'star_aligner' is set. Discard reads with map quality below this.
+star_aligner | Flag, off by default. Off assumes the bwa-tars + REDUX pipeline; set it to apply STAR-aligner settings (raises the default single_map_qual to 255).
 enriched_gene_ids | By default includes: ENSG00000265150;ENSG00000258486;ENSG00000202198;ENSG00000266037;ENSG00000263740;ENSG00000265735
 
 ### Optional output files:
@@ -236,7 +227,7 @@ Each fragment is assigned to a 'category' based on the set of transcripts that i
 
 Any fragment which does not contain a splice junction, is wholly contained within the bounds of a gene, and with fragment size <= maximum insert size distribution is also allowed to map to an ‘UNSPLICED’ transcript of that gene.
 
-Note that reads which are partially exonic, but marginally overhang an exon boundary or are soft clipped at or beyond an exon boundary have special treatment. This is particularly relevant for reads that have an overhang of 1 or 2 bases which will not be mapped by STAR with default parameters. If the overhanging section can be uniquely mapped either to the reference or to the other side of only a single known spliced junction, then the fragment is deemed to be supporting that splice junction or in the case of supporting just the reference is deemed to be supporting the UNSPLICED transcript.  If multiple mappings are possible or the fragment length < unclipped read length (indicating likely adapter sequence) it is truncated at the exon boundary.  If no mapping is possible then the fragment is treated as not supporting any known transcript.
+Note that reads which are partially exonic, but marginally overhang an exon boundary or are soft clipped at or beyond an exon boundary have special treatment. This is particularly relevant for reads that have an overhang of 1 or 2 bases which may not be mapped across the exon boundary by the aligner. If the overhanging section can be uniquely mapped either to the reference or to the other side of only a single known spliced junction, then the fragment is deemed to be supporting that splice junction or in the case of supporting just the reference is deemed to be supporting the UNSPLICED transcript.  If multiple mappings are possible or the fragment length < unclipped read length (indicating likely adapter sequence) it is truncated at the exon boundary.  If no mapping is possible then the fragment is treated as not supporting any known transcript.
 
 ### 5. Fit abundance estimate per transcript