Difference of running on Local and on Terra

[shengqh]

I wrote a WDL for STAR-Fusion: https://github.com/shengqh/warp/blob/develop/pipelines/vumc_biostatistics/rnaseq/VUMCStarFusion.wdl. I tested the workflow on both local cluster and Terra using same raw FASTQ file, same reference GRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play.tar.gz and same parameters. The weird thing is, in the STAR step, the local cluster can generate correct Chimeric.out.junction while Terra one generate a file with header and comment only.

Local cluster result:

>tail -n 3 Chimeric.out.junction
chr2    84906088        +       chr2    235882419       -       0       0       1       A00499:1296:H7HFJDSXF:4:2671:11641:20306:TAGAGAAGG      84906002        8S86M24S    235882395       24M94S  1       118     96      103     103     1       GRPundef
# 2.7.11a   /usr/local/bin/STAR --genomeDir /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/genome_dir/ctat_genome_lib_build_dir/ref_genome.fa.star.idx --outReadsUnmapped None --chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutJunctionFormat 1 --alignSJDBoverhangMin 10 --alignMatesGapMax 100000 --alignIntronMax 100000 --alignSJstitchMismatchNmax 5 -1 5 5 --runThreadN 12 --outSAMstrandField intronMotif --outSAMunmapped Within --alignInsertionFlush Right --alignSplicedMateMapLminOverLmate 0 --alignSplicedMateMapLmin 30 --outSAMtype BAM Unsorted --readFilesIn /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/TL-25-ZZ9D50NI4G_T_RSQ1_1.fastq.gz /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/TL-25-ZZ9D50NI4G_T_RSQ1_3.fastq.gz --outSAMattrRGline ID:GRPundef --chimMultimapScoreRange 3 --chimScoreJunctionNonGTAG -4 --chimMultimapNmax 20 --chimOutType Junctions WithinBAM --chimNonchimScoreDropMin 10 --peOverlapNbasesMin 12 --peOverlapMMp 0.1 --genomeLoad NoSharedMemory --twopassMode None --readFilesCommand "gunzip -c" --quantMode GeneCounts
# Nreads 37358191       NreadsUnique 34233825   NreadsMulti 2176332

Terra result:

>zcat TL-25-ZZ9D50NI4G_Chimeric.out.junction.gz
chr_donorA      brkpt_donorA    strand_donorA   chr_acceptorB   brkpt_acceptorB strand_acceptorB        junction_type   repeat_left_lenA        repeat_right_lenB  read_name        start_alnA      cigar_alnA      start_alnB      cigar_alnB      num_chim_aln    max_poss_aln_score      non_chim_aln_score      this_chim_aln_scorebestall_chim_aln_score   PEmerged_bool   readgrp
# 2.7.11a   /usr/local/bin/STAR --genomeDir /mnt/disks/cromwell_root/genome_dir/ctat_genome_lib_build_dir/ref_genome.fa.star.idx --outReadsUnmapped None --chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutJunctionFormat 1 --alignSJDBoverhangMin 10 --alignMatesGapMax 100000 --alignIntronMax 100000 --alignSJstitchMismatchNmax 5 -1 5 5 --runThreadN 12 --outSAMstrandField intronMotif --outSAMunmapped Within --alignInsertionFlush Right --alignSplicedMateMapLminOverLmate 0 --alignSplicedMateMapLmin 30 --outSAMtype BAM Unsorted --readFilesIn /mnt/disks/cromwell_root/TL-25-ZZ9D50NI4G_T_RSQ1_1.fastq.gz /mnt/disks/cromwell_root/TL-25-ZZ9D50NI4G_T_RSQ1_3.fastq.gz --outSAMattrRGline ID:GRPundef --chimMultimapScoreRange 3 --chimScoreJunctionNonGTAG -4 --chimMultimapNmax 20 --chimOutType Junctions WithinBAM --chimNonchimScoreDropMin 10 --peOverlapNbasesMin 12 --peOverlapMMp 0.1 --genomeLoad NoSharedMemory --twopassMode None --readFilesCommand "gunzip -c" --quantMode GeneCounts
# Nreads 37358191       NreadsUnique 34233825   NreadsMulti 2176332

After removing all the path difference between two command lines, they are identical. I am wondering what would cause this problem. Thanks.

Best,

Tiger

[brianjohnhaas]

