Add support for RNA CRAM input

[oneillkza]

Description of feature

I'm just trying to run oncoanalyser on some of our WGTA data, and all of our bams are archived to cram as a matter of policy. However, it seems oncoanalyser doesn't like the RNA-seq to be input as CRAM:

ERROR ~ found tumor RNA CRAM input for POG1342_WGTS but RNA CRAM input is not supported

I'm in the process of uncramming to try again, but it would be helpful not to need that extra step.

Would it be possible to please request support of CRAM input for RNA?

[luan-n-nguyen]

Hi there, this restriction will be removed in the upcoming version 3.0.0

[oneillkza]

Thanks @luan-n-nguyen , that's good to know!

Do you have a rough timeline of when the 3.0.0 release is expected?

[scwatts]

@luan-n-nguyen just for awareness, the specific reason RNA CRAM inputs were disallowed despite working with ISOFOX is that SAGE append encountered an error with some data:

07:10:10.432 [ERROR] task execution error: java.util.concurrent.ExecutionException: htsjdk.samtools.SAMException: Exception counting mismatches for read ENST00000598249.6|ENSG00000160460.16|OTTHUMG00000182592.6|OTTHUMT00000463877.2|SPTBN4-212|SPTBN4|8910|UTR5:1-320|CDS:321-8015|UTR3:8016-8910|_3638_4187_0:0:0_0:0:0_f7 2/2 250b aligned to chr19:40520115-40529050.
java.util.concurrent.ExecutionException: htsjdk.samtools.SAMException: Exception counting mismatches for read ENST00000598249.6|ENSG00000160460.16|OTTHUMG00000182592.6|OTTHUMT00000463877.2|SPTBN4-212|SPTBN4|8910|UTR5:1-320|CDS:321-8015|UTR3:8016-8910|_3638_4187_0:0:0_0:0:0_f7 2/2 250b aligned to chr19:40520115-40529050.
        at java.base/java.util.concurrent.FutureTask.report(FutureTask.java:122)
        at java.base/java.util.concurrent.FutureTask.get(FutureTask.java:191)
        at com.hartwig.hmftools.common.perf.TaskExecutor.checkTaskCompletion(TaskExecutor.java:100)
        at com.hartwig.hmftools.common.perf.TaskExecutor.executeTasks(TaskExecutor.java:77)
        at com.hartwig.hmftools.sage.append.SageAppendApplication.run(SageAppendApplication.java:226)
        at com.hartwig.hmftools.sage.append.SageAppendApplication.main(SageAppendApplication.java:323)
Caused by: htsjdk.samtools.SAMException: Exception counting mismatches for read ENST00000598249.6|ENSG00000160460.16|OTTHUMG00000182592.6|OTTHUMT00000463877.2|SPTBN4-212|SPTBN4|8910|UTR5:1-320|CDS:321-8015|UTR3:8016-8910|_3638_4187_0:0:0_0:0:0_f7 2/2 250b aligned to chr19:40520115-40529050.
        at htsjdk.samtools.util.SequenceUtil.countMismatches(SequenceUtil.java:490)
        at htsjdk.samtools.util.SequenceUtil.calculateSamNmTag(SequenceUtil.java:635)
        at htsjdk.samtools.util.SequenceUtil.calculateSamNmTag(SequenceUtil.java:619)
        at com.hartwig.hmftools.sage.common.NumberEvents.rawNM(NumberEvents.java:55)
        at com.hartwig.hmftools.sage.common.NumberEvents.calc(NumberEvents.java:18)
        at com.hartwig.hmftools.sage.evidence.ReadContextEvidence.processReadRecord(ReadContextEvidence.java:299)
        at com.hartwig.hmftools.sage.evidence.ReadContextEvidence.processReadRecord(ReadContextEvidence.java:227)
        at com.hartwig.hmftools.sage.evidence.ReadContextEvidence.processReadRecord(ReadContextEvidence.java:223)
        at com.hartwig.hmftools.common.bam.BamSlicer.slice(BamSlicer.java:86)
        at com.hartwig.hmftools.sage.common.SamSlicer.slice(SamSlicer.java:28)
        at com.hartwig.hmftools.sage.evidence.ReadContextEvidence.collectEvidence(ReadContextEvidence.java:135)
        at com.hartwig.hmftools.sage.pipeline.EvidenceStage.findEvidence(EvidenceStage.java:59)
        at com.hartwig.hmftools.sage.append.RegionAppendTask.call(RegionAppendTask.java:108)
        at com.hartwig.hmftools.sage.append.RegionAppendTask.call(RegionAppendTask.java:38)
        at java.base/java.util.concurrent.FutureTask.run(FutureTask.java:317)
        at java.base/java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1144)
        at java.base/java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:642)
        at java.base/java.lang.Thread.run(Thread.java:1583)
