From ad30f46a1fab28f7ab85a098f075a6fd04cbc520 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Fri, 19 Jun 2026 12:37:02 +1000 Subject: [PATCH 1/8] Esvee: removed private method only --- .../hmftools/esvee/assembly/types/JunctionAssembly.java | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java b/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java index ef5f89b38e..a2fe9158f6 100644 --- a/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java +++ b/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java @@ -851,18 +851,13 @@ public SagaMatchBySequence sagaMatch() return mSagaMatch; } - private void setSagaMatch(final SagaMatchBySequence match) - { - mSagaMatch = match; - } - public boolean matchToSaga(final SagaSequenceMatcher sagaMatcher) { int junctionOffset = mJunction.isForward() ? refBaseLength() : baseLength() - refBaseLength(); List junctionInfos = List.of(new SagaJunctionInfo(junctionOffset)); boolean lowerJunctionOverlap = hasLineSequence(); SagaMatchBySequence match = sagaMatcher.matchBySequence(bases(), junctionInfos, lowerJunctionOverlap, true); - setSagaMatch(match); + mSagaMatch = match; return match != null; } From 238af7a24d8a3714761c14543db42c14998a0006 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Mon, 22 Jun 2026 15:31:42 +1000 Subject: [PATCH 2/8] Esvee: removed redundant code --- .../esvee/assembly/types/JunctionAssembly.java | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java b/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java index a2fe9158f6..c442982bb9 100644 --- a/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java +++ b/esvee/src/main/java/com/hartwig/hmftools/esvee/assembly/types/JunctionAssembly.java @@ -11,14 +11,14 @@ import static com.hartwig.hmftools.esvee.assembly.AssemblyUtils.readQualFromJunction; import static com.hartwig.hmftools.esvee.assembly.IndelBuilder.convertedIndelCrossesJunction; import static com.hartwig.hmftools.esvee.assembly.IndelBuilder.findInsertedBases; -import static com.hartwig.hmftools.esvee.common.CommonUtils.aboveMinQual; -import static com.hartwig.hmftools.esvee.common.CommonUtils.belowMinQual; -import static com.hartwig.hmftools.esvee.common.SvConstants.MIN_INDEL_LENGTH; -import static com.hartwig.hmftools.esvee.common.SvConstants.MIN_VARIANT_LENGTH; import static com.hartwig.hmftools.esvee.assembly.types.AssemblyOutcome.UNSET; import static com.hartwig.hmftools.esvee.assembly.types.RepeatInfo.findRepeats; import static com.hartwig.hmftools.esvee.assembly.types.SupportType.INDEL; import static com.hartwig.hmftools.esvee.assembly.types.SupportType.JUNCTION; +import static com.hartwig.hmftools.esvee.common.CommonUtils.aboveMinQual; +import static com.hartwig.hmftools.esvee.common.CommonUtils.belowMinQual; +import static com.hartwig.hmftools.esvee.common.SvConstants.MIN_INDEL_LENGTH; +import static com.hartwig.hmftools.esvee.common.SvConstants.MIN_VARIANT_LENGTH; import java.util.List; import java.util.Set; @@ -166,12 +166,12 @@ public JunctionAssembly( public int mergedAssemblyCount() { return mMergedAssemblies; } public void addMergedAssembly() { ++mMergedAssemblies; } - public int junctionIndex() { return mJunctionIndex; }; - public void setJunctionIndex(int index) { mJunctionIndex = index; }; + public int junctionIndex() { return mJunctionIndex; } + public void setJunctionIndex(int index) { mJunctionIndex = index; } // eg 21 bases, junction index at 10 (so 0-9 = 10 before, 11-20 = 10 after), note: doesn't count the junction base - public int lowerDistanceFromJunction() { return mJunctionIndex; }; - public int upperDistanceFromJunction() { return mBases.length - mJunctionIndex - 1; }; + public int lowerDistanceFromJunction() { return mJunctionIndex; } + public int upperDistanceFromJunction() { return mBases.length - mJunctionIndex - 1; } public int refBaseLength() { return (mJunction.isForward() ? lowerDistanceFromJunction() : upperDistanceFromJunction()) + 1; } From 04babcbb976691c556a68ce48bb2be43638af4f6 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Wed, 30 Sep 2026 06:41:52 +1000 Subject: [PATCH 3/8] Isofox: - removed BamReadCounter - updated retained intron finder to write TSV - updated read-me for Tars+Redux changes --- isofox/README.md | 32 +-- .../com/hartwig/hmftools/isofox/Isofox.java | 27 --- .../hmftools/isofox/IsofoxFunction.java | 1 - .../com/hartwig/hmftools/isofox/TaskType.java | 3 +- .../isofox/common/BamReadCounter.java | 214 ------------------ .../isofox/novel/RetainedIntronFinder.java | 33 ++- .../isofox/results/ResultsWriter.java | 7 +- 7 files changed, 43 insertions(+), 274 deletions(-) delete mode 100644 isofox/src/main/java/com/hartwig/hmftools/isofox/common/BamReadCounter.java diff --git a/isofox/README.md b/isofox/README.md index 1afdaea08a..584c454224 100644 --- a/isofox/README.md +++ b/isofox/README.md @@ -14,17 +14,19 @@ For transcript abundance, Isofox uses a similar methodology to several previous 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 -We recommend to mark duplicates in your pipeline. They are included in gene and transcript expression data (to avoid bias against highly expressed genes) but excluded from novel splice junction analysis. - -We find that 6 genes in particular (RN7SL2, RN7SL1, RN7SL3, RN7SL4P, RN7SL5P & RN7SK) are highly expressed across our cohort and at variable rates - in extreme samples these can account for >75% of all transcripts. Isofox excludes these genes from our GC bias calculations and to determine a normalisation factor for "adjusted TPM" so that they don't dominate expression differences. For any given sample, AdjustedTPM = rawTPM x constant with the constant determined by the normalisation (which excludes the 6 genes and also limits all other genes to 1% contribution). The adjusted TPMs no longer sum to 1M transcripts, but should be more comparable across samples. We suggest to use the adjusted TPM for expression analysis. - -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 duplicates +Duplicates are marked by Redux and are counted towards transcript expression. ### A note on alignment and multi-mapping 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. -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). +Isofox from v2.1 onwards only supports both alignments from bwa-mem2 + 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. +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. +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: @@ -115,7 +117,7 @@ write_read_data | Write data on each BAM read, only recommended with restricted write_exon_data | Write data on transcript exon covered by a supporting fragment, only recommended with restricted genes file ### Memory Usage and Threading -ISOFOX takes ~10 mins to process a 7GB BAM with 120M reads / 60M fragments using 10 cores, with maximum memory usage of 10GB, and ~30 mins to process a 35GB BAM with 440M reads / 200M fragments using 10 cores, with maximum memory usage of 25GB. +ISOFOX takes ~5 mins to process a 20GB BAM with 120M reads using 32 cores. Recommend 64GB memory. ### Example Usage Running all functions: @@ -346,7 +348,7 @@ Each chimeric junction, novel splice junction and retained intron for each sampl ### Summary -Generated file: sample_id.isf.summary.csv +Generated file: sample_id.isf.summary.tsv Field | Description ---|--- @@ -360,12 +362,12 @@ ReadLength | Raw read length of fragments FragLength5th | 5th percentile of genic intronic fragment lengths (from 1M fragments sampled with a max of 1000 per gene) FragLength50th | 50th percentile of genic intronic fragment lengths (from 1M fragments sampled with a max of 1000 per gene) FragLength95th | 95th percentile of genic intronic fragment lengths (from 1M fragments sampled with a max of 1000 per gene) -EnrichedGenePercent | % of fragments supporting one of the following 6 genes: (RN7SL2, RN7SL1,RN7SL3,RN7SL4P,RN7SL5P & RN7SK) MedianGCRatio | Median GC ratio excluding the 6 highly enriched genes +ForwardStrandPercent | Percent of fragments in the forward strand direction, ie F1R2 and not F2R1 ### Gene Level Data -Generated file: sample_id.isf.gene_data.csv +Generated file: sample_id.isf.gene_data.tsv Field | Description ---|--- @@ -382,7 +384,7 @@ TPM | TPM for gene excluding unspliced fragments ### Transcript Level Data -Generated file: sample_id.isf.trans_data.csv +Generated file: sample_id.isf.trans_data.tsv Field | Description ---|--- @@ -406,7 +408,7 @@ UniqueNonSJFragments | Count of fragments uniquely supporting transcript but wit ### Fragment length distribution -Generated file: sample_id.isf.frag_length.csv +Generated file: sample_id.isf.frag_length.tsv Field | Description ---|--- @@ -415,7 +417,7 @@ Count | Count of fragments with specified fragment length ### Alternate Splice Junctions -Generated file: sample_id.isf.alt_splice_junc.csv +Generated file: sample_id.isf.alt_splice_junc.tsv Field | Description ---|--- @@ -441,7 +443,7 @@ OverlappingGenes | List of all genes which overlap the novel splice junction ### Retained Introns -Generated file: sample_id.isf.retained_intron.csv +Generated file: sample_id.isf.retained_intron.tsv Field | Description ---|--- diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/Isofox.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/Isofox.java index aaa448dbb6..ff3c17d13b 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/Isofox.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/Isofox.java @@ -14,7 +14,6 @@ import static com.hartwig.hmftools.isofox.IsofoxConstants.PRIORITISED_CHROMOSOMES; import static com.hartwig.hmftools.isofox.IsofoxFunction.FUSIONS; import static com.hartwig.hmftools.isofox.IsofoxFunction.NEO_EPITOPES; -import static com.hartwig.hmftools.isofox.IsofoxFunction.READ_COUNTS; import static com.hartwig.hmftools.isofox.TaskType.APPLY_GC_ADJUSTMENT; import static com.hartwig.hmftools.isofox.TaskType.TRANSCRIPT_COUNTS; import static com.hartwig.hmftools.isofox.WriteType.FRAG_LENGTH; @@ -47,7 +46,6 @@ import com.hartwig.hmftools.isofox.adjusts.FragmentSizeCalcs; import com.hartwig.hmftools.isofox.adjusts.GcRatioCounts; import com.hartwig.hmftools.isofox.adjusts.GcTranscriptCalculator; -import com.hartwig.hmftools.isofox.common.BamReadCounter; import com.hartwig.hmftools.isofox.common.FragmentTypeCounts; import com.hartwig.hmftools.isofox.common.PerformanceTracking; import com.hartwig.hmftools.isofox.expression.ExpectedCountsCache; @@ -130,13 +128,6 @@ public boolean runAnalysis() return true; } - if(mConfig.runFunction(READ_COUNTS)) - { - boolean status = countBamReads(chrGeneMap); - mResultsWriter.close(); - return status; - } - // BAM processing for the key routines - novel junctions, fusions and gene expression if(!allocateBamFragments(chrGeneMap)) return false; @@ -438,24 +429,6 @@ private void calcFragmentLengths(final Map> chrGeneMap) } } - private boolean countBamReads(final Map> chrGeneMap) - { - ISF_LOGGER.info("basic BAM read counts"); - - List taskList = Lists.newArrayList(); - List> callableList = Lists.newArrayList(); - - for(Map.Entry> entry : chrGeneMap.entrySet()) - { - BamReadCounter bamReaderTask = new BamReadCounter(mConfig, mResultsWriter); - bamReaderTask.initialise(entry.getKey(), entry.getValue()); - taskList.add(bamReaderTask); - callableList.add(bamReaderTask); - } - - return TaskExecutor.executeTasks(callableList, mConfig.Threads); - } - public static void main(@NotNull final String[] args) { ConfigBuilder configBuilder = new ConfigBuilder(APP_NAME); diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/IsofoxFunction.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/IsofoxFunction.java index dfec8d5aa4..1e6242fdb3 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/IsofoxFunction.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/IsofoxFunction.java @@ -9,7 +9,6 @@ public enum IsofoxFunction RETAINED_INTRONS, FUSIONS, STATISTICS, - READ_COUNTS, NEO_EPITOPES; public static final List DEFAULT_FUNCTIONS = List.of(TRANSCRIPT_COUNTS, ALT_SPLICE_JUNCTIONS, FUSIONS); diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/TaskType.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/TaskType.java index 2c47f31d70..883f0d0fa6 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/TaskType.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/TaskType.java @@ -6,6 +6,5 @@ public enum TaskType TRANSCRIPT_COUNTS, GENERATE_GC_COUNTS, GENERATE_EXPECTED_COUNTS, - APPLY_GC_ADJUSTMENT, - BAM_READ_COUNTER; + APPLY_GC_ADJUSTMENT; } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/BamReadCounter.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/BamReadCounter.java deleted file mode 100644 index 2ea822deaf..0000000000 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/BamReadCounter.java +++ /dev/null @@ -1,214 +0,0 @@ -package com.hartwig.hmftools.isofox.common; - -import static java.lang.Math.max; -import static java.lang.Math.min; - -import static com.hartwig.hmftools.common.bam.SamRecordUtils.SUPPLEMENTARY_ATTRIBUTE; -import static com.hartwig.hmftools.common.region.BaseRegion.positionWithin; -import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.createBufferedWriter; -import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_END; -import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_PAIR; -import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_START; -import static com.hartwig.hmftools.isofox.ChromosomeTaskExecutor.findNextOverlappingGenes; -import static com.hartwig.hmftools.isofox.IsofoxConfig.ISF_LOGGER; -import static com.hartwig.hmftools.isofox.WriteType.READ; -import static com.hartwig.hmftools.isofox.common.FragmentType.CHIMERIC; -import static com.hartwig.hmftools.isofox.common.FragmentType.DUPLICATE; -import static com.hartwig.hmftools.isofox.common.FragmentType.TOTAL; - -import java.io.BufferedWriter; -import java.io.File; -import java.io.IOException; -import java.util.List; -import java.util.concurrent.Callable; - -import com.google.common.collect.Lists; -import com.hartwig.hmftools.common.bam.BamSlicer; -import com.hartwig.hmftools.common.bam.SupplementaryReadData; -import com.hartwig.hmftools.common.gene.GeneData; -import com.hartwig.hmftools.common.region.ChrBaseRegion; -import com.hartwig.hmftools.isofox.IsofoxConfig; -import com.hartwig.hmftools.isofox.results.ResultsWriter; - -import htsjdk.samtools.SAMFlag; -import htsjdk.samtools.SAMRecord; -import htsjdk.samtools.SamReader; -import htsjdk.samtools.SamReaderFactory; - -// simple BAM read counter, used for experimental purposes at the moment -public class BamReadCounter implements Callable -{ - private final IsofoxConfig mConfig; - - private final SamReader mSamReader; - private final BamSlicer mBamSlicer; - private final ResultsWriter mResultsWriter; - - private final int[] mCurrentGenesRange; - private int mTotalReadCount; - private int mCurrentGeneReadCount; - private final FragmentTypeCounts mFragmentTypeCounts; - private String mChromosome; - private final List mGeneDataList; - private String mCurrentGenes; - private final int[] mMaqQualFrequencies; - - public BamReadCounter(final IsofoxConfig config, final ResultsWriter resultsWriter) - { - mConfig = config; - mSamReader = mConfig.BamFile != null ? - SamReaderFactory.makeDefault().referenceSequence(mConfig.RefGenomeFile).open(new File(mConfig.BamFile)) : null; - - mBamSlicer = new BamSlicer(0, true, true, true); - mResultsWriter = resultsWriter; - - mGeneDataList = Lists.newArrayList(); - mChromosome = ""; - - mCurrentGenesRange = new int[SE_PAIR]; - mTotalReadCount = 0; - mCurrentGeneReadCount = 0; - mCurrentGenes = ""; - mFragmentTypeCounts = new FragmentTypeCounts(); - mMaqQualFrequencies = new int[4]; - } - - public void initialise(final String chromosome, final List geneDataList) - { - mChromosome = chromosome; - mGeneDataList.clear(); - mGeneDataList.addAll(geneDataList); - } - - @Override - public Void call() - { - processBam(); - return null; - } - - private void processBam() - { - if(mGeneDataList.isEmpty()) - return; - - // walk through each chromosome, taking groups of overlapping genes together - ISF_LOGGER.info("processing reads for chromosome({}) geneCount({})", mChromosome, mGeneDataList.size()); - - final List overlappingGenes = Lists.newArrayList(); - int currentGeneIndex = 0; - int nextLogCount = 100; - - while(currentGeneIndex < mGeneDataList.size()) - { - currentGeneIndex = findNextOverlappingGenes(mGeneDataList, currentGeneIndex, overlappingGenes); - - // if(overlappingGenes.stream().anyMatch(x -> mConfig.Filters.EnrichedGeneIds.contains(x.GeneId))) - // continue; - - mCurrentGenesRange[SE_START] = 0; - mCurrentGenesRange[SE_END] = 0; - - for(int i = 0; i < overlappingGenes.size(); ++i) - { - GeneData geneData = overlappingGenes.get(i); - - mCurrentGenesRange[SE_START] = i == 0 ? geneData.GeneStart : min(geneData.GeneStart, mCurrentGenesRange[SE_START]); - mCurrentGenesRange[SE_END] = i == 0 ? geneData.GeneEnd : max(geneData.GeneEnd, mCurrentGenesRange[SE_END]); - } - - if(currentGeneIndex >= nextLogCount) - { - nextLogCount += 100; - ISF_LOGGER.debug("chromosome({}) processed {} genes, totalReads({})", - mChromosome, currentGeneIndex, mTotalReadCount); - } - - mCurrentGenes = overlappingGenes.get(0).GeneId; - mCurrentGeneReadCount = 0; - - final List regions = Lists.newArrayList(new ChrBaseRegion(mChromosome, mCurrentGenesRange)); - - mBamSlicer.slice(mSamReader, regions, this::processBamRead); - } - - ISF_LOGGER.info("chromosome({}) processing complete: total({}) duplicates({}) chimeric({}) mapQuals(0={} 1={} 2={} 3={})", - mChromosome, mTotalReadCount, mFragmentTypeCounts.typeCount(DUPLICATE), mFragmentTypeCounts.typeCount(CHIMERIC), - mMaqQualFrequencies[0], mMaqQualFrequencies[1], mMaqQualFrequencies[2], mMaqQualFrequencies[3]); - } - - private void processBamRead(final SAMRecord record) - { - if(!positionWithin(record.getStart(), mCurrentGenesRange[SE_START], mCurrentGenesRange[SE_END])) - return; - - if(mConfig.Filters.skipRead(record.getContig(), record.getAlignmentStart())) - return; - - ++mTotalReadCount; - ++mCurrentGeneReadCount; - mFragmentTypeCounts.addCount(TOTAL); - - if(record.getDuplicateReadFlag()) - mFragmentTypeCounts.addCount(DUPLICATE); - - if((record.getFlags() & SAMFlag.SUPPLEMENTARY_ALIGNMENT.intValue()) != 0) - mFragmentTypeCounts.addCount(CHIMERIC); - - if(record.getMappingQuality() <= 3) - { - mMaqQualFrequencies[record.getMappingQuality()]++; - } - - if(mConfig.GeneReadLimit > 0 && mCurrentGeneReadCount > mConfig.GeneReadLimit) - { - mBamSlicer.haltProcessing(); - ISF_LOGGER.info("chromosome({}) gene({}) halting processing after {} reads", mChromosome, mCurrentGenes, mCurrentGeneReadCount); - } - - if(mConfig.writeType(READ)) - writeReadData(mResultsWriter.getReadDataWriter(), record, mCurrentGenes); - } - - public static BufferedWriter createReadDataWriter(final IsofoxConfig config) - { - try - { - final String outputFileName = config.formOutputFile("read_data.csv"); - - BufferedWriter writer = createBufferedWriter(outputFileName, false); - writer.write("GeneId,ReadId,Chromosome,PosStart,PosEnd,Cigar,Flags,InsertSize"); - writer.write(",MateChr,MatePosStart,FirstInPair,ReadReversed,Duplicate,Supplementary,SuppData"); - writer.newLine(); - return writer; - } - catch(IOException e) - { - ISF_LOGGER.error("failed to create read data writer: {}", e.toString()); - return null; - } - } - - private static synchronized void writeReadData(final BufferedWriter writer, final SAMRecord record, final String geneId) - { - try - { - SupplementaryReadData suppData = SupplementaryReadData.extractAlignment(record.getStringAttribute(SUPPLEMENTARY_ATTRIBUTE)); - - writer.write(String.format("%s,%s,%s,%d,%d,%s,%d,%d,%s,%d", - geneId, record.getReadName(), record.getContig(), record.getAlignmentStart(), record.getAlignmentEnd(), - record.getCigarString(), record.getFlags(), record.getInferredInsertSize(), - record.getMateReferenceName(), record.getMateAlignmentStart())); - - writer.write(String.format(",%s,%s,%s,%s,%s,%s", - record.getFirstOfPairFlag(), record.getReadNegativeStrandFlag(), record.getDuplicateReadFlag(), - record.getSupplementaryAlignmentFlag(), suppData != null ? suppData.asDelimStr() : "N/A")); - - writer.newLine(); - } - catch(IOException e) - { - ISF_LOGGER.error("failed to write read data file: {}", e.toString()); - } - } -} diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java index 31f4fe47ea..c3f8093fdc 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java @@ -2,6 +2,11 @@ import static java.lang.Math.max; +import static com.hartwig.hmftools.common.utils.file.CommonFields.FLD_CHROMOSOME; +import static com.hartwig.hmftools.common.utils.file.CommonFields.FLD_GENE_ID; +import static com.hartwig.hmftools.common.utils.file.CommonFields.FLD_GENE_NAME; +import static com.hartwig.hmftools.common.utils.file.CommonFields.FLD_POSITION; +import static com.hartwig.hmftools.common.utils.file.FileDelimiters.TSV_DELIM; import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.createBufferedWriter; import static com.hartwig.hmftools.common.region.BaseRegion.positionWithin; import static com.hartwig.hmftools.isofox.IsofoxConfig.ISF_LOGGER; @@ -13,6 +18,7 @@ import java.io.IOException; import java.util.List; import java.util.Map; +import java.util.StringJoiner; import java.util.stream.Collectors; import com.google.common.collect.Lists; @@ -221,11 +227,14 @@ public static BufferedWriter createWriter(final IsofoxConfig config) { try { - String outputFileName = config.formOutputFile("retained_intron.csv"); + String outputFileName = config.formOutputFile("retained_intron.tsv"); BufferedWriter writer = createBufferedWriter(outputFileName, false); - writer.write("GeneId,GeneName,Chromosome,Strand,Position"); - writer.write(",Type,FragCount,SplicedFragCount,TotalDepth,TranscriptInfo"); + + StringJoiner sj = new StringJoiner(TSV_DELIM); + sj.add(FLD_GENE_ID).add(FLD_GENE_NAME).add(FLD_CHROMOSOME).add("Strand").add(FLD_POSITION); + sj.add("Type").add("FragCount").add("SplicedFragCount").add("TotalDepth").add("TranscriptInfo"); + writer.write(sj.toString()); writer.newLine(); return writer; } @@ -260,16 +269,22 @@ private synchronized static void writeRetainedIntrons( if(!gene.getTranscripts().stream().anyMatch(x -> retIntron.regions().stream().anyMatch(y -> y.hasTransId(x.TransId)))) continue; - writer.write(String.format("%s,%s,%s,%d", - gene.Gene.GeneId, gene.Gene.GeneName, - gene.Gene.Chromosome, gene.Gene.Strand)); + StringJoiner sj = new StringJoiner(TSV_DELIM); + sj.add(gene.Gene.GeneId); + sj.add(gene.Gene.GeneName); + sj.add(gene.Gene.Chromosome); + sj.add(String.valueOf(gene.Gene.Strand)); + sj.add(String.valueOf(retIntron.position())); + sj.add(String.valueOf(retIntron.type(gene.Gene.forwardStrand()))); + sj.add(String.valueOf(retIntron.getFragmentCount())); + sj.add(String.valueOf(retIntron.getSplicedFragmentCount())); int readDepth = max(retIntron.getDepth(), retIntron.getFragmentCount()); + sj.add(String.valueOf(readDepth)); - writer.write(String.format(",%d,%s,%d,%d,%d,%s", - retIntron.position(), retIntron.type(gene.Gene.forwardStrand()), retIntron.getFragmentCount(), - retIntron.getSplicedFragmentCount(), readDepth, retIntron.transcriptInfo())); + sj.add(retIntron.transcriptInfo()); + writer.write(sj.toString()); writer.newLine(); } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/results/ResultsWriter.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/results/ResultsWriter.java index 5af1086e47..b2265e40d4 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/results/ResultsWriter.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/results/ResultsWriter.java @@ -29,7 +29,6 @@ import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.createBufferedWriter; import static com.hartwig.hmftools.isofox.IsofoxConfig.ISF_LOGGER; import static com.hartwig.hmftools.isofox.IsofoxFunction.ALT_SPLICE_JUNCTIONS; -import static com.hartwig.hmftools.isofox.IsofoxFunction.READ_COUNTS; import static com.hartwig.hmftools.isofox.IsofoxFunction.RETAINED_INTRONS; import static com.hartwig.hmftools.isofox.IsofoxFunction.TRANSCRIPT_COUNTS; import static com.hartwig.hmftools.isofox.common.FragmentType.ALT; @@ -72,7 +71,6 @@ import com.hartwig.hmftools.common.rna.RnaStatistics; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.adjusts.FragmentSizeCalcs; -import com.hartwig.hmftools.isofox.common.BamReadCounter; import com.hartwig.hmftools.isofox.common.FragmentTypeCounts; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.GeneReadData; @@ -172,10 +170,7 @@ private void initialiseExternalWriters() if(mConfig.writeType(READ)) { - if(mConfig.runFunction(READ_COUNTS)) - mReadDataWriter = BamReadCounter.createReadDataWriter(mConfig); - else - mReadDataWriter = createReadDataWriter(mConfig); + mReadDataWriter = createReadDataWriter(mConfig); } if(mConfig.writeType(SPLICE_SITE)) From 7f013420ae9f66a16ebf7a7fc27018a3ca1a3328 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Wed, 30 Sep 2026 08:39:32 +1000 Subject: [PATCH 4/8] Isofox: fixed no setting genic data for chimeric primary reads --- .../hmftools/isofox/FragmentAllocator.java | 40 +++++++++++++------ .../isofox/common/FragmentTracker.java | 2 + 2 files changed, 30 insertions(+), 12 deletions(-) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java index 412a955087..87be263a37 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java @@ -220,10 +220,7 @@ public void processBam(final GeneCollection geneCollection, final ChrBaseRegion mBamSlicer.slice(mSamReader, geneRegion, this::processSamRecord); if(mChimericReads.enabled()) - { - mChimericReads.postProcessChimericReads(mBaseDepth, mFragmentReads); - processChimericNovelJunctions(); - } + processIncompleteReads(); ISF_LOGGER.trace("genes({}) bamReadCount({}) depth(bases={} perc={} max={})", mCurrentGenes.geneNames(), mGeneReadCount, mBaseDepth.basesWithDepth(), @@ -329,7 +326,7 @@ private void markGeneDataRegions(final Read read) private void handleSupplementaryRead(final Read read) { - if(read.isDuplicate() || read.isMateUnmapped()) + if(read.isMateUnmapped()) return; mBaseDepth.processRead(read.getMappedRegionCoords()); @@ -337,13 +334,7 @@ private void handleSupplementaryRead(final Read read) markGeneDataRegions(read); mChimericReads.addSupplementaryRead(read); - if(mReadDataWriter != null && mConfig.writeType(WriteType.READ)) - { - List overlapGenes = mCurrentGenes.findGenesCoveringRange( - read.alignmentStart(), read.alignmentEnd(), true); - - writeReadData(mReadDataWriter, overlapGenes, read, CHIMERIC, 0); - } + writeChimericReadData(read); } private void processFragmentReads(final Read read1, final Read read2) @@ -852,6 +843,21 @@ private void processIntronicReads( } } + private void processIncompleteReads() + { + // now slicing has completed, process unpaired primaries (supps have been handled already) + for(Object readObject : mFragmentReads.readMap().values()) + { + Read read = (Read)readObject; + markGeneDataRegions(read); + mBaseDepth.processRead(read.getMappedRegionCoords()); + writeChimericReadData(read); + } + + mChimericReads.postProcessChimericReads(mBaseDepth, mFragmentReads); + processChimericNovelJunctions(); + } + private void processChimericNovelJunctions() { // examine chimeric reads to see if they can instead be handled as novel alternate splicing @@ -932,6 +938,16 @@ public void registerKnownFusionPairs(final EnsemblDataCache geneTransCache) mChimericReads.registerKnownFusionPairs(geneTransCache); } + private void writeChimericReadData(final Read read) + { + if(mReadDataWriter == null) + return; + + List overlapGenes = mCurrentGenes.findGenesCoveringRange( + read.alignmentStart(), read.alignmentEnd(), true); + + writeReadData(mReadDataWriter, overlapGenes, read, CHIMERIC, 0); + } @VisibleForTesting public void processReadRecords(final GeneCollection geneCollection, final List reads) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/FragmentTracker.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/FragmentTracker.java index cbe7c54983..bf8c9fa326 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/FragmentTracker.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/FragmentTracker.java @@ -15,6 +15,8 @@ public FragmentTracker() mReadMap = Maps.newHashMap(); } + public Map readMap() { return mReadMap; } + public List getValues() { return mReadMap.values().stream().collect(Collectors.toList()); } public int readsCount() { return mReadMap.size(); } From f53d6531d7a79262b98c5ba622ca69d743f8d99b Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Wed, 30 Sep 2026 10:54:10 +1000 Subject: [PATCH 5/8] Isofox: added unit test for previous issue --- .../hmftools/isofox/FragmentAllocator.java | 6 ++ .../isofox/fusion/ChimericReadTest.java | 2 +- .../isofox/fusion/FusionFiltersTest.java | 98 ++++++++++++++++++- 3 files changed, 104 insertions(+), 2 deletions(-) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java index 87be263a37..f480d1081d 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java @@ -968,6 +968,12 @@ public void processReadRecords(final GeneCollection geneCollection, final List processRead(x)); } + @VisibleForTesting + public void postSliceProcessReads() + { + processIncompleteReads(); + } + @VisibleForTesting public final FragmentTracker getFragmentTracker() { return mFragmentReads; } diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java index 03cc496d77..fbe115f234 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java @@ -475,7 +475,7 @@ public void testMultiGeneChimericRead() assertFalse(ChimericUtils.setHasMultipleKnownSpliceGenes(Lists.newArrayList(read), knownPairGeneIds)); - // if the genes are known then treat this as chimeroc + // if the genes are known then treat this as chimeric knownPairGeneIds.add(new String[] {GENE_ID_1, GENE_ID_2}); assertTrue(ChimericUtils.setHasMultipleKnownSpliceGenes(Lists.newArrayList(read), knownPairGeneIds)); } diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionFiltersTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionFiltersTest.java index dd7234c6d7..ec354dd3ff 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionFiltersTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionFiltersTest.java @@ -4,6 +4,7 @@ import static com.hartwig.hmftools.isofox.TestUtils.ALT_SJ_COHORT_CACHE; import static com.hartwig.hmftools.isofox.TestUtils.CHR_1; import static com.hartwig.hmftools.isofox.TestUtils.CHR_2; +import static com.hartwig.hmftools.isofox.TestUtils.GENE_ID_1; import static com.hartwig.hmftools.isofox.TestUtils.GENE_ID_3; import static com.hartwig.hmftools.isofox.TestUtils.GENE_ID_5; import static com.hartwig.hmftools.isofox.TestUtils.addTestGenes; @@ -15,14 +16,18 @@ import static com.hartwig.hmftools.isofox.TestUtils.createSupplementaryReadPair; import static com.hartwig.hmftools.isofox.TestUtils.populateRefGenome; import static com.hartwig.hmftools.isofox.TestUtils.setReadFirstSecondInPair; +import static com.hartwig.hmftools.isofox.fusion.FusionJunctionType.KNOWN; import static com.hartwig.hmftools.isofox.fusion.FusionTestUtils.createGeneDataCache; import static junit.framework.TestCase.assertEquals; import static junit.framework.TestCase.assertTrue; +import java.util.Collection; +import java.util.Collections; import java.util.List; import java.util.Map; import java.util.Set; +import java.util.stream.Collectors; import com.google.common.collect.Lists; import com.google.common.collect.Maps; @@ -206,7 +211,7 @@ public void testFusionHardFilters() readId1, gc5, gc3, 10200, 10219, 20281, 20300, createCigar(20, 20, 0), createCigar(0, 20, 20), true); - // TOD): should be setting both primary and supp to have the same strandedness + // TODO: should be setting both primary and supp to have the same strandedness readPair1[0].setStrand(true, false); bamReader1.processReadRecords(gc5, Lists.newArrayList(read1, readPair1[0])); @@ -304,4 +309,95 @@ public void testFusionHardFilters() assertEquals(1, finderChr2.getFusionCandidates().values().stream().mapToInt(x -> x.size()).sum()); } + @Test + public void testFusionSplitPrimaries() + { + // a fusion formed from the primaries being in different gene collections + EnsemblDataCache geneTransCache = createGeneDataCache(); + + addTestGenes(geneTransCache); + addTestTranscripts(geneTransCache); + + IsofoxConfig config = createIsofoxConfig(); + config.Functions.clear(); + config.Functions.add(FUSIONS); + config.Fusions.MinHardFilterFrags = 1; + + populateRefGenome(config.RefGenome); + + int gcId = 0; + + // the mate read + GeneCollection gc1 = createGeneCollection(geneTransCache, gcId++, Lists.newArrayList(geneTransCache.getGeneDataById(GENE_ID_1))); + + // the primary for the fusion + GeneCollection gc3 = createGeneCollection(geneTransCache, gcId++, Lists.newArrayList(geneTransCache.getGeneDataById(GENE_ID_3))); + + // the supp for the fusion + GeneCollection gc5 = createGeneCollection(geneTransCache, gcId, Lists.newArrayList(geneTransCache.getGeneDataById(GENE_ID_5))); + + FragmentAllocator bamReader = new FragmentAllocator(config, geneTransCache, ALT_SJ_COHORT_CACHE, new ResultsWriter(config)); + + FusionTaskManager fusionTaskManager = new FusionTaskManager(config, geneTransCache); + + FusionFinder finderChr1 = fusionTaskManager.createFusionFinder(gc3.chromosome()); + + int readId = 1; + + // the mate read + Read mateRead = createMappedRead(readId, gc1, 1200, 1240, createCigar(0, 41, 0)); + setReadFirstSecondInPair(mateRead, false); + + Read[] readPair = createSupplementaryReadPair( + readId, gc3, gc5, 20261, 20300, 10600, 10619, + createCigar(0, 40, 20), createCigar(40, 20, 00), true); + + bamReader.processReadRecords(gc1, Lists.newArrayList(mateRead)); + bamReader.postSliceProcessReads(); + + finderChr1.processNewChimericReadGroups(gc1, bamReader.getBaseDepth(), bamReader.getChimericReadTracker().fusionReadGroupMap()); + + bamReader.clearCache(); + + bamReader.processReadRecords(gc3, Lists.newArrayList(readPair[0])); + bamReader.postSliceProcessReads(); + + finderChr1.processNewChimericReadGroups(gc3, bamReader.getBaseDepth(), bamReader.getChimericReadTracker().fusionReadGroupMap()); + + Map> chrIncompleteReadsGroups = finderChr1.extractIncompleteReadGroups( + gc3.chromosome(), bamReader.getChimericReadTracker().getHardFilteredReadIds()); + + // note the chromosome is the one the imcomplete groups link to + List interChromosomalGroups = fusionTaskManager.addIncompleteReadGroup( + gc3.chromosome(), chrIncompleteReadsGroups, bamReader.getChimericReadTracker().getHardFilteredReadIds()); + + finderChr1.processInterChromosomalReadGroups(interChromosomalGroups); + + assertTrue(finderChr1.getFusionCandidates().isEmpty()); + + // now chromosome 2 + FragmentAllocator bamReader2 = new FragmentAllocator(config, geneTransCache, ALT_SJ_COHORT_CACHE, new ResultsWriter(config)); + bamReader2.processReadRecords(gc5, Lists.newArrayList(readPair[1])); + bamReader2.postSliceProcessReads(); + + FusionFinder finderChr2 = fusionTaskManager.createFusionFinder(gc3.chromosome()); + + finderChr2.processNewChimericReadGroups(gc5, bamReader2.getBaseDepth(), bamReader2.getChimericReadTracker().fusionReadGroupMap()); + + chrIncompleteReadsGroups = finderChr2.extractIncompleteReadGroups( + gc5.chromosome(), bamReader2.getChimericReadTracker().getHardFilteredReadIds()); + + interChromosomalGroups = fusionTaskManager.addIncompleteReadGroup( + gc5.chromosome(), chrIncompleteReadsGroups, bamReader2.getChimericReadTracker().getHardFilteredReadIds()); + + finderChr2.processInterChromosomalReadGroups(interChromosomalGroups); + + assertEquals(1, finderChr2.getFusionCandidates().values().stream().mapToInt(x -> x.size()).sum()); + + List fusions = finderChr2.getFusionCandidates().values().stream().findFirst().orElse(Collections.emptyList()); + FusionReadData fusionReadData = fusions.get(0); + + assertEquals(KNOWN, fusionReadData.junctionTypes()[0]); + assertEquals(KNOWN, fusionReadData.junctionTypes()[1]); + } } From 7cb017aa8f1d02136cf0f4a039f69101dfef6229 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Wed, 30 Sep 2026 12:37:35 +1000 Subject: [PATCH 6/8] Isofox: fusion passing filter test for short locals now requires dup-del orientations --- .../com/hartwig/hmftools/isofox/fusion/PassingFusions.java | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/PassingFusions.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/PassingFusions.java index d53f579e34..150c8f44a6 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/PassingFusions.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/PassingFusions.java @@ -149,8 +149,9 @@ public List findPassingFusions(final List allFusions) private static boolean isShortLocalFusion(final FusionData fusion) { - return fusion.Chromosomes[SE_START].equals(fusion.Chromosomes[SE_END]) && - abs(fusion.JunctionPositions[SE_END] - fusion.JunctionPositions[SE_START]) < LOCAL_FUSION_THRESHOLD; + return fusion.Chromosomes[SE_START].equals(fusion.Chromosomes[SE_END]) + && fusion.JunctionOrientations[SE_START] != fusion.JunctionOrientations[SE_END] + && abs(fusion.JunctionPositions[SE_END] - fusion.JunctionPositions[SE_START]) < LOCAL_FUSION_THRESHOLD; } private boolean isPassingFusion(final FusionData fusion) From 7c0c78651c5af0903479e73a8297b11cf9b0db34 Mon Sep 17 00:00:00 2001 From: Charles Shale Date: Mon, 5 Oct 2026 18:01:05 +1100 Subject: [PATCH 7/8] Isofox: allow INV and BND discordant fragment support without requiring soft-clips --- .../hmftools/isofox/fusion/ChimericUtils.java | 21 +++++++++++++++ .../isofox/fusion/FusionConstants.java | 2 ++ .../isofox/fusion/FusionFragmentBuilder.java | 15 ++++++++--- .../isofox/fusion/FusionReadData.java | 26 +++++++++++++++---- 4 files changed, 56 insertions(+), 8 deletions(-) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java index c6b848ae97..6ef2e48931 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java @@ -1,11 +1,14 @@ package com.hartwig.hmftools.isofox.fusion; +import static java.lang.Math.abs; + import static com.hartwig.hmftools.common.fusion.FusionCommon.FS_DOWN; import static com.hartwig.hmftools.common.fusion.FusionCommon.FS_UP; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_END; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_START; import static com.hartwig.hmftools.common.genome.region.Orientation.ORIENT_REV; import static com.hartwig.hmftools.common.genome.region.Orientation.ORIENT_FWD; +import static com.hartwig.hmftools.isofox.fusion.FusionConstants.CHIMERIC_SHORT_INV_MIN_LENGTH; import static com.hartwig.hmftools.isofox.fusion.FusionUtils.aboveJunctionSoftClipThreshold; import static htsjdk.samtools.CigarOperator.N; @@ -22,6 +25,24 @@ public final class ChimericUtils { public static boolean isInversion(final List reads) { + // allow discordant fragments if not short + if(reads.size() == 2) + { + if(reads.stream().anyMatch(x -> x.isSupplementaryAlignment())) + return false; + + Read read1 = reads.get(0); + Read read2 = reads.get(1); + + if(read1.chromosome().equals(read2.chromosome()) && read1.orientation() == read2.orientation() + && abs(read1.alignmentStart() - read2.alignmentStart()) >= CHIMERIC_SHORT_INV_MIN_LENGTH) + { + return true; + } + + return false; + } + // an inversion must a) be same chromosome b) have supplementary alignment c) have same orientations around the chimeric junction if(!reads.stream().anyMatch(x -> x.hasSuppAlignment()) || reads.size() != 3) return false; diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionConstants.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionConstants.java index d383a23575..61053f7e03 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionConstants.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionConstants.java @@ -9,6 +9,8 @@ public class FusionConstants public static final int DEFAULT_HARD_FILTER_MIN_FRAGS = 2; + public static final int CHIMERIC_SHORT_INV_MIN_LENGTH = 2000; + public static final int HIGH_LOG_COUNT = 10000; public static final int FILTER_COHORT_LIMIT_KNOWN = 5; diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java index c363a83707..eb52c6a0e0 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java @@ -315,9 +315,18 @@ private static void setNonJunctionData(final FusionFragment fragment) { fragment.geneCollections()[SE_START] = fragment.geneCollections()[SE_END] = geneCollections.get(0); - // orientation could be set based on the orientations and positions of the reads.. do this when the junction data is set - fragment.orientations()[SE_START] = ORIENT_FWD; - fragment.orientations()[SE_END] = ORIENT_REV; + if(fragment.reads().size() == 2 && fragment.reads().get(0).Orientation == fragment.reads().get(1).Orientation) + { + fragment.setType(DISCORDANT); + fragment.orientations()[SE_END] = fragment.orientations()[SE_START] = fragment.reads().get(0).Orientation; + } + else + { + // orientation could be set based on the orientations and positions of the reads.. do this when the junction data is set + fragment.orientations()[SE_START] = ORIENT_FWD; + fragment.orientations()[SE_END] = ORIENT_REV; + } + return; } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadData.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadData.java index c0308819b3..b3b4f002a2 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadData.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadData.java @@ -462,6 +462,19 @@ public boolean canAddDiscordantFragment(final FusionFragment fragment, int maxFr int impliedFragmentLength = fragment.reads().get(0).ReadBaseLength * 2; + boolean isDelDupFragment = true; + + if(fragment.reads().size() >= 2) + { + FusionRead read1 = fragment.reads().get(0); + FusionRead read2 = fragment.reads().get(1); + + if(!read1.Chromosome.equals(read2.Chromosome) || read1.Orientation == read2.Orientation) + { + isDelDupFragment = false; + } + } + for(int se = SE_START; se <= SE_END; ++se) { final List fragmentRefs = fragment.getTransExonRefs()[se]; @@ -495,11 +508,14 @@ else if(fragment.regionMatchTypes()[se] == INTRON) int fragmentPosition = read.getCoordsBoundary(switchIndex(se)); // cannot be on the wrong side of the junction - if((mJunctionOrientations[se] == 1) == (fragmentPosition > mJunctionPositions[se])) + if((mJunctionOrientations[se] == ORIENT_FWD) == (fragmentPosition > mJunctionPositions[se])) { - // check if a mis-mapping explains the over-hang - if(!softClippedReadSupportsJunction(read, se)) - return false; + if(isDelDupFragment || read.isSoftClipped(se)) + { + // check if a mis-mapping explains the over-hang + if(!softClippedReadSupportsJunction(read, se)) + return false; + } } // of the fusion is unspliced or the fragment is intronic, then measure the genomic distance vs permitted fragment length @@ -512,7 +528,7 @@ else if(fragment.regionMatchTypes()[se] == INTRON) } } - if(impliedFragmentLength > maxFragmentDistance) + if(isDelDupFragment && impliedFragmentLength > maxFragmentDistance) return false; return true; From 6dfc4531266d480e58bb4e317221983c256e675c Mon Sep 17 00:00:00 2001 From: shiv-hartwig Date: Tue, 6 Oct 2026 10:30:51 +1100 Subject: [PATCH 8/8] Isofox: add Fragment class and single-end read support (AUS418) --- .../hmftools/isofox/FragmentAllocator.java | 246 ++++------ .../isofox/adjusts/FragmentSizeCalcs.java | 35 +- .../hmftools/isofox/common/Fragment.java | 223 +++++++++ .../isofox/common/GeneRegionFilters.java | 2 +- .../hartwig/hmftools/isofox/common/Read.java | 7 +- .../isofox/common/ReadTranscriptUtils.java | 26 -- .../expression/ExpressionReadTracker.java | 6 +- .../isofox/fusion/ChimericReadGroup.java | 7 +- .../isofox/fusion/ChimericReadTracker.java | 45 +- .../hmftools/isofox/fusion/ChimericUtils.java | 8 +- .../hmftools/isofox/fusion/FusionFinder.java | 2 +- .../isofox/fusion/FusionFragmentBuilder.java | 21 +- .../isofox/fusion/FusionReadGroup.java | 2 +- .../isofox/novel/AltSpliceJunctionFinder.java | 22 +- .../isofox/novel/RetainedIntronFinder.java | 11 +- .../isofox/novel/SpliceSiteCounter.java | 10 +- .../hartwig/hmftools/isofox/FragmentTest.java | 440 ++++++++++++++++++ .../hmftools/isofox/NovelJunctionsTest.java | 17 +- .../isofox/TransClassificationTest.java | 16 +- .../isofox/fusion/ChimericReadTest.java | 33 +- .../isofox/fusion/FusionDataTest.java | 21 +- 21 files changed, 899 insertions(+), 301 deletions(-) create mode 100644 isofox/src/main/java/com/hartwig/hmftools/isofox/common/Fragment.java create mode 100644 isofox/src/test/java/com/hartwig/hmftools/isofox/FragmentTest.java diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java index f480d1081d..caa5295c87 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/FragmentAllocator.java @@ -26,41 +26,34 @@ import static com.hartwig.hmftools.isofox.common.FragmentType.UNSPLICED; import static com.hartwig.hmftools.isofox.IsofoxFunction.FUSIONS; import static com.hartwig.hmftools.isofox.common.Read.findOverlappingRegions; -import static com.hartwig.hmftools.isofox.common.ReadTranscriptUtils.getUniqueValidRegion; import static com.hartwig.hmftools.isofox.common.ReadTranscriptUtils.markRegionBases; -import static com.hartwig.hmftools.isofox.common.ReadTranscriptUtils.validTranscriptType; import static com.hartwig.hmftools.isofox.common.ReadUtils.consensusDuplicateCount; -import static com.hartwig.hmftools.isofox.common.ReadUtils.trimAdapterBases; -import static com.hartwig.hmftools.isofox.common.RegionMatchType.EXON_INTRON; -import static com.hartwig.hmftools.isofox.common.CommonUtils.deriveCommonRegions; -import static com.hartwig.hmftools.isofox.common.TransMatchType.OTHER_TRANS; import static com.hartwig.hmftools.isofox.common.TransMatchType.SPLICE_JUNCTION; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_END; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_START; -import static com.hartwig.hmftools.isofox.fusion.ChimericUtils.isRealignedFragmentCandidate; import static com.hartwig.hmftools.isofox.results.ResultsWriter.writeReadData; import java.io.BufferedWriter; import java.io.File; import java.io.IOException; import java.util.List; -import java.util.Map; import java.util.Set; import java.util.StringJoiner; import java.util.stream.Collectors; import com.google.common.annotations.VisibleForTesting; import com.google.common.collect.Lists; -import com.google.common.collect.Sets; import com.hartwig.hmftools.common.ensemblcache.EnsemblDataCache; import com.hartwig.hmftools.common.gene.ExonData; import com.hartwig.hmftools.common.gene.GeneData; import com.hartwig.hmftools.common.gene.TranscriptData; import com.hartwig.hmftools.common.bam.BamSlicer; +import com.hartwig.hmftools.common.genome.region.Orientation; import com.hartwig.hmftools.common.region.BaseRegion; import com.hartwig.hmftools.common.region.ChrBaseRegion; import com.hartwig.hmftools.isofox.common.BaseDepth; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.FragmentMatchType; import com.hartwig.hmftools.isofox.common.FragmentTracker; import com.hartwig.hmftools.isofox.common.GeneCollection; @@ -70,11 +63,11 @@ import com.hartwig.hmftools.isofox.common.ReadTranscriptUtils; import com.hartwig.hmftools.isofox.common.RegionReadData; import com.hartwig.hmftools.isofox.common.TransExonRef; -import com.hartwig.hmftools.isofox.common.TransMatchType; import com.hartwig.hmftools.isofox.expression.CategoryCountsData; import com.hartwig.hmftools.isofox.adjusts.GcRatioCounts; import com.hartwig.hmftools.isofox.expression.ExpressionReadTracker; import com.hartwig.hmftools.isofox.fusion.ChimericReadTracker; +import com.hartwig.hmftools.isofox.fusion.ChimericUtils; import com.hartwig.hmftools.isofox.novel.AltSjCohortCache; import com.hartwig.hmftools.isofox.novel.AltSpliceJunctionFinder; import com.hartwig.hmftools.isofox.novel.RetainedIntronFinder; @@ -302,12 +295,18 @@ private void processRead(final Read read) private boolean checkFragmentRead(final Read read) { + if(!read.isReadPaired()) + { + processFragmentReads(new Fragment(read)); + return true; + } + // check if the 2 reads from a fragment exist and if so handle them a pair, returning true Read otherRead = mFragmentReads.checkRead(read); if(otherRead != null) { - processFragmentReads(read, otherRead); + processFragmentReads(new Fragment(read, otherRead)); return true; } @@ -337,49 +336,39 @@ private void handleSupplementaryRead(final Read read) writeChimericReadData(read); } - private void processFragmentReads(final Read read1, final Read read2) + private void processFragmentReads(final Fragment fragment) { - /* process the pair of reads from a fragment: + /* process fragment: - fully outside the gene (due to the buffer used, ignore - read through a gene ie start or end outside - purely intronic - chimeric an inversion or translocation - supporting 1 or more transcripts - - both reads fully with an exon - if exon has only 1 transcript then consider unambiguous - - both reads within 2 exons (including spanning intermediary ones) and/or either exon at the boundary + - fragment fully with an exon - if exon has only 1 transcript then consider unambiguous + - fragment within 2 exons (including spanning intermediary ones) and/or either exon at the boundary - not supporting any transcript - eg alternative splice sites or unspliced reads */ - trimAdapterBases(read1, read2); - - markGeneDataRegions(read1); - markGeneDataRegions(read2); - - int fragmentCount = 1; - - boolean isConsensusRead = read1.isConsensusRead() || read2.isConsensusRead(); + fragment.trimAdapterBases(); - if(isConsensusRead) + for(Read read : fragment.reads()) { - // consensus reads from Redux - these are primaries where duplicate reads are also expected - // to avoid the additional count from the artificially created primary, skip these for any logic which is expression related, - // (note: only expression uses duplicates) - int duplicateCount = consensusDuplicateCount(read1.bamRecord()); - fragmentCount += duplicateCount; + markGeneDataRegions(read); } - int numLoci = min(read1.numLoci(), read2.numLoci()); + int fragmentCount = fragment.fragmentCount(); + int numLoci = fragment.minNumLoci(); boolean isMultiMapped = numLoci > 1; - List altLoci = read1.numLoci() <= read2.numLoci() ? read1.altLoci() : read2.altLoci(); + List altLoci = fragment.altLoci(); - boolean isChimeric = mChimericReads.isChimeric(read1, read2, isMultiMapped); + boolean isChimeric = mChimericReads.isChimeric(fragment, isMultiMapped); if(mStatsOnly) { if(isChimeric) mCurrentGenes.addCount(CHIMERIC, 1); - else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty()) + else if(fragment.isFullyIntronic()) mCurrentGenes.addCount(UNSPLICED, fragmentCount); else mCurrentGenes.addCount(TRANS_SUPPORTING, fragmentCount); @@ -387,7 +376,7 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() return; } - List commonMappings = deriveCommonRegions(read1.getMappedRegionCoords(), read2.getMappedRegionCoords()); + List commonMappings = fragment.mergedMappings(); mBaseDepth.processRead(commonMappings); @@ -398,7 +387,7 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() if(!isMultiMapped) { if(mChimericReads.enabled()) - mChimericReads.addChimericReadPair(read1, read2); + mChimericReads.addChimericFragment(fragment); else mCurrentGenes.addCount(CHIMERIC, 1); } @@ -406,11 +395,13 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() if(mReadDataWriter != null && mConfig.writeType(WriteType.READ)) { List overlapGenes = mCurrentGenes.findGenesCoveringRange( - min(read1.alignmentStart(), read2.alignmentStart()), - max(read1.alignmentEnd(), read2.alignmentEnd()), true); + fragment.minAlignmentStart(), + fragment.maxAlignmentEnd(), true); - writeReadData(mReadDataWriter, overlapGenes, read1, CHIMERIC, 0); - writeReadData(mReadDataWriter, overlapGenes, read2, CHIMERIC, 0); + for(Read read : fragment.reads()) + { + writeReadData(mReadDataWriter, overlapGenes, read, CHIMERIC, 0); + } } return; @@ -419,9 +410,9 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() if(mRunFusions) { // reads with sufficient soft-clipping and not mapped to an adjacent region are candidates for fusion re-alignment - if(isRealignedFragmentCandidate(read1) || isRealignedFragmentCandidate(read2)) + if(fragment.reads().stream().anyMatch(ChimericUtils::isRealignedFragmentCandidate)) { - mChimericReads.addRealignmentCandidates(read1, read2); + mChimericReads.addRealignmentCandidates(fragment); } if(mFusionsOnly) @@ -429,28 +420,24 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() } if(numLoci > 1 && altLoci != null && mMultiMapLociWriter != null) - recordMultiMapLoci(read1, read2, altLoci); + recordMultiMapLoci(fragment, altLoci); - int readPosMin = min(read1.alignmentStart(), read2.alignmentStart()); - int readPosMax = max(read1.alignmentEnd(), read2.alignmentEnd()); + int readPosMin = fragment.minAlignmentStart(); + int readPosMax = fragment.maxAlignmentEnd(); List overlapGenes = mCurrentGenes.findGenesCoveringRange(readPosMin, readPosMax, true); - if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty()) + if(fragment.isFullyIntronic()) { // fully intronic read in every transcript and gene - processIntronicReads(overlapGenes, read1, read2, fragmentCount, isMultiMapped); + processIntronicReads(overlapGenes, fragment, fragmentCount, isMultiMapped); return; } - Map firstReadTransTypes = read1.getTranscriptClassifications(); - Map secondReadTransTypes = read2.getTranscriptClassifications(); - // first find valid transcripts in both reads List validTranscripts = Lists.newArrayList(); - Set invalidTranscripts = Sets.newHashSet(); - List validRegions = getUniqueValidRegion(read1, read2); + List validRegions = fragment.uniqueValidRegions(); if(mConfig.RunValidations) { @@ -466,40 +453,18 @@ else if(read1.getMappedRegions().isEmpty() && read2.getMappedRegions().isEmpty() // track splice site info if(mConfig.writeType(SPLICE_SITE)) { - mSpliceSiteCounter.registerSpliceSiteSupport( - read1.getMappedRegionCoords(), read2.getMappedRegionCoords(), mCurrentGenes.getExonRegions()); + mSpliceSiteCounter.registerSpliceSiteSupport(fragment, mCurrentGenes.getExonRegions()); } - for(Map.Entry entry : firstReadTransTypes.entrySet()) + for(int transId : fragment.validTypeTranscripts()) { - int transId = entry.getKey(); - - if(validTranscriptType(entry.getValue()) && secondReadTransTypes.containsKey(transId) - && validTranscriptType(secondReadTransTypes.get(transId))) - { - int calcFragmentLength = calcFragmentLength(transId, read1, read2); - boolean validFragmentLength = calcFragmentLength > 0 && calcFragmentLength <= mConfig.MaxFragmentLength; + int calcFragmentLength = calcFragmentLength(transId, fragment); - if(validFragmentLength) - { - validTranscripts.add(transId); - } - else - { - invalidTranscripts.add(transId); - } - } - else - { - invalidTranscripts.add(transId); - } + if(calcFragmentLength > 0 && calcFragmentLength <= mConfig.MaxFragmentLength) + validTranscripts.add(transId); } - for(Integer transId : secondReadTransTypes.keySet()) - { - if(!validTranscripts.contains(transId)) - invalidTranscripts.add(transId); - } + Set invalidTranscripts = fragment.invalidTranscripts(validTranscripts); FragmentType fragmentType = UNSPLICED; @@ -509,14 +474,14 @@ && validTranscriptType(secondReadTransTypes.get(transId))) // no valid transcripts but record against the gene further information about these reads boolean checkRetainedIntrons = false; - if(read1.containsSplit() || read2.containsSplit()) + if(fragment.containsSplit()) { fragmentType = ALT; if(mAltSpliceJunctionFinder.enabled()) { mAltSpliceJunctionFinder.evaluateFragmentReads( - overlapGenes, read1, read2, invalidTranscripts.stream().collect(Collectors.toList())); + overlapGenes, fragment.reads(), invalidTranscripts.stream().collect(Collectors.toList())); } checkRetainedIntrons = true; @@ -526,23 +491,7 @@ && validTranscriptType(secondReadTransTypes.get(transId))) // look for alternative splicing from long reads involving more than one region and not spanning into an intron for(int transId : invalidTranscripts) { - List regions = read1.getMappedRegions().entrySet().stream() - .filter(x -> x.getKey().hasTransId(transId)) - .filter(x -> x.getValue() != EXON_INTRON) - .map(x -> x.getKey()).collect(Collectors.toList());; - - List regions2 = read2.getMappedRegions().entrySet().stream() - .filter(x -> x.getKey().hasTransId(transId)) - .filter(x -> x.getValue() != EXON_INTRON) - .map(x -> x.getKey()).collect(Collectors.toList()); - - for(RegionReadData region : regions2) - { - if(!regions.contains(region)) - regions.add(region); - } - - if(regions.size() > 1) + if(fragment.spansMultipleRegions(transId)) { fragmentType = ALT; break; @@ -553,7 +502,7 @@ && validTranscriptType(secondReadTransTypes.get(transId))) } if(checkRetainedIntrons && mRetainedIntronFinder.enabled()) - mRetainedIntronFinder.evaluateFragmentReads(read1, read2); + mRetainedIntronFinder.evaluateFragmentReads(fragment); if(fragmentType == UNSPLICED) { @@ -566,15 +515,7 @@ && validTranscriptType(secondReadTransTypes.get(transId))) fragmentType = TRANS_SUPPORTING; // first mark any invalid trans as 'other' meaning it doesn't require any further classification since a valid trans exists - firstReadTransTypes.entrySet().stream() - .filter(x -> validTranscriptType(x.getValue())) - .filter(x -> !validTranscripts.contains(x.getKey())) - .forEach(x -> x.setValue(OTHER_TRANS)); - - secondReadTransTypes.entrySet().stream() - .filter(x -> validTranscriptType(x.getValue())) - .filter(x -> !validTranscripts.contains(x.getKey())) - .forEach(x -> x.setValue(OTHER_TRANS)); + fragment.setOtherTranscripts(validTranscripts); if(mConfig.RunValidations) { @@ -603,13 +544,16 @@ && validTranscriptType(secondReadTransTypes.get(transId))) if(supportedGeneIsForward == null) { - supportedGeneIsForward = findGeneStrand(read1, validTranscripts); + for(Read read : fragment.reads()) + { + supportedGeneIsForward = findGeneStrand(read, validTranscripts); - if(supportedGeneIsForward == null) - supportedGeneIsForward = findGeneStrand(read2, validTranscripts); + if (supportedGeneIsForward != null) + break; + } } - if(read1.getTranscriptClassification(transId) == SPLICE_JUNCTION || read2.getTranscriptClassification(transId) == SPLICE_JUNCTION) + if(fragment.hasTranscriptClassification(transId, SPLICE_JUNCTION)) { transMatchType = FragmentMatchType.SPLICED; comboTransMatchType = FragmentMatchType.SPLICED; @@ -630,40 +574,13 @@ else if(regionCount > 1) mCurrentGenes.addTranscriptReadMatch(transId, isUniqueTrans, transMatchType); // separately record discordant reads spanning 2+ exons - if(!read1.containsSplit() && !read2.containsSplit()) + if(!fragment.containsSplit() && fragment.readsInDifferentExons(transId)) { - boolean hasExonRankMatch = false; - - for(RegionReadData region1 : read1.getMappedRegions().keySet()) - { - TransExonRef transExonRef1 = region1.getTransExonRefs().stream().filter(x -> x.TransId == transId).findFirst().orElse(null); - - if(transExonRef1 == null) - continue; - - for(RegionReadData region2 : read2.getMappedRegions().keySet()) - { - TransExonRef transExonRef2 = region2.getTransExonRefs().stream().filter(x -> x.TransId == transId).findFirst().orElse(null); - - if(transExonRef2 != null && transExonRef1.ExonRank == transExonRef2.ExonRank) - { - hasExonRankMatch = true; - break; - } - } - - if(hasExonRankMatch) - break; - } - - if(!hasExonRankMatch) // region1 != region2 && region1.getExonRank(transId) != region2.getExonRank(transId) - { - mCurrentGenes.addTranscriptReadMatch(transId, DISCORDANT); - } + mCurrentGenes.addTranscriptReadMatch(transId, DISCORDANT); } // keep track of which regions have been allocated from this fragment as a whole, so not counting each read separately - mExpressionReadTracker.processValidTranscript(transId, List.of(read1, read2), isUniqueTrans); + mExpressionReadTracker.processValidTranscript(transId, fragment.reads(), isUniqueTrans); } mExpressionReadTracker.processUnsplicedGenes( @@ -672,12 +589,11 @@ else if(regionCount > 1) if(supportedGeneIsForward != null) { // track fragment strandedness - boolean firstIsForward = read1.isFirstOfPair() ? !read1.isReadReversed() : !read2.isReadReversed(); - boolean secondIsForward = !read1.isFirstOfPair() ? !read1.isReadReversed() : !read2.isReadReversed(); + Orientation fragmentOrientation = fragment.orientation(); - if(firstIsForward != secondIsForward) + if(fragmentOrientation != null) { - if(firstIsForward == supportedGeneIsForward) + if(fragmentOrientation.isForward() == supportedGeneIsForward) mCurrentGenes.addCount(FORWARD_STRAND, fragmentCount); else mCurrentGenes.addCount(REVERSE_STRAND, fragmentCount); @@ -689,8 +605,10 @@ else if(regionCount > 1) if(mReadDataWriter != null && mConfig.writeType(WriteType.READ)) { - writeReadData(mReadDataWriter, overlapGenes, read1, fragmentType, validTranscripts.size()); - writeReadData(mReadDataWriter, overlapGenes, read2, fragmentType, validTranscripts.size()); + for(Read read : fragment.reads()) + { + writeReadData(mReadDataWriter, overlapGenes, read, fragmentType, validTranscripts.size()); + } } } @@ -715,13 +633,13 @@ private Boolean findGeneStrand(final Read read, final List transcripts) return null; } - private int calcFragmentLength(int transId, final Read read1, final Read read2) + private int calcFragmentLength(int transId, final Fragment fragment) { TranscriptData transData = mCurrentGenes.getTranscripts().stream().filter(x -> x.TransId == transId).findFirst().orElse(null); if(transData == null) return -1; - return ReadTranscriptUtils.calcFragmentLength(transData, read1, read2); + return ReadTranscriptUtils.calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); } private boolean reachedGeneReadLimit() @@ -756,19 +674,19 @@ private boolean altOverlapsExon(final GeneData gene, final ChrBaseRegion altRegi return false; } - private void recordMultiMapLoci(final Read read1, final Read read2, final List altLoci) + private void recordMultiMapLoci(final Fragment fragment, final List altLoci) { if(mMultiMapLociWriter == null) return; // record a multi-mapped fragment's primary alignment plus each XA alternate locus (opt-in WriteType.MULTI_MAP_LOCI); // InGeneCollection flags whether the locus falls within the gene collection currently being processed - int fragStart = min(read1.alignmentStart(), read2.alignmentStart()); - int fragEnd = max(read1.alignmentEnd(), read2.alignmentEnd()); - boolean primarySpliced = read1.containsSplit() || read2.containsSplit(); + int fragStart = fragment.minAlignmentStart(); + int fragEnd = fragment.maxAlignmentEnd(); + boolean primarySpliced = fragment.containsSplit(); writeMultiMapLocus( - mMultiMapLociWriter, mCurrentGenes.id(), read1.id(), "PRIMARY", read1.chromosome(), fragStart, fragEnd, + mMultiMapLociWriter, mCurrentGenes.id(), fragment.id(), "PRIMARY", fragment.chromosome(), fragStart, fragEnd, primarySpliced, mCurrentGenes.geneNames(), true); int[] bounds = mCurrentGenes.regionBounds(); @@ -779,7 +697,7 @@ private void recordMultiMapLoci(final Read read1, final Read read2, final List genes, final Read read1, final Read read2, int fragmentCount, boolean multiMapped) + final List genes, final Fragment fragment, int fragmentCount, boolean multiMapped) { - if(read1.containsSplit() || read2.containsSplit()) + if(fragment.containsSplit()) { mCurrentGenes.addCount(ALT, 1); // does not count duplicates since not expression related if(mAltSpliceJunctionFinder.enabled()) - mAltSpliceJunctionFinder.evaluateFragmentReads(genes, read1, read2, Lists.newArrayList()); + mAltSpliceJunctionFinder.evaluateFragmentReads(genes, fragment.reads(), Lists.newArrayList()); return; } - mExpressionReadTracker.processIntronicReads(genes, read1, read2, fragmentCount, multiMapped); + mExpressionReadTracker.processIntronicReads(genes, fragment, fragmentCount, multiMapped); mCurrentGenes.addCount(UNSPLICED, fragmentCount); if(mReadDataWriter != null && mConfig.writeType(WriteType.READ)) { - writeReadData(mReadDataWriter, genes, read1, UNSPLICED, 0); - writeReadData(mReadDataWriter, genes, read2, UNSPLICED, 0); + for(Read read : fragment.reads()) + { + writeReadData(mReadDataWriter, genes, read, UNSPLICED, 0); + } } } @@ -902,7 +822,7 @@ else if(reads.size() == 3) int readPosMax = max(read1.alignmentEnd(), read2.alignmentEnd()); List overlapGenes = mCurrentGenes.findGenesCoveringRange(readPosMin, readPosMax, false); - mAltSpliceJunctionFinder.evaluateFragmentReads(overlapGenes, read1, read2, invalidTrans); + mAltSpliceJunctionFinder.evaluateFragmentReads(overlapGenes, List.of(read1, read2), invalidTrans); } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/adjusts/FragmentSizeCalcs.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/adjusts/FragmentSizeCalcs.java index f2fd811a8e..c6079b04b7 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/adjusts/FragmentSizeCalcs.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/adjusts/FragmentSizeCalcs.java @@ -1,19 +1,19 @@ package com.hartwig.hmftools.isofox.adjusts; -import static java.lang.Math.abs; import static java.lang.Math.max; import static java.lang.Math.min; import static java.lang.Math.round; +import static com.hartwig.hmftools.common.bam.SamRecordUtils.inferredInsertSizeAbs; import static com.hartwig.hmftools.common.bam.SamRecordUtils.mateNegativeStrand; import static com.hartwig.hmftools.common.region.BaseRegion.positionWithin; import static com.hartwig.hmftools.common.region.BaseRegion.positionsOverlap; -import static com.hartwig.hmftools.common.utils.file.FileDelimiters.ITEM_DELIM; -import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.closeBufferedWriter; -import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.createBufferedWriter; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_END; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_PAIR; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_START; +import static com.hartwig.hmftools.common.utils.file.FileDelimiters.ITEM_DELIM; +import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.closeBufferedWriter; +import static com.hartwig.hmftools.common.utils.file.FileWriterUtils.createBufferedWriter; import static com.hartwig.hmftools.isofox.ChromosomeTaskExecutor.findNextOverlappingGenes; import static com.hartwig.hmftools.isofox.IsofoxConfig.ISF_LOGGER; import static com.hartwig.hmftools.isofox.IsofoxConstants.SINGLE_MAP_QUALITY; @@ -30,7 +30,6 @@ import com.hartwig.hmftools.common.gene.GeneData; import com.hartwig.hmftools.common.gene.TranscriptData; import com.hartwig.hmftools.common.perf.PerformanceCounter; -import com.hartwig.hmftools.common.region.BaseRegion; import com.hartwig.hmftools.common.region.ChrBaseRegion; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.WriteType; @@ -259,10 +258,13 @@ private void processBamRead(@NotNull final SAMRecord read) mMaxReadLength = max(mMaxReadLength, read.getReadLength()); - final SAMRecord otherRead = (SAMRecord) mFragmentTracker.checkRead(read.getReadName(), read); + if(read.getReadPairedFlag()) + { + final SAMRecord otherRead = (SAMRecord) mFragmentTracker.checkRead(read.getReadName(), read); - if(otherRead == null) - return; + if(otherRead == null) + return; + } addFragmentLength(read, mFragmentLengths); @@ -274,18 +276,26 @@ private void processBamRead(@NotNull final SAMRecord read) private boolean isCandidateRecord(final SAMRecord record) { - int fragmentLength = abs(record.getInferredInsertSize()); + int fragmentLength = inferredInsertSizeAbs(record); if(fragmentLength > MAX_FRAGMENT_LENGTH) return false; // ignore translocations and inversions - if(!record.getMateReferenceName().equals(record.getReferenceName()) || mateNegativeStrand(record) == record.getReadNegativeStrandFlag()) + if(record.getReadPairedFlag() + && (!record.getMateReferenceName().equals(record.getReferenceName()) || + mateNegativeStrand(record) == record.getReadNegativeStrandFlag())) + { return false; + } // ignore split and soft-clipped reads above the read length if(record.getCigar().containsOperator(CigarOperator.N) || !record.getCigar().containsOperator(CigarOperator.M)) return false; + // unpaired reads must be fully aligned + if(!record.getReadPairedFlag() && record.getCigar().containsOperator(CigarOperator.S)) + return false; + int readLength = max(mMaxReadLength, mConfig.ReadLength); if(readLength > 0 && fragmentLength > readLength && record.getCigar().containsOperator(CigarOperator.S)) @@ -296,8 +306,7 @@ private boolean isCandidateRecord(final SAMRecord record) } // both reads must fall in the current gene - int otherStartPos = record.getMateAlignmentStart(); - if(!positionWithin(otherStartPos, mCurrentGenesRange[SE_START], mCurrentGenesRange[SE_END])) + if(record.getReadPairedFlag() && !positionWithin(record.getMateAlignmentStart(), mCurrentGenesRange[SE_START], mCurrentGenesRange[SE_END])) return false; // reads cannot cover any part of an exon @@ -315,7 +324,7 @@ private boolean isCandidateRecord(final SAMRecord record) private void addFragmentLength(final SAMRecord record, final List fragmentLengths) { - int fragmentLength = getLengthBucket(abs(record.getInferredInsertSize())); + int fragmentLength = getLengthBucket(inferredInsertSizeAbs(record)); if(fragmentLength == 0) return; diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Fragment.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Fragment.java new file mode 100644 index 0000000000..a24036f34e --- /dev/null +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Fragment.java @@ -0,0 +1,223 @@ +package com.hartwig.hmftools.isofox.common; + +import static com.hartwig.hmftools.isofox.common.CommonUtils.deriveCommonRegions; +import static com.hartwig.hmftools.isofox.common.ReadTranscriptUtils.validTranscriptType; +import static com.hartwig.hmftools.isofox.common.ReadUtils.consensusDuplicateCount; +import static com.hartwig.hmftools.isofox.common.RegionMatchType.EXON_INTRON; +import static com.hartwig.hmftools.isofox.common.RegionMatchType.validExonMatch; +import static com.hartwig.hmftools.isofox.common.RegionReadData.NO_EXON; +import static com.hartwig.hmftools.isofox.common.TransMatchType.OTHER_TRANS; + +import java.util.List; +import java.util.Map; +import java.util.Set; + +import com.google.common.collect.Lists; +import com.google.common.collect.Sets; +import com.hartwig.hmftools.common.genome.region.Orientation; +import com.hartwig.hmftools.common.region.BaseRegion; + +public class Fragment +{ + private List mReads; + + public Fragment(final Read read) + { + mReads = List.of(read); + } + + public Fragment(final Read read1, final Read read2) + { + mReads = List.of(read1, read2); + } + + public List reads() + { + return mReads; + } + + public String id() { return mReads.get(0).id(); } + + public String chromosome() { return mReads.get(0).chromosome(); } + + // consensus reads from Redux - these are primaries where duplicate reads are also expected + // to avoid the additional count from the artificially created primary, skip these for any logic which is expression related, + // (note: only expression uses duplicates) + public int fragmentCount() + { + return 1 + consensusDuplicateCount(mReads.get(0).bamRecord()); + } + + public int minNumLoci() + { + return mReads.stream().mapToInt(Read::numLoci).min().getAsInt(); + } + + public List altLoci() + { + Read minLociRead = mReads.get(0); + + for(Read read : mReads) + { + if(read.numLoci() < minLociRead.numLoci()) + minLociRead = read; + } + + return minLociRead.altLoci(); + } + + public boolean isFullyIntronic() + { + return mReads.stream().allMatch(x -> x.getMappedRegions().isEmpty()); + } + + public boolean containsSplit() + { + return mReads.stream().anyMatch(Read::containsSplit); + } + + public int minAlignmentStart() + { + return mReads.stream().mapToInt(Read::alignmentStart).min().getAsInt(); + } + + public int maxAlignmentEnd() + { + return mReads.stream().mapToInt(Read::alignmentEnd).max().getAsInt(); + } + + public List mergedMappings() + { + if(mReads.size() == 1) + return mReads.get(0).getMappedRegionCoords(); + + return deriveCommonRegions(mReads.get(0).getMappedRegionCoords(), mReads.get(1).getMappedRegionCoords()); + } + + public List uniqueValidRegions() + { + List validRegions = Lists.newArrayList(); + + for(Read read : mReads) + { + for(Map.Entry entry : read.getMappedRegions().entrySet()) + { + if(validExonMatch(entry.getValue()) && !validRegions.contains(entry.getKey())) + validRegions.add(entry.getKey()); + } + } + + return validRegions; + } + + public boolean spansMultipleRegions(int transId) + { + List regions = Lists.newArrayList(); + + for(Read read : mReads) + { + for(Map.Entry entry : read.getMappedRegions().entrySet()) + { + RegionReadData region = entry.getKey(); + + if(region.hasTransId(transId) && entry.getValue() != EXON_INTRON && !regions.contains(region)) + { + regions.add(region); + + if(regions.size() > 1) + return true; + } + } + } + + return false; + } + + public List validTypeTranscripts() + { + List transIds = Lists.newArrayList(); + + for(int transId : mReads.get(0).getTranscriptClassifications().keySet()) + { + if(mReads.stream().allMatch(x -> validTranscriptType(x.getTranscriptClassification(transId)))) + transIds.add(transId); + } + + return transIds; + } + + public Set invalidTranscripts(final List validTranscripts) + { + Set transIds = Sets.newHashSet(); + + for(Read read : mReads) + { + for(int transId : read.getTranscriptClassifications().keySet()) + { + if(!validTranscripts.contains(transId)) + transIds.add(transId); + } + } + + return transIds; + } + + public void setOtherTranscripts(final List validTranscripts) + { + for(Read read : mReads) + { + read.getTranscriptClassifications().entrySet().stream() + .filter(x -> validTranscriptType(x.getValue())) + .filter(x -> !validTranscripts.contains(x.getKey())) + .forEach(x -> x.setValue(OTHER_TRANS)); + } + } + + public boolean hasTranscriptClassification(int transId, final TransMatchType type) + { + return mReads.stream().anyMatch(x -> x.getTranscriptClassification(transId) == type); + } + + public void trimAdapterBases() + { + if(mReads.size() == 2) + ReadUtils.trimAdapterBases(mReads.get(0), mReads.get(1)); + } + + public boolean readsInDifferentExons(int transId) + { + if(mReads.size() < 2) + return false; + + for(RegionReadData region1 : mReads.get(0).getMappedRegions().keySet()) + { + int exonRank = region1.getExonRank(transId); + + if(exonRank == NO_EXON) + continue; + + for(RegionReadData region2 : mReads.get(1).getMappedRegions().keySet()) + { + if(region2.getExonRank(transId) == exonRank) + return false; + } + } + + return true; + } + + public Orientation orientation() + { + Read read1 = mReads.get(0); + + if(mReads.size() == 1) + return read1.orientation(); + + Read read2 = mReads.get(1); + + if(read1.orientation() == read2.orientation()) + return null; + + return read1.isFirstOfPair() ? read1.orientation() : read2.orientation(); + } +} \ No newline at end of file diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/GeneRegionFilters.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/GeneRegionFilters.java index 9ab98bd9c0..0edd5a1396 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/GeneRegionFilters.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/GeneRegionFilters.java @@ -81,7 +81,7 @@ public boolean skipRead(final SAMRecord read, boolean checkMateAndSupp) if(checkMateAndSupp) { // simple, non-cigar aware read end - if(!read.getMateUnmappedFlag()) + if(read.getReadPairedFlag() && !read.getMateUnmappedFlag()) { int mateReadStart = read.getMateAlignmentStart(); if(skipRead(read.getMateReferenceName(), mateReadStart, mateReadStart, true)) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Read.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Read.java index d2d84f4b72..f49aeb850a 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Read.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/Read.java @@ -10,6 +10,7 @@ import static com.hartwig.hmftools.common.bam.SamRecordUtils.firstInPair; import static com.hartwig.hmftools.common.bam.SamRecordUtils.inferredInsertSize; import static com.hartwig.hmftools.common.bam.SamRecordUtils.mateNegativeStrand; +import static com.hartwig.hmftools.common.bam.SamRecordUtils.mateUnmapped; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_END; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_PAIR; import static com.hartwig.hmftools.common.sv.StartEndIterator.SE_START; @@ -149,10 +150,10 @@ private void setBoundaries() public boolean isReadReversed() { return mRecord.getReadNegativeStrandFlag(); } public boolean isFirstOfPair() { return firstInPair(mRecord); } public boolean isDuplicate() { return mRecord.getDuplicateReadFlag(); } - public boolean isTranslocation() { return !chromosome().equals(mateChromosome()); } + public boolean isTranslocation() { return isReadPaired() && !chromosome().equals(mateChromosome()); } public boolean isMateNegStrand() { return mateNegativeStrand(mRecord); } - public boolean isMateUnmapped() { return mRecord.getMateUnmappedFlag(); } - public boolean isInversion() { return isReadReversed() == isMateNegStrand(); } + public boolean isMateUnmapped() { return mateUnmapped(mRecord); } + public boolean isInversion() { return isReadPaired() && isReadReversed() == isMateNegStrand(); } public boolean isSupplementaryAlignment() { return mRecord.getSupplementaryAlignmentFlag(); } public SupplementaryReadData supplementaryData() { return mSupplementaryData; } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/ReadTranscriptUtils.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/ReadTranscriptUtils.java index 5484c8e0ab..54247a6312 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/common/ReadTranscriptUtils.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/common/ReadTranscriptUtils.java @@ -347,25 +347,6 @@ private static boolean shortClipWithinExon(final Read read, int se, final List= transcriptBoundary : readPositionBoundary + clipLength <= transcriptBoundary; } - public static List getUniqueValidRegion(final Read read1, final Read read2) - { - List regions = read1.getMappedRegions().entrySet().stream() - .filter(x -> validExonMatch(x.getValue())) - .map(x -> x.getKey()).collect(Collectors.toList()); - - List regions2 = read2.getMappedRegions().entrySet().stream() - .filter(x -> validExonMatch(x.getValue())) - .map(x -> x.getKey()).collect(Collectors.toList()); - - for(RegionReadData region : regions2) - { - if(!regions.contains(region)) - regions.add(region); - } - - return regions; - } - public static boolean validTranscriptType(TransMatchType transType) { return transType == EXONIC || transType == SPLICE_JUNCTION; @@ -448,13 +429,6 @@ public static void markRegionBases(final List readCoords, final Regi } } - public static int calcFragmentLength(final TranscriptData transData, final Read read1, final Read read2) - { - int minReadPos = min(read1.alignmentStart(), read2.alignmentStart()); - int maxReadPos = max(read1.alignmentEnd(), read2.alignmentEnd()); - return calcFragmentLength(transData, minReadPos, maxReadPos); - } - public static int calcFragmentLength(final TranscriptData transData, final int minReadPos, final int maxReadPos) { // calculate fragment length within this transcript assuming it has been spliced diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/expression/ExpressionReadTracker.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/expression/ExpressionReadTracker.java index 578e0a05e3..2b297c95ec 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/expression/ExpressionReadTracker.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/expression/ExpressionReadTracker.java @@ -6,7 +6,6 @@ import static com.hartwig.hmftools.isofox.adjusts.GcRatioCounts.calcGcRatioFromReadRegions; import static com.hartwig.hmftools.isofox.common.FragmentType.MULTI_MAPPED; import static com.hartwig.hmftools.isofox.common.RegionMatchType.validExonMatch; -import static com.hartwig.hmftools.isofox.common.CommonUtils.deriveCommonRegions; import static com.hartwig.hmftools.isofox.common.TransMatchType.SPLICE_JUNCTION; import java.util.List; @@ -16,6 +15,7 @@ import com.hartwig.hmftools.common.region.BaseRegion; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.adjusts.GcRatioCounts; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.FragmentMatchType; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.GeneReadData; @@ -89,7 +89,7 @@ public void processUnsplicedGenes( } public void processIntronicReads( - final List genes, final Read read1, final Read read2, int fragmentCount, boolean multiMapped) + final List genes, final Fragment fragment, int fragmentCount, boolean multiMapped) { if(!mEnabled) return; @@ -100,7 +100,7 @@ public void processIntronicReads( { CategoryCountsData catCounts = getCategoryCountsData(Lists.newArrayList(), unsplicedGeneIds); - List readRegions = deriveCommonRegions(read1.getMappedRegionCoords(), read2.getMappedRegionCoords()); + List readRegions = fragment.mergedMappings(); addGcCounts(catCounts, readRegions, fragmentCount, multiMapped); } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadGroup.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadGroup.java index 828a54404b..63335ed2d1 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadGroup.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadGroup.java @@ -3,6 +3,7 @@ import java.util.List; import com.google.common.collect.Lists; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.Read; public class ChimericReadGroup @@ -18,11 +19,9 @@ public ChimericReadGroup(final Read read) mComplete = readGroupComplete(); } - public ChimericReadGroup(final Read read1, final Read read2) + public ChimericReadGroup(final Fragment fragment) { - mReads = Lists.newArrayListWithCapacity(2); - mReads.add(read1); - mReads.add(read2); + mReads = Lists.newArrayList(fragment.reads()); mComplete = readGroupComplete(); } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTracker.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTracker.java index 4b2f7f67c3..8b04c57097 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTracker.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTracker.java @@ -40,6 +40,7 @@ import com.hartwig.hmftools.common.bam.SupplementaryReadData; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.common.BaseDepth; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.FragmentTracker; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.Read; @@ -123,14 +124,14 @@ public JunctionRacFragments extractJunctionRacFragments() public void setChimericPosDataWriter(final BufferedWriter writer) { mChimericPosDataWriter = writer; } - public boolean isChimeric(final Read read1, final Read read2, boolean isMultiMapped) + public boolean isChimeric(final Fragment fragment, boolean isMultiMapped) { - if(read1.isChimeric() || read2.isChimeric() || !read1.withinGeneCollection() || !read2.withinGeneCollection()) + if(fragment.reads().stream().anyMatch(x->x.isChimeric() || !x.withinGeneCollection())) return true; - if(!isMultiMapped && enabled() && (read1.containsSplit() || read2.containsSplit())) + if(!isMultiMapped && enabled() && fragment.containsSplit()) { - return setHasMultipleKnownSpliceGenes(Lists.newArrayList(read1, read2), mKnownPairGeneIds); + return setHasMultipleKnownSpliceGenes(fragment.reads(), mKnownPairGeneIds); } return false; @@ -206,45 +207,47 @@ private void clear(boolean full) } } - public void addRealignmentCandidates(final Read read1, final Read read2) + public void addRealignmentCandidates(final Fragment fragment) { - if(read1.isDuplicate() || read2.isDuplicate()) + if(fragment.reads().stream().anyMatch(Read::isDuplicate)) return; - mCandidateRealignedGroups.add(new ChimericReadGroup(read1, read2)); + mCandidateRealignedGroups.add(new ChimericReadGroup(fragment)); } - public void addChimericReadPair(final Read read1, final Read read2) + public void addChimericFragment(final Fragment fragment) { - if(read1.isDuplicate() || read2.isDuplicate() || inImmuneRegion(read1) || inImmuneRegion(read2)) + if(fragment.reads().stream().anyMatch(x -> x.isDuplicate() || inImmuneRegion(x))) return; // populate transcript info for intronic reads since it will be used in fusion matching - addIntronicTranscriptData(read1); - addIntronicTranscriptData(read2); + for(Read read : fragment.reads()) + { + addIntronicTranscriptData(read); + } // add the pair when it's clear there aren't others with the same ID in the map - if(mConfig.RunValidations && mChimericReadMap.containsKey(read1.id())) + if(mConfig.RunValidations && mChimericReadMap.containsKey(fragment.id())) { // shouldn't occur - ISF_LOGGER.error("overriding chimeric read({})", read1.id()); + ISF_LOGGER.error("overriding chimeric read({})", fragment.id()); - ChimericReadGroup existingGroup = mChimericReadMap.get(read1.id()); + ChimericReadGroup existingGroup = mChimericReadMap.get(fragment.id()); for(Read read : existingGroup.reads()) { ISF_LOGGER.error("existing read: {}", read); } - ISF_LOGGER.error("new read: {}", read1); - ISF_LOGGER.error("new read: {}", read2); - - existingGroup.addRead(read1); - existingGroup.addRead(read2); + for(Read read : fragment.reads()) + { + ISF_LOGGER.error("new read: {}", read); + existingGroup.addRead(read); + } } else { - ChimericReadGroup readGroup = new ChimericReadGroup(read1, read2); + ChimericReadGroup readGroup = new ChimericReadGroup(fragment); if(mChimericPosDataWriter != null) { @@ -253,7 +256,7 @@ public void addChimericReadPair(final Read read1, final Read read2) if(!mConfig.Fusions.WriteChimericOnly) { - mChimericReadMap.put(read1.id(), readGroup); + mChimericReadMap.put(fragment.id(), readGroup); } } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java index 6ef2e48931..2cee283adf 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/ChimericUtils.java @@ -26,11 +26,8 @@ public final class ChimericUtils public static boolean isInversion(final List reads) { // allow discordant fragments if not short - if(reads.size() == 2) + if(reads.size() == 2 && reads.stream().noneMatch(x -> x.isSupplementaryAlignment())) { - if(reads.stream().anyMatch(x -> x.isSupplementaryAlignment())) - return false; - Read read1 = reads.get(0); Read read2 = reads.get(1); @@ -44,7 +41,8 @@ && abs(read1.alignmentStart() - read2.alignmentStart()) >= CHIMERIC_SHORT_INV_MI } // an inversion must a) be same chromosome b) have supplementary alignment c) have same orientations around the chimeric junction - if(!reads.stream().anyMatch(x -> x.hasSuppAlignment()) || reads.size() != 3) + if(reads.stream().noneMatch(x -> x.hasSuppAlignment()) || + reads.stream().filter(x -> x.isSupplementaryAlignment()).count() != 1) return false; byte existingOrient = 0; diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFinder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFinder.java index 432322e35c..9960b4c784 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFinder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFinder.java @@ -305,7 +305,7 @@ private void processReadGroups(final List readGroups, boolean i ISF_LOGGER.info("chr({}) processed {} {} chimeric read groups", mChromosome, readGroupCount, scope); } - if(readGroup.reads().stream().anyMatch(x -> mConfig.Filters.skipRead(x.MateChromosome, x.MatePosStart))) + if(readGroup.reads().stream().anyMatch(x -> x.isReadPaired() && mConfig.Filters.skipRead(x.MateChromosome, x.MatePosStart))) { ++mExcludedFilteredCount; continue; diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java index eb52c6a0e0..0a0db9a61b 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionFragmentBuilder.java @@ -116,17 +116,26 @@ public static void setFragmentProperties(final FusionFragment fragment) } // set single junction info for candidate realignable fragments - if(fragment.reads().size() == 2 - && fragment.reads().stream().anyMatch(x -> hasCandidateJunctionSoftClips(x)) - && fragment.reads().stream().noneMatch(x -> x.spansGeneCollections())) + if(fragment.reads().stream().anyMatch(x -> hasCandidateJunctionSoftClips(x)) + && fragment.reads().stream().noneMatch(x -> x.spansGeneCollections())) { - FusionRead read1 = fragment.reads().get(0); - FusionRead read2 = fragment.reads().get(1); - if(read1.GeneCollections[0] == read2.GeneCollections[1]) + if(fragment.reads().size() == 1 && !fragment.reads().get(0).isReadPaired()) { setSingleSoftClipJunctionData(fragment); return; } + + if(fragment.reads().size() == 2) + { + FusionRead read1 = fragment.reads().get(0); + FusionRead read2 = fragment.reads().get(1); + + if(read1.GeneCollections[0] == read2.GeneCollections[1]) + { + setSingleSoftClipJunctionData(fragment); + return; + } + } } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadGroup.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadGroup.java index e97ea7e397..8cfc0f2474 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadGroup.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/fusion/FusionReadGroup.java @@ -53,7 +53,7 @@ public String findOtherChromosome(final String chromosome) { for(FusionRead read : mReads) { - if(!read.MateChromosome.equals(chromosome)) + if(read.isReadPaired() && !read.MateChromosome.equals(chromosome)) return read.MateChromosome; if(read.SuppData != null) diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/AltSpliceJunctionFinder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/AltSpliceJunctionFinder.java index 716dc1acb4..c9317f6830 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/AltSpliceJunctionFinder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/AltSpliceJunctionFinder.java @@ -100,12 +100,12 @@ public void setGeneData(final GeneCollection genes) } public void evaluateFragmentReads( - final List genes, final Read read1, final Read read2, final List relatedTransIds) + final List genes, final List reads, final List relatedTransIds) { if(!mEnabled) return; - if(read1.isDuplicate() || read2.isDuplicate() || genes.isEmpty() || read1.isMultiMapped() || read2.isMultiMapped()) + if(genes.isEmpty() || reads.stream().anyMatch(x -> x.isDuplicate() || x.isMultiMapped())) return; // exclude SJs too far outside known transcripts @@ -113,21 +113,31 @@ public void evaluateFragmentReads( genes.stream().mapToInt(x -> x.Gene.GeneStart).min().orElse(0) - MAX_NOVEL_SJ_DISTANCE, genes.stream().mapToInt(x -> x.Gene.GeneStart).max().orElse(0) + MAX_NOVEL_SJ_DISTANCE }; - if(!positionsWithin(read1.alignmentStart(), read1.alignmentEnd(), geneBounds[SE_START], geneBounds[SE_END]) - || !positionsWithin(read2.alignmentStart(), read2.alignmentEnd(), geneBounds[SE_START], geneBounds[SE_END])) + if(reads.stream().anyMatch(x -> !positionsWithin(x.alignmentStart(), x.alignmentEnd(), + geneBounds[SE_START], geneBounds[SE_END]))) { return; } // at least one of the reads must fall within a gene final List candidateGenes = genes.stream() - .filter(x -> positionsWithin(read1.alignmentStart(), read1.alignmentEnd(), x.Gene.GeneStart,x.Gene.GeneEnd) - || positionsWithin(read2.alignmentStart(), read2.alignmentEnd(), x.Gene.GeneStart,x.Gene.GeneEnd)) + .filter(x -> reads.stream().anyMatch(y -> positionsWithin(y.alignmentStart(), y.alignmentEnd(), x.Gene.GeneStart, x.Gene.GeneEnd))) .collect(Collectors.toList()); if(candidateGenes.isEmpty()) return; + if(reads.size() == 1) + { + if(isCandidate(reads.get(0))) + registerAltSpliceJunction(candidateGenes, reads.get(0), relatedTransIds); + + return; + } + + Read read1 = reads.get(0); + Read read2 = reads.get(1); + if(isCandidateCircular(read1, read2)) { registerAltSpliceJunction(candidateGenes, read1, read2, relatedTransIds); diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java index c3f8093fdc..ccc555296f 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/RetainedIntronFinder.java @@ -24,6 +24,7 @@ import com.google.common.collect.Lists; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.common.BaseDepth; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.GeneReadData; import com.hartwig.hmftools.isofox.common.Read; @@ -59,12 +60,12 @@ public void setGeneData(final GeneCollection genes) public final List getRetainedIntrons() { return mRetainedIntrons; } - public void evaluateFragmentReads(final Read read1, final Read read2) + public void evaluateFragmentReads(final Fragment fragment) { if(!mEnabled) return; - if(read1.isDuplicate() || read2.isDuplicate() || read1.isMultiMapped() || read2.isMultiMapped()) + if(fragment.reads().stream().anyMatch(x -> x.isDuplicate() || x.isMultiMapped())) return; // reads must span an exon boundary without being exonic in another transcript @@ -78,10 +79,8 @@ public void evaluateFragmentReads(final Read read1, final Read read2) List retIntrons = Lists.newArrayList(); - for(int i = 0; i <= 1; ++i) + for(Read read : fragment.reads()) { - Read read = (i == 0) ? read1 : read2; - if(read.containsSplit()) { splicedTrans.addAll(read.getTranscriptClassifications().entrySet().stream() @@ -112,7 +111,7 @@ public void evaluateFragmentReads(final Read read1, final Read read2) { if(retIntron1.regions().stream().anyMatch(x -> retIntron2.regions().contains(x))) { - ISF_LOGGER.trace("reads({}) support the same exon from exon-intron reads", read1.id()); + ISF_LOGGER.trace("reads({}) support the same exon from exon-intron reads", fragment.id()); return; } } diff --git a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/SpliceSiteCounter.java b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/SpliceSiteCounter.java index 100833b029..da68718ace 100644 --- a/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/SpliceSiteCounter.java +++ b/isofox/src/main/java/com/hartwig/hmftools/isofox/novel/SpliceSiteCounter.java @@ -16,7 +16,9 @@ import com.google.common.collect.Sets; import com.hartwig.hmftools.common.region.BaseRegion; import com.hartwig.hmftools.isofox.IsofoxConfig; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.GeneCollection; +import com.hartwig.hmftools.isofox.common.Read; import com.hartwig.hmftools.isofox.common.RegionReadData; public class SpliceSiteCounter @@ -38,13 +40,15 @@ public SpliceSiteCounter(final BufferedWriter writer) public void clear() { mSiteCounts.clear(); } public void registerSpliceSiteSupport( - final List readMappedCoords1, final List readMappedCoords2, final List allRegions) + final Fragment fragment, final List allRegions) { final Set traversedSites = Sets.newHashSet(); final Set supportedSites = Sets.newHashSet(); - registerSpliceSiteSupport(readMappedCoords1, allRegions, traversedSites, supportedSites); - registerSpliceSiteSupport(readMappedCoords2, allRegions, traversedSites, supportedSites); + for(Read read : fragment.reads()) + { + registerSpliceSiteSupport(read.getMappedRegionCoords(), allRegions, traversedSites, supportedSites); + } traversedSites.forEach(x -> addCount(x, SPLICE_SITE_TRAVERSED)); supportedSites.forEach(x -> addCount(x, SPLICE_SITE_SUPPORT)); diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/FragmentTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/FragmentTest.java new file mode 100644 index 0000000000..41887304fc --- /dev/null +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/FragmentTest.java @@ -0,0 +1,440 @@ +package com.hartwig.hmftools.isofox; + +import static com.hartwig.hmftools.common.bam.CigarUtils.cigarFromStr; +import static com.hartwig.hmftools.common.bam.SamRecordUtils.CONSENSUS_READ_ATTRIBUTE; +import static com.hartwig.hmftools.common.bam.SamRecordUtils.XA_ATTRIBUTE; +import static com.hartwig.hmftools.common.genome.region.Orientation.FORWARD; +import static com.hartwig.hmftools.common.genome.region.Orientation.REVERSE; +import static com.hartwig.hmftools.common.test.SamRecordTestUtils.createSamRecord; +import static com.hartwig.hmftools.isofox.TestUtils.createRegion; +import static com.hartwig.hmftools.isofox.common.RegionMatchType.EXON_BOUNDARY; +import static com.hartwig.hmftools.isofox.common.RegionMatchType.EXON_INTRON; +import static com.hartwig.hmftools.isofox.common.RegionMatchType.WITHIN_EXON; +import static com.hartwig.hmftools.isofox.common.TransMatchType.ALT; +import static com.hartwig.hmftools.isofox.common.TransMatchType.EXONIC; +import static com.hartwig.hmftools.isofox.common.TransMatchType.OTHER_TRANS; +import static com.hartwig.hmftools.isofox.common.TransMatchType.SPLICE_JUNCTION; +import static com.hartwig.hmftools.isofox.common.TransMatchType.UNKNOWN; +import static com.hartwig.hmftools.isofox.common.TransMatchType.UNSPLICED; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertFalse; +import static org.junit.Assert.assertNull; +import static org.junit.Assert.assertThrows; +import static org.junit.Assert.assertTrue; + +import java.util.List; +import java.util.Map; +import java.util.Set; + +import com.hartwig.hmftools.common.region.BaseRegion; +import com.hartwig.hmftools.isofox.common.Fragment; +import com.hartwig.hmftools.isofox.common.Read; +import com.hartwig.hmftools.isofox.common.RegionReadData; + +import org.junit.Test; + +import htsjdk.samtools.SAMRecord; + +public class FragmentTest +{ + @Test + public void testSingleEndFragment() + { + Read read = new Read(createRecord(100, "5S10M100N10M5S")); + Fragment fragment = new Fragment(read); + + assertFalse(read.isReadPaired()); + assertEquals("fragment", fragment.id()); + assertEquals("1", fragment.chromosome()); + assertEquals(List.of(read), fragment.reads()); + assertThrows(UnsupportedOperationException.class, () -> fragment.reads().clear()); + assertEquals(1, fragment.fragmentCount()); + assertEquals(100, fragment.minAlignmentStart()); + assertEquals(219, fragment.maxAlignmentEnd()); + assertTrue(fragment.containsSplit()); + assertTrue(fragment.isFullyIntronic()); + assertTrue(fragment.uniqueValidRegions().isEmpty()); + assertFalse(fragment.spansMultipleRegions(1)); + assertFalse(fragment.readsInDifferentExons(1)); + assertEquals(List.of(new BaseRegion(100, 109), new BaseRegion(210, 219)), fragment.mergedMappings()); + } + + @Test + public void testPairedFragmentInEitherOrder() + { + Read first = new Read(createRecord(200, "20M")); + Read second = new Read(createRecord(100, "10S150M5S")); + markAsPair(first, second); + + assertEquals(List.of(first, second), new Fragment(first, second).reads()); + assertEquals(List.of(second, first), new Fragment(second, first).reads()); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertEquals("fragment", fragment.id()); + assertEquals("1", fragment.chromosome()); + assertThrows(UnsupportedOperationException.class, () -> fragment.reads().add(first)); + assertEquals(1, fragment.fragmentCount()); + assertEquals(100, fragment.minAlignmentStart()); + assertEquals(249, fragment.maxAlignmentEnd()); + assertFalse(fragment.containsSplit()); + assertTrue(fragment.isFullyIntronic()); + assertEquals(List.of(new BaseRegion(100, 249)), fragment.mergedMappings()); + } + } + + @Test + public void testConsensusCountIsNotSummedAcrossMates() + { + SAMRecord firstRecord = createRecord(100, "20M"); + SAMRecord secondRecord = createRecord(200, "20M"); + firstRecord.setAttribute(CONSENSUS_READ_ATTRIBUTE, "5;5"); + secondRecord.setAttribute(CONSENSUS_READ_ATTRIBUTE, "5;5"); + Read first = new Read(firstRecord); + Read second = new Read(secondRecord); + + assertEquals(5, new Fragment(first).fragmentCount()); + + markAsPair(first, second); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertEquals(5, fragment.fragmentCount()); + } + } + + @Test + public void testLeastAmbiguousReadDeterminesLoci() + { + Read unique = new Read(createRecord(100, "20M")); + Read ambiguous = createMultiMappedRead("2,+500,20M,0;3,-700,20M,1;"); + Read lessAmbiguous = createMultiMappedRead("4,+900,20M,0;"); + + assertEquals(1, new Fragment(unique).minNumLoci()); + assertNull(new Fragment(unique).altLoci()); + assertEquals(3, new Fragment(ambiguous).minNumLoci()); + assertEquals(2, new Fragment(ambiguous).altLoci().size()); + + markAsPair(ambiguous, unique); + + for(Fragment fragment : List.of(new Fragment(ambiguous, unique), new Fragment(unique, ambiguous))) + { + assertEquals(1, fragment.minNumLoci()); + assertNull(fragment.altLoci()); + } + + markAsPair(ambiguous, lessAmbiguous); + + for(Fragment fragment : List.of(new Fragment(ambiguous, lessAmbiguous), new Fragment(lessAmbiguous, ambiguous))) + { + assertEquals(2, fragment.minNumLoci()); + assertEquals(1, fragment.altLoci().size()); + assertEquals(900, fragment.altLoci().get(0).Region.start()); + } + } + + @Test + public void testEqualLocusCountsKeepFirstSuppliedRead() + { + Read first = createMultiMappedRead("2,+500,20M,0;"); + Read second = createMultiMappedRead("3,+700,20M,0;"); + markAsPair(first, second); + + Fragment firstSupplied = new Fragment(first, second); + Fragment secondSupplied = new Fragment(second, first); + + assertEquals(2, firstSupplied.minNumLoci()); + assertEquals(500, firstSupplied.altLoci().get(0).Region.start()); + assertEquals(2, secondSupplied.minNumLoci()); + assertEquals(700, secondSupplied.altLoci().get(0).Region.start()); + } + + @Test + public void testMergedMappings() + { + Read split = new Read(createRecord(100, "10M100N10M")); + Read overlapping = new Read(createRecord(105, "10M100N15M")); + + assertEquals(List.of(new BaseRegion(100, 109), new BaseRegion(210, 219)), new Fragment(split).mergedMappings()); + assertEquals(List.of(new BaseRegion(105, 114), new BaseRegion(215, 229)), new Fragment(overlapping).mergedMappings()); + + markAsPair(split, overlapping); + + for(Fragment fragment : List.of(new Fragment(split, overlapping), new Fragment(overlapping, split))) + { + assertEquals(List.of(new BaseRegion(100, 114), new BaseRegion(210, 229)), fragment.mergedMappings()); + } + + Read splitRead = new Read(createRecord(100, "10M100N10M")); + Read downstream = new Read(createRecord(300, "20M")); + markAsPair(splitRead, downstream); + + for(Fragment fragment : List.of(new Fragment(splitRead, downstream), new Fragment(downstream, splitRead))) + { + assertTrue(fragment.containsSplit()); + assertEquals(List.of(new BaseRegion(100, 109), new BaseRegion(210, 219), new BaseRegion(300, 319)), fragment.mergedMappings()); + } + + Read upstream = new Read(createRecord(100, "30M")); + Read adjacent = new Read(createRecord(120, "30M")); + markAsPair(upstream, adjacent); + + for(Fragment fragment : List.of(new Fragment(upstream, adjacent), new Fragment(adjacent, upstream))) + { + assertEquals(List.of(new BaseRegion(100, 149)), fragment.mergedMappings()); + } + } + + @Test + public void testUniqueValidRegionsDeduplicatesAndSkipsExonIntron() + { + Read first = new Read(createRecord(100, "10M100N10M")); + Read second = new Read(createRecord(100, "20M")); + RegionReadData shared = createRegion("gene", 1, 1, "1", 100, 199); + RegionReadData exonIntron = createRegion("gene", 1, 2, "1", 210, 299); + RegionReadData unrelated = createRegion("other", 2, 1, "1", 300, 399); + first.getMappedRegions().putAll(Map.of(shared, WITHIN_EXON, exonIntron, EXON_INTRON, unrelated, EXON_BOUNDARY)); + second.getMappedRegions().put(shared, EXON_BOUNDARY); + + assertFalse(new Fragment(first).isFullyIntronic()); + assertEquals(Set.of(shared, unrelated), Set.copyOf(new Fragment(first).uniqueValidRegions())); + + markAsPair(first, second); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertFalse(fragment.isFullyIntronic()); + assertEquals(2, fragment.uniqueValidRegions().size()); + assertEquals(Set.of(shared, unrelated), Set.copyOf(fragment.uniqueValidRegions())); + } + } + + @Test + public void testSpansMultipleRegionsCombinesBothReads() + { + RegionReadData exon1 = createRegion("gene", 1, 1, "1", 100, 199); + RegionReadData exon2 = createRegion("gene", 1, 2, "1", 210, 299); + + Read exonIntronRead = new Read(createRecord(100, "10M100N10M")); + exonIntronRead.getMappedRegions().putAll(Map.of(exon1, WITHIN_EXON, exon2, EXON_INTRON)); + assertFalse(new Fragment(exonIntronRead).spansMultipleRegions(1)); + + Read bothExonsRead = new Read(createRecord(100, "10M100N10M")); + bothExonsRead.getMappedRegions().putAll(Map.of(exon1, WITHIN_EXON, exon2, EXON_BOUNDARY)); + assertTrue(new Fragment(bothExonsRead).spansMultipleRegions(1)); + assertFalse(new Fragment(bothExonsRead).spansMultipleRegions(99)); + + Read sameExonFirst = new Read(createRecord(100, "20M")); + Read sameExonSecond = new Read(createRecord(120, "20M")); + sameExonFirst.getMappedRegions().put(exon1, WITHIN_EXON); + sameExonSecond.getMappedRegions().put(exon1, WITHIN_EXON); + markAsPair(sameExonFirst, sameExonSecond); + + for(Fragment fragment : List.of(new Fragment(sameExonFirst, sameExonSecond), new Fragment(sameExonSecond, sameExonFirst))) + { + assertFalse(fragment.spansMultipleRegions(1)); + } + + Read exon1Read = new Read(createRecord(100, "20M")); + Read exon2Read = new Read(createRecord(210, "20M")); + exon1Read.getMappedRegions().put(exon1, WITHIN_EXON); + exon2Read.getMappedRegions().put(exon2, EXON_BOUNDARY); + assertFalse(new Fragment(exon1Read).spansMultipleRegions(1)); + + markAsPair(exon1Read, exon2Read); + + for(Fragment fragment : List.of(new Fragment(exon1Read, exon2Read), new Fragment(exon2Read, exon1Read))) + { + assertTrue(fragment.spansMultipleRegions(1)); + } + } + + @Test + public void testPairedTranscriptSupportIsAnIntersection() + { + Read first = new Read(createRecord(100, "20M")); + Read second = new Read(createRecord(200, "20M")); + first.getTranscriptClassifications().putAll(Map.of(1, EXONIC, 2, SPLICE_JUNCTION, 3, ALT, 4, EXONIC, 5, EXONIC)); + second.getTranscriptClassifications().putAll(Map.of(1, SPLICE_JUNCTION, 2, EXONIC, 3, EXONIC, 5, ALT, 6, SPLICE_JUNCTION)); + + assertEquals(Set.of(1, 2, 4, 5), Set.copyOf(new Fragment(first).validTypeTranscripts())); + assertEquals(Set.of(3), new Fragment(first).invalidTranscripts(List.of(1, 2, 4, 5))); + + markAsPair(first, second); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertEquals(2, fragment.validTypeTranscripts().size()); + assertEquals(Set.of(1, 2), Set.copyOf(fragment.validTypeTranscripts())); + assertEquals(Set.of(3, 4, 5, 6), fragment.invalidTranscripts(List.of(1, 2))); + assertEquals(Set.of(2, 3, 4, 5, 6), fragment.invalidTranscripts(List.of(1))); + assertTrue(fragment.hasTranscriptClassification(1, EXONIC)); + assertTrue(fragment.hasTranscriptClassification(1, SPLICE_JUNCTION)); + assertTrue(fragment.hasTranscriptClassification(6, SPLICE_JUNCTION)); + assertFalse(fragment.hasTranscriptClassification(1, ALT)); + assertFalse(fragment.hasTranscriptClassification(99, EXONIC)); + } + + second.getTranscriptClassifications().clear(); + assertTrue(new Fragment(second).validTypeTranscripts().isEmpty()); + assertTrue(new Fragment(second).invalidTranscripts(List.of()).isEmpty()); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertTrue(fragment.validTypeTranscripts().isEmpty()); + assertEquals(Set.of(1, 2, 3, 4, 5), fragment.invalidTranscripts(List.of())); + } + } + + @Test + public void testSetOtherTranscriptsPreservesAcceptedAndInvalidTypes() + { + for(boolean paired : List.of(false, true)) + { + Read first = new Read(createRecord(100, "20M")); + Read second = new Read(createRecord(200, "20M")); + + if(paired) + markAsPair(first, second); + + Fragment fragment = paired ? new Fragment(first, second) : new Fragment(first); + + for(Read read : fragment.reads()) + { + read.getTranscriptClassifications().putAll(Map.of( + 1, EXONIC, 2, SPLICE_JUNCTION, 3, ALT, 4, UNSPLICED, 5, UNKNOWN, 6, OTHER_TRANS, 7, EXONIC)); + } + + fragment.setOtherTranscripts(List.of(1)); + + for(Read read : fragment.reads()) + { + assertEquals(Map.of(1, EXONIC, 2, OTHER_TRANS, 3, ALT, 4, UNSPLICED, 5, UNKNOWN, 6, OTHER_TRANS, 7, OTHER_TRANS), + read.getTranscriptClassifications()); + } + + assertEquals(List.of(1), fragment.validTypeTranscripts()); + } + } + + @Test + public void testDifferentExonsCompareRanksWithinTheRequestedTranscript() + { + Read first = new Read(createRecord(100, "20M")); + Read second = new Read(createRecord(300, "20M")); + first.getMappedRegions().put(createRegion("gene", 1, 1, "1", 100, 199), WITHIN_EXON); + first.getMappedRegions().put(createRegion("other", 2, 2, "1", 200, 299), WITHIN_EXON); + second.getMappedRegions().put(createRegion("gene", 1, 2, "1", 300, 399), WITHIN_EXON); + + assertFalse(new Fragment(first).readsInDifferentExons(1)); + + markAsPair(first, second); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertTrue(fragment.readsInDifferentExons(1)); + } + + second.getMappedRegions().put(createRegion("gene", 1, 1, "1", 100, 199), WITHIN_EXON); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertFalse(fragment.readsInDifferentExons(1)); + } + } + + @Test + public void testOrientationUsesFirstOfPairRatherThanPositionOrListOrder() + { + for(boolean reversed : List.of(false, true)) + { + Read first = new Read(createRecord(200, "20M")); + Read second = new Read(createRecord(100, "20M")); + first.bamRecord().setReadNegativeStrandFlag(reversed); + second.bamRecord().setReadNegativeStrandFlag(!reversed); + + assertEquals(reversed ? REVERSE : FORWARD, new Fragment(first).orientation()); + + markAsPair(first, second); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertEquals(reversed ? REVERSE : FORWARD, fragment.orientation()); + } + + second.bamRecord().setReadNegativeStrandFlag(reversed); + + for(Fragment fragment : List.of(new Fragment(first, second), new Fragment(second, first))) + { + assertNull(fragment.orientation()); + } + } + } + + @Test + public void testAdapterTrimmingRequiresAMate() + { + for(int order : List.of(0, 1)) + { + Read first = new Read(createRecord(105, "10S20M10S")); + Read second = new Read(createRecord(100, "10S20M10S")); + second.bamRecord().setReadNegativeStrandFlag(true); + String firstBases = first.readBases(); + String secondBases = second.readBases(); + + for(Read read : List.of(first, second)) + { + new Fragment(read).trimAdapterBases(); + assertEquals("10S20M10S", read.cigarStr()); + assertEquals(40, read.baseLength()); + } + + markAsPair(first, second); + + Fragment fragment = order == 0 ? new Fragment(first, second) : new Fragment(second, first); + fragment.trimAdapterBases(); + + assertEquals("10S20M5S", first.cigarStr()); + assertEquals("5S20M10S", second.cigarStr()); + assertEquals(firstBases.substring(0, 35), first.readBases()); + assertEquals(secondBases.substring(5), second.readBases()); + assertEquals(129, first.unclippedEnd()); + assertEquals(95, second.unclippedStart()); + assertEquals(List.of(new BaseRegion(100, 124)), fragment.mergedMappings()); + } + } + + private static SAMRecord createRecord(int start, final String cigar) + { + int readLength = cigarFromStr(cigar).getReadLength(); + String bases = "ACGT".repeat((readLength + 3) / 4).substring(0, readLength); + SAMRecord record = createSamRecord("fragment", "1", start, bases, cigar, "*", 0, false, false, null); + record.setFlags(0); + record.setMateAlignmentStart(0); + record.setInferredInsertSize(0); + return record; + } + + private static Read createMultiMappedRead(final String xaTag) + { + SAMRecord record = createRecord(100, "20M"); + record.setAttribute(XA_ATTRIBUTE, xaTag); + return new Read(record); + } + + private static void markAsPair(final Read first, final Read second) + { + for(Read read : List.of(first, second)) + { + Read mate = read == first ? second : first; + SAMRecord record = read.bamRecord(); + record.setReadPairedFlag(true); + record.setFirstOfPairFlag(read == first); + record.setSecondOfPairFlag(read == second); + record.setMateReferenceName(mate.chromosome()); + record.setMateAlignmentStart(mate.alignmentStart()); + record.setMateNegativeStrandFlag(mate.isReadReversed()); + } + } +} diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/NovelJunctionsTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/NovelJunctionsTest.java index b0800fd730..7cc86a5ea4 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/NovelJunctionsTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/NovelJunctionsTest.java @@ -36,6 +36,7 @@ import com.hartwig.hmftools.common.gene.TranscriptData; import com.hartwig.hmftools.common.test.MockRefGenome; import com.hartwig.hmftools.isofox.adjusts.FragmentSize; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.GeneReadData; import com.hartwig.hmftools.isofox.common.Read; @@ -329,7 +330,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(0, riFinder.getRetainedIntrons().size()); @@ -340,7 +341,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(1, riFinder.getRetainedIntrons().size()); RetainedIntron retIntron = riFinder.getRetainedIntrons().get(0); @@ -355,7 +356,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(2, riFinder.getRetainedIntrons().size()); retIntron = riFinder.getRetainedIntrons().get(1); @@ -372,7 +373,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(2, riFinder.getRetainedIntrons().size()); retIntron = riFinder.getRetainedIntrons().get(1); @@ -391,7 +392,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(2, riFinder.getRetainedIntrons().size()); @@ -401,7 +402,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(2, riFinder.getRetainedIntrons().size()); @@ -414,7 +415,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertTrue(riFinder.getRetainedIntrons().isEmpty()); @@ -426,7 +427,7 @@ public void testRetainedIntrons() processOverlappingRegions(read1, gene.findOverlappingRegions(read1)); processOverlappingRegions(read2, gene.findOverlappingRegions(read2)); - riFinder.evaluateFragmentReads(read1, read2); + riFinder.evaluateFragmentReads(new Fragment(read1, read2)); assertEquals(2, riFinder.getRetainedIntrons().size()); diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/TransClassificationTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/TransClassificationTest.java index 329d72783f..2b46c44ad5 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/TransClassificationTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/TransClassificationTest.java @@ -30,6 +30,7 @@ import com.google.common.collect.Lists; import com.hartwig.hmftools.common.gene.ExonData; import com.hartwig.hmftools.common.gene.TranscriptData; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.FragmentType; import com.hartwig.hmftools.isofox.common.FragmentTypeCounts; import com.hartwig.hmftools.isofox.common.GeneCollection; @@ -214,21 +215,24 @@ public void testFragmentLengthCalcs() Read read1 = createReadRecord(1, CHR_1, 1010, 1029, REF_BASE_STR_1, createCigar(0, 20, 0)); Read read2 = createReadRecord(1, CHR_1, 1170, 1199, REF_BASE_STR_1, createCigar(0, 30, 0)); - int fragLength = calcFragmentLength(transData, read1, read2); + Fragment fragment = new Fragment(read1, read2); + int fragLength = calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); assertEquals(190, fragLength); // spanning 2 exons, both exonic read1 = createReadRecord(1, CHR_1, 1170, 1189, REF_BASE_STR_1, createCigar(0, 20, 0)); read2 = createReadRecord(1, CHR_1, 2010, 2019, REF_BASE_STR_1, createCigar(0, 10, 0)); - fragLength = calcFragmentLength(transData, read1, read2); + fragment = new Fragment(read1, read2); + fragLength = calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); assertEquals(31 + 20, fragLength); // spanning 3 exons, both exonic read1 = createReadRecord(1, CHR_1, 1170, 1189, REF_BASE_STR_1, createCigar(0, 20, 0)); read2 = createReadRecord(1, CHR_1, 4510, 4519, REF_BASE_STR_1, createCigar(0, 10, 0)); - fragLength = calcFragmentLength(transData, read1, read2); + fragment = new Fragment(read1, read2); + fragLength = calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); assertEquals(31 + 501 + 20, fragLength); // with 2 split reads @@ -237,7 +241,8 @@ public void testFragmentLengthCalcs() read2 = createReadRecord( 1, CHR_1, 2491, 4509, REF_BASE_STR_1, createCigar(0, 10, 1999, 10, 0)); - fragLength = calcFragmentLength(transData, read1, read2); + fragment = new Fragment(read1, read2); + fragLength = calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); assertEquals(10 + 501 + 10, fragLength); // with 2 split reads skipping an exon @@ -250,7 +255,8 @@ public void testFragmentLengthCalcs() read2 = createReadRecord( 1, CHR_1, 4991, 5509, REF_BASE_STR_1, createCigar(0, 10, 499, 10, 0)); - fragLength = calcFragmentLength(transData, read1, read2); + fragment = new Fragment(read1, read2); + fragLength = calcFragmentLength(transData, fragment.minAlignmentStart(), fragment.maxAlignmentEnd()); assertEquals(10 + 501 + 501 + 10, fragLength); } diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java index fbe115f234..b0e0e1f885 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/ChimericReadTest.java @@ -39,6 +39,7 @@ import com.hartwig.hmftools.common.gene.GeneData; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.common.BaseDepth; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.FragmentTracker; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.Read; @@ -78,7 +79,7 @@ public void testBasicReads() Read read2 = createMappedRead(readId, gc1, 1081, 1100, createCigar(0, 20, 20)); read2.setSuppAlignment(TEST_SUPP_DATA); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertEquals(1, chimericRT.fusionReadGroupMap().size()); @@ -98,7 +99,7 @@ public void testBasicReads() read1.setFlag(FIRST_OF_PAIR, true); read2 = createMappedRead(readId, gc1, 1066, 1100, createCigar(0, 35, 5)); - chimericRT.addRealignmentCandidates(read1, read2); + chimericRT.addRealignmentCandidates(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); @@ -143,7 +144,7 @@ public void testSameGeneCollection() Read read2 = createMappedRead(readId, gc1, 1500, 1519, createCigar(20, 20, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertTrue(chimericRT.fusionReadGroupMap().isEmpty()); @@ -160,7 +161,7 @@ public void testSameGeneCollection() //read2 = createMappedRead(readId, gc1, 10400, 10419, createCigar(20, 20, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertEquals(1, chimericRT.fusionReadGroupMap().size()); @@ -177,7 +178,7 @@ public void testSameGeneCollection() read2 = createMappedRead(readId, gc1, 10450, 10469, createCigar(20, 20, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertTrue(chimericRT.fusionReadGroupMap().isEmpty()); @@ -199,7 +200,7 @@ public void testSameGeneCollection() chimericRT.initialise(gc2); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertTrue(chimericRT.fusionReadGroupMap().isEmpty()); @@ -223,7 +224,7 @@ public void testSameGeneCollection() read2.setGeneCollection(SE_START, gc2.id(), false); read2.setGeneCollection(SE_END, gc2.id(), false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertEquals(1, chimericRT.fusionReadGroupMap().size()); @@ -263,7 +264,7 @@ public void testPrePosGeneReads() Read read2 = createMappedRead(readId, gc1, 2000, 2019, createCigar(20, 20, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); assertTrue(chimericRT.fusionReadGroupMap().isEmpty()); @@ -306,8 +307,8 @@ public void testPrePosGeneReads() fragTracker.checkRead(read1); fragTracker.checkRead(read4); - chimericRT.addChimericReadPair(read2, read3); - chimericRT.addChimericReadPair(read5, read6); + chimericRT.addChimericFragment(new Fragment(read2, read3)); + chimericRT.addChimericFragment(new Fragment(read5, read6)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); @@ -328,7 +329,7 @@ public void testPrePosGeneReads() chimericRT.initialise(gc2); baseDepth.initialise(gc2.regionBounds()); - chimericRT.addChimericReadPair(read2, read3); + chimericRT.addChimericFragment(new Fragment(read2, read3)); fragTracker.checkRead(read4); fragTracker.checkRead(read6); @@ -373,14 +374,14 @@ public void testJunctionPositionTracking() Read read2 = createMappedRead(readId, gc1, 1050, 1089, createCigar(0, 40, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); read1 = createMappedRead(++readId, gc1, 1081, 1100, createCigar(0, 20, 3)); read1.setFlag(FIRST_OF_PAIR, true); read2 = createMappedRead(readId, gc1, 1050, 1089, createCigar(0, 40, 0)); read2.setStrand(true, false); - chimericRT.addRealignmentCandidates(read1, read2); + chimericRT.addRealignmentCandidates(new Fragment(read1, read2)); // another one split to the 3rd gene read1 = createMappedRead(++readId, gc1, 1281, 20419, createCigar(0, 20, 19100, 20, 0)); @@ -389,7 +390,7 @@ public void testJunctionPositionTracking() read2 = createMappedRead(readId, gc1, 1050, 1089, createCigar(0, 40, 0)); read2.setStrand(true, false); - chimericRT.addChimericReadPair(read1, read2); + chimericRT.addChimericFragment(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); @@ -406,7 +407,7 @@ public void testJunctionPositionTracking() read2 = createMappedRead(readId, gc1, 10210, 10249, createCigar(0, 40, 0)); read2.setStrand(true, false); - chimericRT.addRealignmentCandidates(read1, read2); + chimericRT.addRealignmentCandidates(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); @@ -423,7 +424,7 @@ public void testJunctionPositionTracking() read2 = createMappedRead(readId, gc1, 20410, 20449, createCigar(0, 40, 0)); read2.setStrand(true, false); - chimericRT.addRealignmentCandidates(read1, read2); + chimericRT.addRealignmentCandidates(new Fragment(read1, read2)); chimericRT.postProcessChimericReads(baseDepth, fragTracker); diff --git a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionDataTest.java b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionDataTest.java index 8873df53a7..837daebdbb 100644 --- a/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionDataTest.java +++ b/isofox/src/test/java/com/hartwig/hmftools/isofox/fusion/FusionDataTest.java @@ -51,6 +51,7 @@ import com.hartwig.hmftools.common.fusion.KnownFusionData; import com.hartwig.hmftools.isofox.IsofoxConfig; import com.hartwig.hmftools.isofox.common.BaseDepth; +import com.hartwig.hmftools.isofox.common.Fragment; import com.hartwig.hmftools.isofox.common.GeneCollection; import com.hartwig.hmftools.isofox.common.Read; @@ -289,7 +290,7 @@ public void testInterChromosomalFusion() Read read2 = createMappedRead(readId, gc3, 20250, 20289, createCigar(0, 40, 0)); read2.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read1, read2), ORIENT_FWD, 20300); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read1, read2)), ORIENT_FWD, 20300); // RAC fragment for GC5 junctionBases = config.RefGenome.getBaseString(gc5.chromosome(), 20298, 20300) @@ -299,7 +300,7 @@ public void testInterChromosomalFusion() read2 = createMappedRead(readId, gc5, 10210, 10259, createCigar(0, 40, 0)); read2.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read1, read2), ORIENT_REV, 10200); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read1, read2)), ORIENT_REV, 10200); // 1 intronic discordant read Read[] discordantReads = createReadPair(++readId, gc3, gc5, 20150, 20189, 10320, 10359, @@ -401,7 +402,7 @@ public void testSoftClippedFragmentRealignment() Read read5 = createMappedRead(readId, gc1, 1051, 1090, createCigar(0, 40, 0)); read5.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read4, read4), ORIENT_FWD, 1100); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read4, read4)), ORIENT_FWD, 1100); // a soft-clipped read matching 2 bases into the ref due to homology with the other side of the fusion junction junctionBases = config.RefGenome.getBaseString(gc1.chromosome(), 1091, 1100) @@ -411,7 +412,7 @@ public void testSoftClippedFragmentRealignment() Read read7 = createMappedRead(readId, gc2, 10210, 10249, createCigar(0, 40, 0)); read7.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read6, read7), ORIENT_REV, 10200); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read6, read7)), ORIENT_REV, 10200); BaseDepth baseDepth = new BaseDepth(); List completeGroups = finder.processNewChimericReadGroups(gc1, baseDepth, readGroups1); @@ -484,7 +485,7 @@ public void testLocalDelFusion() // readGroups1.put(read3.id(), new ReadGroup(read3, read4)); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read3, read4), ORIENT_FWD, 1100); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read3, read4)), ORIENT_FWD, 1100); junctionBases = config.RefGenome.getBaseString(gc1.chromosome(), 1091, 1100) + config.RefGenome.getBaseString(gc2.chromosome(), 10200, 10229); @@ -495,7 +496,7 @@ public void testLocalDelFusion() readGroups1.put(read5.id(), createGroup(read5, read6)); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read5, read6), ORIENT_REV, 10200); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read5, read6)), ORIENT_REV, 10200); // and a discordant fragment read3 = createMappedRead(++readId, gc1, 1050, 1089, createCigar(0, 40, 0)); @@ -566,7 +567,7 @@ public void testLocalSplitReadDelFusions() Read read4 = createMappedRead(readId, gc1, 1051, 1090, createCigar(0, 40, 0)); read4.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read3, read4), ORIENT_FWD, 1100); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read3, read4)), ORIENT_FWD, 1100); // then on GC2 junctionBases = config.RefGenome.getBaseString(gc1.chromosome(), 1091, 1100) @@ -575,7 +576,7 @@ public void testLocalSplitReadDelFusions() Read read6 = createMappedRead(readId, gc2, 10210, 10249, createCigar(0, 40, 0)); read6.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read5, read6), ORIENT_REV, 10200); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read5, read6)), ORIENT_REV, 10200); // and a discordant fragment Read read7 = createMappedRead(++readId, gc1, 1050, 1089, createCigar(0, 40, 0)); @@ -665,7 +666,7 @@ public void testNonGenicFusions() Read read4 = createMappedRead(readId, gc1, 551, 590, createCigar(0, 40, 0)); read4.setStrand(true, false); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read3, read4), ORIENT_FWD, 600); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read3, read4)), ORIENT_FWD, 600); junctionBases = config.RefGenome.getBaseString(gc1.chromosome(), 591, 600) + config.RefGenome.getBaseString(gc2.chromosome(), 10200, 10229); @@ -675,7 +676,7 @@ public void testNonGenicFusions() read6.setStrand(true, false); readGroups1.put(read5.id(), createGroup(read5 ,read6)); - addRacReadGroup(racFragmentCache, new ChimericReadGroup(read5, read6), ORIENT_REV, 10200); + addRacReadGroup(racFragmentCache, new ChimericReadGroup(new Fragment(read5, read6)), ORIENT_REV, 10200); // and a discordant fragment read3 = createMappedRead(++readId, gc1, 550, 589, createCigar(0, 40, 0));