Hi, That is very peculiar. We have a WDL here in case you want to try that: https://github.com/STAR-Fusion/STAR-Fusion/blob/master/WDL/star_fusion_workflow.wdl We also have Terra set up pre-configured for either hg19 or hg38 in case you want to try that and compare: https://app.terra.bio/#workspaces/ctat-firecloud/STAR-Fusion/workflows hope this helps

…

On Thu, Dec 18, 2025 at 12:31 PM Quanhu Sheng ***@***.***> wrote: *shengqh* created an issue (STAR-Fusion/STAR-Fusion#437) <#437> I wrote a WDL for STAR-Fusion: https://github.com/shengqh/warp/blob/develop/pipelines/vumc_biostatistics/rnaseq/VUMCStarFusion.wdl. I tested the workflow on both local cluster and Terra using same raw FASTQ file, same reference GRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play.tar.gz and same parameters. The weird thing is, in the STAR step, the local cluster can generate correct Chimeric.out.junction while Terra one generate a file with header and comment only. Local cluster result: >tail -n 3 Chimeric.out.junction chr2 84906088 + chr2 235882419 - 0 0 1 A00499:1296:H7HFJDSXF:4:2671:11641:20306:TAGAGAAGG 84906002 8S86M24S 235882395 24M94S 1 118 96 103 103 1 GRPundef # 2.7.11a /usr/local/bin/STAR --genomeDir /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/genome_dir/ctat_genome_lib_build_dir/ref_genome.fa.star.idx --outReadsUnmapped None --chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutJunctionFormat 1 --alignSJDBoverhangMin 10 --alignMatesGapMax 100000 --alignIntronMax 100000 --alignSJstitchMismatchNmax 5 -1 5 5 --runThreadN 12 --outSAMstrandField intronMotif --outSAMunmapped Within --alignInsertionFlush Right --alignSplicedMateMapLminOverLmate 0 --alignSplicedMateMapLmin 30 --outSAMtype BAM Unsorted --readFilesIn /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/TL-25-ZZ9D50NI4G_T_RSQ1_1.fastq.gz /cromwell-executions/VUMCStarFusion/4267711a-d4b4-42b9-88de-29980b41062d/call-STARFusion/execution/TL-25-ZZ9D50NI4G_T_RSQ1_3.fastq.gz --outSAMattrRGline ID:GRPundef --chimMultimapScoreRange 3 --chimScoreJunctionNonGTAG -4 --chimMultimapNmax 20 --chimOutType Junctions WithinBAM --chimNonchimScoreDropMin 10 --peOverlapNbasesMin 12 --peOverlapMMp 0.1 --genomeLoad NoSharedMemory --twopassMode None --readFilesCommand "gunzip -c" --quantMode GeneCounts # Nreads 37358191 NreadsUnique 34233825 NreadsMulti 2176332 Terra result: >zcat TL-25-ZZ9D50NI4G_Chimeric.out.junction.gz chr_donorA brkpt_donorA strand_donorA chr_acceptorB brkpt_acceptorB strand_acceptorB junction_type repeat_left_lenA repeat_right_lenB read_name start_alnA cigar_alnA start_alnB cigar_alnB num_chim_aln max_poss_aln_score non_chim_aln_score this_chim_aln_scorebestall_chim_aln_score PEmerged_bool readgrp # 2.7.11a /usr/local/bin/STAR --genomeDir /mnt/disks/cromwell_root/genome_dir/ctat_genome_lib_build_dir/ref_genome.fa.star.idx --outReadsUnmapped None --chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutJunctionFormat 1 --alignSJDBoverhangMin 10 --alignMatesGapMax 100000 --alignIntronMax 100000 --alignSJstitchMismatchNmax 5 -1 5 5 --runThreadN 12 --outSAMstrandField intronMotif --outSAMunmapped Within --alignInsertionFlush Right --alignSplicedMateMapLminOverLmate 0 --alignSplicedMateMapLmin 30 --outSAMtype BAM Unsorted --readFilesIn /mnt/disks/cromwell_root/TL-25-ZZ9D50NI4G_T_RSQ1_1.fastq.gz /mnt/disks/cromwell_root/TL-25-ZZ9D50NI4G_T_RSQ1_3.fastq.gz --outSAMattrRGline ID:GRPundef --chimMultimapScoreRange 3 --chimScoreJunctionNonGTAG -4 --chimMultimapNmax 20 --chimOutType Junctions WithinBAM --chimNonchimScoreDropMin 10 --peOverlapNbasesMin 12 --peOverlapMMp 0.1 --genomeLoad NoSharedMemory --twopassMode None --readFilesCommand "gunzip -c" --quantMode GeneCounts # Nreads 37358191 NreadsUnique 34233825 NreadsMulti 2176332 After removing all the path difference between two command lines, they are identical. I am wondering what would cause this problem. Thanks. Best, Tiger — Reply to this email directly, view it on GitHub <#437>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABZRKX5UY24E6R3MFAZJG5D4CLQHJAVCNFSM6AAAAACPORKQYOVHI2DSMVQWIX3LMV43ASLTON2WKOZTG42DIMBWHE3TCNA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- -- Brian J. Haas The Broad Institute http://broadinstitute.org/~bhaas <http://broad.mit.edu/~bhaas>

[shengqh]

Hi,

That is very peculiar.

We have a WDL here in case you want to try that:
https://github.com/STAR-Fusion/STAR-Fusion/blob/master/WDL/star_fusion_workflow.wdl

We also have Terra set up pre-configured for either hg19 or hg38 in case
you want to try that and compare:
https://app.terra.bio/#workspaces/ctat-firecloud/STAR-Fusion/workflows

hope this helps
…

I actually borrow a lot of code from your wonderful workflow to fit my data. I have _1.fastq.gz and _3.fastq.gz in my tar.gz file which forbid me to use your workflow directly. I implemented a new workflow:

https://github.com/shengqh/warp/blob/develop/pipelines/vumc_biostatistics/rnaseq/VUMCStarFusion.wdl

1. Add Untar task to assign the fastq to fastq1 and fastq2
2. Call your star_fusion task to generate fusion result
3. Add MoveOrCopyFiles task to copy result to target gcp folder if it is defined.

I still didn't figure out why your workflow worked but my previous one didn't. Anyway, I tested the new workflow on 8 samples with docker "trinityctat/starfusion:1.15.1" and reference "GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play.tar.gz". Since the workflow worked on all samples, I will move on.

BTW, when I was testing the code below, it didn't work.

    if [[ ! -z "~{fastq_pair_tar_gz}" ]]; then
        # untar the fq pair
        mv ~{fastq_pair_tar_gz} reads.tar.gz
        tar xvf reads.tar.gz

        left_fq=(*_1.fastq* *_1.fq*)
    
        if [[ ! -z "*_2.fastq*" ]] || [[ ! -z "*_2.fq*" ]]; then
             right_fq=(*_2.fastq* *_2.fq*)
        fi
    else
        left_fq="~{left_fq}"
        right_fq="~{right_fq}"
    fi


    if [[ -z "${left_fq[0]}" && -z "${right_fq[0]}" ]]; then
        echo "Error, not finding fastq files here"
        ls -ltr
        exit 1
    fi

I modified it in my testing code and it worked. Hope it helps.

set -ex
shopt -s nullglob

if [[ ! -z "~{fastq_pair_tar_gz}" ]]; then
  # untar the fq pair
  mv ~{fastq_pair_tar_gz} reads.tar.gz
  tar xzvf reads.tar.gz
  rm reads.tar.gz

  # left reads
  if compgen -G "*_1.fastq*" > /dev/null; then
    left_fq=(*_1.fastq*)
  elif compgen -G "*_1.fq*" > /dev/null; then
    left_fq=(*_1.fq*)
  fi

  # right reads (_2 preferred, fallback to _3)
  if compgen -G "*_2.fastq*" > /dev/null; then
    right_fq=(*_2.fastq*)
  elif compgen -G "*_2.fq*" > /dev/null; then
    right_fq=(*_2.fq*)
  elif compgen -G "*_3.fastq*" > /dev/null; then
    right_fq=(*_3.fastq*)
  elif compgen -G "*_3.fq*" > /dev/null; then
    right_fq=(*_3.fq*)
  fi
else
  left_fq="~{left_fq}"
  right_fq="~{right_fq}"
fi

# sanity check
if [[ -z "${left_fq[0]:-}" || -z "${right_fq[0]:-}" ]]; then
  echo "Error: fastq files not found"
  ls -ltr
  exit 1
fi

echo "left_fq:  ${left_fq[@]}"
echo "right_fq: ${right_fq[@]}"