Caused by: java.lang.ArrayIndexOutOfBoundsException: Index -6151 out of bounds for length 2061
        at htsjdk.samtools.util.SequenceUtil.countMismatches(SequenceUtil.java:482)
        ... 17 more

[luan-n-nguyen]

Thanks, will do a test run with RNA CRAM anyway to confirm OA doesn't break

[luan-n-nguyen]

Thanks @luan-n-nguyen , that's good to know!

Do you have a rough timeline of when the 3.0.0 release is expected?

We still have quite some testing to do, but i estimate around 2 months

[scwatts]

Appreciate you taking a look at this @luan-n-nguyen!

[luan-n-nguyen]

@scwatts Seems like the Sage append 'Exception counting mismatches' error happens when i use the below samtools command from the OA readme to create the CRAM:

samtools view \
--cram \
--output-fmt-option version=3.0 \
--output-fmt-option embed_ref=1 \
--output-fmt-option reference=/opt/resources/oncoanalyser/genome/38/GRCh38_masked_exclusions_alts_hlas.fasta \
--output /data/samples/colo_mini/cram/38/COLO829_RNA.cram \
--threads 4 \
--write-index \
/data/samples/colo_mini/bam/38/COLO829_RNA.bam

However, adding this samtools arg fixes the issue (i'll add this to the OA docs):

--output-fmt-option store_nm=1

This prevents Sage from calling htsjdk to recompute the NM tag at NumberEvents.java#L60. I'm guessing the NM tag recompute doesn't work properly with RNA data due to splicing (STAR probably has some special way to compute the NM tag).

Tested/confirmed working with the following OA run setup:

WORKING_DIR=/data/experiments/oncoanalyser/20260504_oa_dev_3.0.0-beta/; cd $WORKING_DIR

echo "group_id,subject_id,sample_id,sample_type,sequence_type,filetype,filepath
colo_mini,colo_mini,COLO829T,tumor,dna,bam,/data/samples/colo_mini/bam/38/COLO829T.bam
colo_mini,colo_mini,COLO829R,normal,dna,bam,/data/samples/colo_mini/bam/38/COLO829R.bam
colo_mini,colo_mini,COLO829T-RNA,tumor,rna,cram,/data/samples/colo_mini/cram/38/COLO829_RNA.cram
" > ${WORKING_DIR}/sample_sheets/sample_sheet.colo_mini.tumor_normal_rna.bam.38.csv

OUT_DIR=${WORKING_DIR}/output/colo_mini.tumor_normal_rna/
nextflow run /opt/tools/oncoanalyser \
-profile docker \
-config ${WORKING_DIR}/config/compute_resources.config \
-config /opt/resources/oncoanalyser/genome/38/genome.config \
--genome GRCh38_hmf \
--ref_data_hmf_data_path /opt/resources/oncoanalyser/ref_data/3.0.0--4/38/hmftools/ \
--input ${WORKING_DIR}/sample_sheets/sample_sheet.colo_mini.tumor_normal_rna.bam.38.csv \
--mode wgts \
--outdir ${OUT_DIR} \
-work-dir ${OUT_DIR}/work/ \
-resume