Quantifying expression from RNA-seq data

Takes raw FASTQ files through quality control, adaptor trimming, STAR alignment and duplicate marking to gene-level and transcript-level expression estimates.

Overview

This module follows the GTEx and TOPMed RNA-seq pipeline. Paired-end or single-end fastq.gz files go in; per-sample expression tables and a merged cohort matrix come out.

The recommended order is: fastqc to inspect read quality, fastp_trim_adaptor to remove adaptors if the reports call for it, STAR_align to map reads to the genome, then rnaseqc_call and rsem_call to quantify. Alignment and duplicate marking are combined: picard_qc is declared as the second step of STAR_align, so running the alignment also marks duplicates and collects Picard metrics. The alignment step reimplements the GTEx run_STAR.py wrapper inline rather than calling it, so that read length may differ between samples, and re-sorts the unsorted STAR output with samtools instead of letting STAR sort in memory.

Gene-level expression comes from RNA-SeQC v2.4.2 against a gene model collapsed so that each gene is a single feature rather than a set of transcripts DeLuca et al. (2012). Its read filters are strict: a read counts only if it is uniquely mapped (mapping quality 255 for STAR BAMs), aligned in a proper pair, carries no more than six non-reference bases, and falls entirely within exon boundaries. Reads overlapping introns are discarded. Exon-level counts are produced as well, and a read spanning several exons contributes to each in proportion to the share of the read it covers. Transcript-level expression comes from RSEM v1.3.0 against a transcript-level GTF.

Strandedness can be given per sample in the manifest as rf, fr, unstranded or strand_missing, following the library types of Signal et al. (2022); see the strand settings reference for what these mean. When it is not supplied, the strand_detected workflow infers it from the gene count table STAR emits. Read length is likewise per sample; a length of zero falls back to 100, which --sjdbOverhang can change.

Three workflows sit outside that main path. bam_to_fastq recovers FASTQ from aligned BAM when only BAM is available. trimmomatic_trim_adaptor is an alternative to fastp. filter_reads filters the alignments for allele-specific work: WASP correction drops reads carrying WASP flags, which mitigates reference-allele bias, and unique-read filtering keeps only alignments at mapping quality 255, following GTEx practice. Either can be applied alone or both together.

When to run it. At the start of the expression branch, on raw sequencing data, and before bulk_expression_QC and bulk_expression_normalization, which take the merged matrices produced here.

Input

• --sample-list – tab-delimited manifest with one line per sample and a header. ID and fq1 are required, fq2 is required for paired-end data, and strand and read_length are optional:

  ID	fq1	fq2	strand	read_length
  SAMPLE_001	SAMPLE_001_R1.fastq.gz	SAMPLE_001_R2.fastq.gz	rf	100
  SAMPLE_002	SAMPLE_002_R1.fastq.gz	SAMPLE_002_R2.fastq.gz	rf	100
  

  strand is rf, fr, unstranded or strand_missing; leaving it out lets strand_detected infer it. read_length of 0 means unknown and falls back to 100. If no sample has a known read length, omit the column entirely rather than filling it with zeros.

• --data-dir – directory holding the FASTQ files named in fq1 and fq2. The manifest carries bare filenames, so this has to point at the directory that actually contains them.

• --cwd – output directory.

Reference data, built by reference_data_preparation:

• --STAR-index – STAR genome index, required by STAR_align.

• --RSEM-index – RSEM reference, required by rsem_call.

• --gtf – gene model. STAR_align and picard_qc take the transcript-level GTF, while rnaseqc_call takes the collapsed gene-level one; passing the wrong one silently changes what is quantified.

• --reference-fasta – the genome FASTA the STAR index was built from. It has to be exactly that file; a different copy of the same assembly makes Picard fail.

• --ref-flat – refFlat annotation for Picard CollectRnaSeqMetrics, which reports how bases distribute across UTRs, introns, intergenic regions and coding exons. reference_data_preparation generates it in its RefFlat_generation step.

• --bam-list – the BAM manifest written by STAR_align, required by rnaseqc_call and rsem_call.

Adaptor trimming, for fastp_trim_adaptor and trimmomatic_trim_adaptor:

• --min-len – reads shorter than this after trimming are discarded; 15 under fastp, 50 under Trimmomatic.

• --window-size – sliding window for quality trimming, default 4.

• --required-quality – quality floor within the window, default 20.

• --leading / --trailing – mean quality below which the leading or trailing window is cut. fastp defaults to 20 and Trimmomatic to 3, so the fastp setting is by far the stricter, and matches cutadapt.

• --fasta-with-adapters-etc – adaptor reference. fastp infers adaptors from the reads and does not need one; Trimmomatic does, for instance TruSeq3-PE.fa from the Trimmomatic repository, chosen on the evidence of the FastQC overrepresented-sequences report.

• --seed-mismatches, --palindrome-clip-threshold, --simple-clip-threshold – Trimmomatic Illumina-clip settings.

Frequently adjusted:

• --sjdbOverhang – splice junction overhang, default 100, i.e. read length minus one.

• --chimSegmentMin – minimum chimeric segment length, default 15; 0 disables chimeric detection.

• --mapping-quality – uniqueness cutoff, default 255.

• --detection-threshold (rnaseqc_call) and --max-frag-len, --estimate-rspd (rsem_call).

• --optical-distance (picard_qc) – default 100; patterned flowcells such as HiSeq X need 2500.

• --zap-raw-bam – delete the raw BAM once duplicates are marked.

• --varVCFfile – variant VCF used to prepare WASP filtering.

• --uncompressed – set when the FASTQ files are not gzipped.

The captured help below lists every option, including the fastp and Trimmomatic trimming settings and the full set of STAR arguments.

Runtime:

• --numThreads, --job-size, --walltime, --mem – threads and cluster resources.

• --container, --entrypoint – software environment.

• --modular-script-dir – location of the drivers, default code/script.

Output

• <sample>_fastqc.html and <sample>_fastqc.zip – per-sample FastQC reports.

• <fastq>.trimmed.fq.gz with <fastq>.trimmed.fq.html and .json – trimmed reads and the fastp report; <sample_list>.trimmed.txt is the manifest of trimmed files. A read whose mate did not survive trimming is discarded by default. trimmomatic_trim_adaptor instead keeps them, writing paired and unpaired FASTQ separately, four files per paired-end sample.

• <sample>.Aligned.sortedByCoord.out.bam and <sample>.Aligned.toTranscriptome.out.bam – STAR genome and transcriptome alignments. filter_reads writes a WASP-filtered ..._wasp_qc...bam alongside.

• <sample_list>.bam_file_list – the BAM manifest that rnaseqc_call and rsem_call consume through --bam-list:

  sample_id	strand	coord_bam_list	BW_list	SJ_list	trans_bam_list
  SAMPLE_001	rf	SAMPLE_001.Aligned.sortedByCoord.out_nowasp_noqc.md.bam	SAMPLE_001.Aligned.sor
  

• <sample>.alignment_summary_metrics and the merged <sample_list>.picard.aggregated_quality.metrics.tsv – Picard duplicate and alignment metrics; the marked BAM flags duplicates with 0x400, decimal 1024:

  Sample PF_READS PF_READS_ALIGNED PCT_PF_READS_ALIGNED PCT_RIBOSOMAL_BASES PCT_CODING_BASES P
  SAMPLE_002 200000 12498 0.06249 0.001012 0.039564 0.289782 0.270124 0.399518 0.329346 0.0184
  

• <sample>.rnaseqc.gene_tpm.gct.gz, .gene_readsCount.gct.gz and .exon_readsCount.gct.gz – RNA-SeQC gene- and exon-level quantifications, merged across samples into <bam_list>.rnaseqc.*.gct.gz, next to <sample>.metrics.tsv and per-transcript coverage statistics.

• <sample>.rsem.genes.results and <sample>.rsem.isoforms.results – RSEM per-sample estimates, with run statistics in <sample>.rsem.stat/ including the .cnt read-assignment counts. These merge into seven files, four from the isoform columns and three from the gene columns, among them <bam_list>.rsem_transcripts_expected_count.txt.gz and the matching TPM and gene-level matrices.

• <sample_list>.multiqc_report.html – MultiQC summary across all of the above.

Each step also writes .stdout and .stderr beside its output. Example results for the toy data are under output/rnaseq/, but the raw FASTQ files that produced them are not distributed in this repository.

Minimal Working Example

These commands use the two demo samples listed in input/rnaseq/protocol_example.rnaseq.fastq.list.txt, with the simplest alignment recipe, no WASP filtering and no unique-read filtering, which is what quantifying expression for eQTL analysis calls for. The main path is fastqc, fastp_trim_adaptor, STAR_align, rnaseqc_call, rsem_call, run in that order; the remaining blocks are optional or alternatives. The raw FASTQ files are not distributed in this repository, so --data-dir is shown as a placeholder and has to point at the directory holding the files the manifest names.

Convert BAM back to FASTQ (optional)

Only needed when the starting material is aligned BAM rather than FASTQ.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb bam_to_fastq \
--cwd output/rnaseq \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/bam

Quality control with fastqc

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb fastqc \
--cwd output/rnaseq/fastqc \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq

Trim adaptors with fastp (optional)

Generates trimmed FASTQ files and a new sample list pointing at them.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb fastp_trim_adaptor \
--cwd output/rnaseq \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq

Trim adaptors with Trimmomatic (alternative)

An alternative to the fastp step above. Run one or the other, not both.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb trimmomatic_trim_adaptor \
--cwd output/rnaseq \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq

Align reads with STAR

Aligns reads and produces the per-run BAM manifest (<sample-list-stem>.bam_file_list) that the downstream RSEM step consumes. The recipe below adds no WASP/unique-filter tags (fastest; suitable for eQTL expression quantification). Other recipes are available by adding --wasp yes and/or --unique yes together with a --varVCFfile.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb STAR_align \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq \
--STAR-index reference_data/STAR_Index \
--gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
--reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
--ref-flat reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat \
--chimSegmentMin 0 \
-J 2 --mem 16G --numThreads 2

Detect strandedness

Infers the library strand from the STAR gene counts, for samples whose manifest leaves strand empty. Run it after STAR_align.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb strand_detected \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq

Mark duplicates and collect Picard metrics

picard_qc is also the second step of STAR_align, so the alignment above has already run it. This command runs it on its own, for instance on BAMs aligned outside this pipeline.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb picard_qc \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq \
--gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
--ref-flat reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat \
--reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta

Gene-level expression with rnaseqc

Uses the collapsed gene-level GTF. Point --bam_list at the manifest produced by STAR_align and keep --cwd consistent with where the BAMs live.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb rnaseqc_call \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq \
--gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf \
--reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
--bam-list output/rnaseq/bam/protocol_example.rnaseq.fastq.list.bam_file_list

Transcript-level expression with RSEM

Uses the transcript-level GTF and the RSEM reference. As above, --bam_list is the STAR_align manifest and --cwd must match the BAM location (the manifest stores bare filenames that are resolved relative to --cwd).

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb rsem_call \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq \
--RSEM-index reference_data/RSEM_Index \
--bam-list output/rnaseq/bam/protocol_example.rnaseq.fastq.list.bam_file_list

Filter reads with WASP

Applies WASP filtering to the alignments, producing the _wasp_qc BAMs used for allele-specific analysis.

Timing: TBD (on toy dataset)

sos run pipeline/RNA_calling.ipynb filter_reads \
--cwd output/rnaseq/bam \
--sample-list input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastq

Command Interface

sos run pipeline/RNA_calling.ipynb -h

usage: sos run pipeline/RNA_calling.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  bam_to_fastq
  fastqc
  fastp_trim_adaptor
  trimmomatic_trim_adaptor
  STAR_align
  strand_detected
  picard_qc
  rnaseqc_call
  rsem_call
  filter_reads

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        The output directory for generated files.
  --sample-list VAL (as path, required)
                        Sample meta data list
  --sample-name  (as list)
                        Sample names to analyze
  --data-dir  path(f"{sample_list:d}")

                        Raw data directory, default to the same directory as
                        sample list
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --java-mem 6G
                        Memory for Java virtual mechine (`picard`)
  --numThreads 8 (as int)
                        Number of threads
  --varVCFfile ''
                        VarVCFfile for data preparation for wasp_filtering
  --[no-]uncompressed (default to False)
                        Whether the fasta/fastq file is compressed or not.

Sections
  bam_to_fastq:
  fastqc:
  fastp_trim_adaptor_1:
    Workflow Options:
      --window-size 4 (as int)
                        sliding window setting
      --required-quality 20 (as int)
      --leading 20 (as int)
                        the mean quality requirement option for cut_front
      --trailing 20 (as int)
                        the mean quality requirement option for cut_tail
      --min-len 15 (as int)
                        reads shorter than length_required will be discarded
      --fasta-with-adapters-etc . (as path)
                        Path to the reference adaptors
  fastp_trim_adaptor_2:
  trimmomatic_trim_adaptor:
    Workflow Options:
      --fasta-with-adapters-etc . (as path)
                        Illumina clip setting Path to the reference adaptors
      --seed-mismatches 2 (as int)
      --palindrome-clip-threshold 30 (as int)
      --simple-clip-threshold 10 (as int)
      --window-size 4 (as int)
                        sliding window setting
      --required-quality 20 (as int)
      --leading 3 (as int)
                        Other settings
      --trailing 3 (as int)
      --min-len 50 (as int)
  STAR_align_1:
    Workflow Options:
      --gtf VAL (as path, required)
                        Reference gene model
      --STAR-index VAL (as path, required)
                        STAR indexing file
      --outFilterMultimapNmax 20 (as int)
      --alignSJoverhangMin 8 (as int)
      --alignSJDBoverhangMin 1 (as int)
      --outFilterMismatchNmax 999 (as int)
      --outFilterMismatchNoverLmax 0.1 (as float)
      --alignIntronMin 20 (as int)
      --alignIntronMax 1000000 (as int)
      --alignMatesGapMax 1000000 (as int)
      --outFilterType BySJout
      --outFilterScoreMinOverLread 0.33 (as float)
      --outFilterMatchNminOverLread 0.33 (as float)
      --limitSjdbInsertNsj 1200000 (as int)
      --outSAMstrandField intronMotif
      --outFilterIntronMotifs None
      --alignSoftClipAtReferenceEnds Yes
      --quantMode TranscriptomeSAM GeneCounts (as list)
      --outSAMattrRGline ID:rg1 SM:sm1 (as list)
      --outSAMattributes NH HI AS nM NM ch (as list)
      --chimSegmentMin 15 (as int)
      --chimJunctionOverhangMin 15 (as int)
      --chimOutType Junctions WithinBAM SoftClip (as list)
      --chimMainSegmentMultNmax 1 (as int)
      --sjdbOverhang 100 (as int)
      --mapping-quality 255 (as int)
  strand_detected_1:
  strand_detected_2:
    Workflow Options:
      --strand ''
  picard_qc, STAR_align_2:
    Workflow Options:
      --gtf VAL (as path, required)
                        Reference gene model
      --ref-flat VAL (as path, required)
                        Path to flat reference file, for computing QC metric
      --reference-fasta VAL (as path, required)
                        The fasta reference file used to generate star index
      --optical-distance 100 (as int)
                        For the patterned flowcell models (HiSeq X), change to
                        2500
      --[no-]zap-raw-bam (default to False)
  STAR_align_3:
  rnaseqc_call_1:
    Workflow Options:
      --bam-list VAL (as path, required)
      --gtf VAL (as path, required)
                        Reference gene model
      --reference-fasta VAL (as path, required)
      --detection-threshold 5 (as int)
      --mapping-quality 255 (as int)
  rnaseqc_call_2:
    Workflow Options:
      --bam-list VAL (as path, required)
  rsem_call_1:
    Workflow Options:
      --bam-list VAL (as path, required)
      --RSEM-index VAL (as path, required)
      --max-frag-len 1000 (as int)
      --[no-]estimate-rspd (default to True)
  rsem_call_2:
    Workflow Options:
      --bam-list VAL (as path, required)
  rsem_call_3, rnaseqc_call_3:
  rsem_call_4, rnaseqc_call_4:
  filter_reads:
    Workflow Options:
      --unique ''
      --wasp ''
      --mapping-quality 255 (as int)

Workflow implementation

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# The output directory for generated files.
parameter: cwd = path("output")
# Sample meta data list
parameter: sample_list = path
# Sample names to analyze
parameter: sample_name = list() #input should be --sample_name sample1 sample2, if multiple samples
# Raw data directory, default to the same directory as sample list
parameter: data_dir = path(f"{sample_list:d}")
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Memory for Java virtual mechine (`picard`)
parameter: java_mem = "6G"
# Number of threads
parameter: numThreads = 8
# VarVCFfile for data preparation for wasp_filtering
parameter: varVCFfile = ""
from sos.utils import expand_size
cwd = path(f'{cwd:a}')
data_dir = path(f'{data_dir:a}')
sample_list = path(f'{sample_list:a}')
modular_script_dir = path(f'{modular_script_dir:a}')
## Whether the fasta/fastq file is compressed or not.
parameter: uncompressed = False
import os
import csv
## FIX: The way to get sample needs to be revamped to 1. Accomodate rf/fr as column 2. accomodate single end read (Only 2 fq/samples)
with open(sample_list) as _fh:
    _reader = csv.DictReader(_fh, delimiter="\t")
    _columns = list(_reader.fieldnames)
    sample_inv = list(_reader)
del _reader, _fh

# Align specific sample, to ensure correct grouping in subsequent steps, perform this step prior to extracting strand, read_length, and sample_id information
if len(sample_name)>0:
    print("Align sample",sample_name, 'only...')
    sample_inv = [r for r in sample_inv if r['ID'] in sample_name]

## Extract strand information if user have specified the strand
strand_inv = []
if "strand" in _columns:
    strand_inv = [r['strand'] for r in sample_inv]
    stop_if(not all([x in ["fr", "rf", "unstranded","strand_missing"] for x in strand_inv ]), msg = "strand columns should only include ``fr``, ``rf``, ``strand_missing`` or ``unstranded``")

## Extract read_length if user have specified read_length
read_length = [0] * len(sample_inv)
if "read_length" in _columns:
    read_length = [int(r['read_length']) for r in sample_inv]

## Read columns follow the documented fq1/fq2 contract; legacy one/two-read-column sample sheets are also accepted
file_columns = [x for x in ["fq1", "fq2"] if x in _columns]
if len(file_columns) == 0:
    legacy_file_columns = [x for x in _columns if x not in ("ID", "strand", "read_length")]
    stop_if(len(legacy_file_columns) == 0 or len(legacy_file_columns) > 2,
        msg = "Sample list must contain fq1 and optional fq2 columns, or a legacy ID + one/two read-column layout after removing strand/read_length")
    file_columns = legacy_file_columns

## Extract sample_id
sample_id = [r['ID'] for r in sample_inv]

## Get the file name for single/paired end data
file_inv = [r[c] for r in sample_inv for c in file_columns]
file_inv = [x for x in file_inv if x is not None and str(x).strip() != ""]

raw_reads = [f'{data_dir}/{x}' for x in file_inv]


for y in raw_reads:
        if not os.path.isfile(y):
            raise ValueError(f"File {y} does not exist")

if len(raw_reads) != len(set(raw_reads)):
        raise ValueError("Duplicated files are found (but should not be allowed) in sample file list")

# Is the RNA-seq data pair-end
is_paired_end = 0 if len(raw_reads) == len(sample_id) else 1 
from sos.utils import env
env.logger.info(f'Input samples are {"paired-end" if is_paired_end else "single-end"} sequences.')

Step 0. Convert BAM back to FASTQ

Recovers paired FASTQ from an aligned BAM, for datasets that arrive as BAM.

[bam_to_fastq]
input: raw_reads, group_by = 1
output: f'{cwd}/{_input:bn}.1.fastq',f'{cwd}/{_input:bn}.2.fastq'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh bam_to_fastq \
        --cwd "${cwd}" \
        --input-bam "${_input}" \
        --sorted-bam "${_output[0]:nn}.sorted.bam" \
        --output-fastq1 "${_output[0]}" \
        --output-fastq2 "${_output[1]}"

Step 1. QC before alignment

fastqc on each FASTQ. Paired-end mates are assessed separately, so a paired sample yields two reports.

[fastqc]
input: raw_reads, group_by =  1
_n = os.path.basename(str(_input[0]))
_base = next((_n[:-len(s)] for s in (".fastq.gz", ".fq.gz", ".fastq", ".fq", ".bam", ".sam", ".cram") if _n.endswith(s)), os.path.splitext(_n)[0])
output: f'{cwd}/{_base}_fastqc.html',f'{cwd}/{_base}_fastqc.zip' 
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh fastqc \
        --cwd "${cwd}" \
        --sample-list "${sample_list}" \
        --data-dir "${data_dir}" \
        --input-fastq "${_input}" \
        --numThreads ${numThreads}

Step 2. Remove adaptors with fastp

Documentation: fastp. fastp infers the adaptor from the reads themselves: assuming a single adaptor confined to read tails, it counts 10-mers over the first million reads, keeps the frequent ones as seeds once low-complexity artefacts are dropped, and extends them by a tree search. No adaptor reference is required.

[fastp_trim_adaptor_1]
# sliding window setting
parameter: window_size = 4
parameter: required_quality = 20
# the mean quality requirement option for cut_front
parameter: leading = 20
# the mean quality requirement option for cut_tail
parameter: trailing = 20
# reads shorter than length_required will be discarded
parameter: min_len = 15
# Path to the reference adaptors
parameter: fasta_with_adapters_etc = path(".")
warn_if(fasta_with_adapters_etc.is_file(),msg = "Use input fasta and adaptor detection of paired-end read was disabled" )

input: raw_reads, group_by = is_paired_end + 1 , group_with = "sample_id"
output: [f'{cwd}/{path(x):bn}.trimmed.fq.gz' for x in _input]
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash:  expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh fastp_trim_adaptor \
        --cwd ${cwd} \
        --paired-end ${is_paired_end} \
        --input-fastq1 ${_input[0]} ${f'--input-fastq2 {_input[1]}' if is_paired_end else ''} \
        --output-fastq1 ${_output[0]} ${f'--output-fastq2 {_output[1]}' if is_paired_end else ''} \
        --adapter-fasta ${fasta_with_adapters_etc} \
        --min-len ${min_len} --window-size ${window_size} \
        --required-quality ${required_quality} --leading ${leading} --trailing ${trailing} \
        --numThreads ${numThreads}

[fastp_trim_adaptor_2]
input: group_by = "all"
output: f'{cwd}/{sample_list:n}.trimmed.txt'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \
        --step fastp_manifest \
        --sample-list "${sample_list}" \
        --is-paired-end ${is_paired_end} \
        --input ${' '.join([f'"{x}"' for x in _input])} \
        --output "${_output}"

Step 2 alternative: remove adaptors with Trimmomatic

Documentation: Trimmomatic. Superseded by fastp, which catches adaptors Trimmomatic cannot detect and needs no adaptor reference. Kept for cases where a known reference is preferred.

[trimmomatic_trim_adaptor]
# Illumina clip setting
# Path to the reference adaptors
parameter: fasta_with_adapters_etc = path(".")
parameter: seed_mismatches = 2
parameter: palindrome_clip_threshold = 30
parameter: simple_clip_threshold = 10
# sliding window setting
parameter: window_size = 4
parameter: required_quality = 20
# Other settings
parameter: leading = 3
parameter: trailing = 3
parameter: min_len = 50
input: raw_reads, group_by = is_paired_end + 1 , group_with = "sample_id"
output: ([ f'{cwd}/{_sample_id}_paired_{_input[0]:bn}.gz', f'{cwd}/{_sample_id}_unpaired_{_input[0]:bn}.gz',  f'{cwd}/{_sample_id}_paired_{_input[1]:bn}.gz',f'{cwd}/{_sample_id}_unpaired_{_input[1]:bn}.gz' ] if is_paired_end else f'{cwd}/{_sample_id}_trimmed_{_input:bn}.gz')
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh trimmomatic_trim_adaptor \
        --cwd ${cwd} \
        --paired-end ${is_paired_end} \
        --java-mem ${java_mem} \
        --input-fastq1 ${_input[0]} ${f'--input-fastq2 {_input[1]}' if is_paired_end else ''} \
        ${f'--output-paired1 {_output[0]} --output-unpaired1 {_output[1]} --output-paired2 {_output[2]} --output-unpaired2 {_output[3]}' if is_paired_end else f'--output-se {_output}'} \
        --adapter-fasta ${fasta_with_adapters_etc} \
        --seed-mismatches ${seed_mismatches} --palindrome-clip-threshold ${palindrome_clip_threshold} \
        --simple-clip-threshold ${simple_clip_threshold} \
        --leading ${leading} --trailing ${trailing} \
        --window-size ${window_size} --required-quality ${required_quality} --min-len ${min_len}

Step 3. Alignment with STAR

Documentation: STAR, and the GTEx run_STAR.py wrapper this step reimplements inline. Produces both the genome-coordinate and the transcriptome BAM.

[STAR_align_1]
# Reference gene model
parameter: gtf = path
# STAR indexing file
parameter: STAR_index = path
gtf = path(f'{gtf:a}')
STAR_index = path(f'{STAR_index:a}')
parameter: outFilterMultimapNmax = 20 
parameter: alignSJoverhangMin = 8 
parameter: alignSJDBoverhangMin = 1 
parameter: outFilterMismatchNmax = 999 
parameter: outFilterMismatchNoverLmax = 0.1 #larger than Yang's group (outFilterMismatchNoverReadLmax 0.04)
parameter: alignIntronMin = 20 
parameter: alignIntronMax = 1000000 
parameter: alignMatesGapMax = 1000000 
parameter: outFilterType =  "BySJout" 
parameter: outFilterScoreMinOverLread = 0.33 
parameter: outFilterMatchNminOverLread = 0.33 
parameter: limitSjdbInsertNsj = 1200000 
parameter: outSAMstrandField = "intronMotif" 
parameter: outFilterIntronMotifs = "None" 
parameter: alignSoftClipAtReferenceEnds = "Yes" 
parameter: quantMode = ["TranscriptomeSAM", "GeneCounts"]
parameter: outSAMattrRGline = ["ID:rg1", "SM:sm1"]
parameter: outSAMattributes = ["NH", "HI", "AS", "nM", "NM", "ch"] #(Yang's group: outSAMattributes NH HI AS nM XS vW). vW added is varVCFfile is set
parameter: chimSegmentMin = 15 
parameter: chimJunctionOverhangMin = 15 
parameter: chimOutType = ["Junctions", "WithinBAM", "SoftClip"]
parameter: chimMainSegmentMultNmax = 1 
parameter: sjdbOverhang = 100
parameter: mapping_quality = 255
if int(mem.replace("G","")) <  40:
    print("Insufficent memory for STAR, changing to 40G")
    star_mem = '40G'
else:
    star_mem = mem
if varVCFfile:
    if not path(varVCFfile).is_file():
        raise FileNotFoundError(f"Cannot find varVCFfile ``{varVCFfile}``")
    print("Adding wasp filter in STAR alignment")
# This option is commented out because it will force the downstream analysis to use 40G, which significantlly slow down the process.
input: raw_reads, group_by = is_paired_end + 1, group_with = {"sample_id","read_length"}
output: cord_bam = f'{cwd}/{_sample_id}.Aligned.sortedByCoord.out.bam',
        trans_bam = f'{cwd}/{_sample_id}.Aligned.toTranscriptome.out.bam'
if _read_length == 0:
    print("Using specified --sjdbOverhang as read length")
else:
    print("Using read length specified in the sample list")
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = star_mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh star_align_1 \
        --cwd "${cwd}" \
        --gtf "${gtf}" \
        --STAR-index "${STAR_index}" \
        --input-reads "${_input:r}" \
        --sample-id "${_sample_id}" \
        --paired-end ${is_paired_end} \
        --read-files-command "${"cat" if uncompressed else "zcat"}" \
        --output-prefix "${_output[0]:nnnn}" \
        --output-cord-bam "${_output["cord_bam"]}" \
        --output-trans-bam "${_output["trans_bam"]}" \
        --outFilterMultimapNmax ${outFilterMultimapNmax} \
        --alignSJoverhangMin ${alignSJoverhangMin} \
        --alignSJDBoverhangMin ${alignSJDBoverhangMin} \
        --outFilterMismatchNmax ${outFilterMismatchNmax} \
        --outFilterMismatchNoverLmax ${outFilterMismatchNoverLmax} \
        --alignIntronMin ${alignIntronMin} \
        --alignIntronMax ${alignIntronMax} \
        --alignMatesGapMax ${alignMatesGapMax} \
        --outFilterType "${outFilterType}" \
        --outFilterScoreMinOverLread ${outFilterScoreMinOverLread} \
        --outFilterMatchNminOverLread ${outFilterMatchNminOverLread} \
        --limitSjdbInsertNsj ${limitSjdbInsertNsj} \
        --outSAMstrandField "${outSAMstrandField}" \
        --outFilterIntronMotifs "${outFilterIntronMotifs}" \
        --alignSoftClipAtReferenceEnds "${alignSoftClipAtReferenceEnds}" \
        --quantMode "${" ".join(quantMode)}" \
        --outSAMattrRGline "${" ".join(outSAMattrRGline)}" \
        --outSAMattributes "${" ".join(outSAMattributes)}" \
        --chimSegmentMin ${chimSegmentMin} \
        --chimJunctionOverhangMin ${chimJunctionOverhangMin} \
        --chimOutType "${" ".join(chimOutType)}" \
        --chimMainSegmentMultNmax ${chimMainSegmentMultNmax} \
        --sjdbOverhang ${sjdbOverhang if _read_length == 0 else _read_length} \
        --var-vcf-file "${varVCFfile}"

[strand_detected_1: shared="step_strand_detected"]
input: output_from("STAR_align_1")["trans_bam"], group_with = "sample_id"
import subprocess
strand_detected = subprocess.run(
    ["Rscript", f"{modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R",
     "--step", "detect_strand", "--input", f"{cwd}/{_sample_id}.ReadsPerGene.out.tab"],
    capture_output=True, text=True, check=True).stdout.strip()
print(f"Detected strand for sample {_sample_id}: {strand_detected}")

[strand_detected_2: shared="strand"]
input: group_by = "all"
parameter: strand = ""
if not strand:
    if len(strand_inv) > 0:
        strand = strand_inv
        for i in range(0,len(strand)):
            if strand[i] == "strand_missing":
                strand[i] = step_strand_detected[i]
        print(f'Using strand specified in the input samples list {strand}, replacing strand_missing with detected strand')
    else:
        warn_if(not all(x is step_strand_detected[0] for x in step_strand_detected), msg = "strands detected are different among samples, please check your protocol, we will use the detected strand for each samples")    
        strand = step_strand_detected
        print(f'Using detected strand for each samples {strand}')
else:
    stop_if(strand not in ["fr", "rf", "unstranded"], msg = "``--strand`` option should be ``fr``, ``rf`` or ``unstranded``")
    print(f'Using ``--strand`` overwrite option for all the samples {strand[0]}') 
    strand = [strand] * len(step_strand_detected)

Step 4. Mark duplicates and collect metrics with Picard

The first QC pass after alignment: collects Picard RNA-seq metrics and rewrites the BAM with duplicates flagged.

[picard_qc, STAR_align_2]
# Reference gene model
parameter: gtf = path
depends: sos_variable('strand')
# Path to flat reference file, for computing QC metric
parameter: ref_flat = path
# The fasta reference file used to generate star index
parameter: reference_fasta = path
gtf = path(f'{gtf:a}')
ref_flat = path(f'{ref_flat:a}')
reference_fasta = path(f'{reference_fasta:a}')
# For the patterned flowcell models (HiSeq X), change to 2500
parameter: optical_distance = 100
parameter: zap_raw_bam = False
picard_strand_dict = {"rf":"SECOND_READ_TRANSCRIPTION_STRAND","fr": "FIRST_READ_TRANSCRIPTION_STRAND","unstranded":"NONE" }
input: output_from('filter_reads'), group_by = 2, group_with = {"sample_id","strand"}
output: picard_metrics = f'{cwd}/{_sample_id}.alignment_summary_metrics',
        picard_rna_metrics = f'{cwd}/{_sample_id}.rna_metrics',
        md_bam = f'{_input[0]:n}.md.bam',
        md_metrics = f'{_input[0]:n}.md.metrics',
        bigwig = f'{_input[0]:n}.md.bw',
        output_summary = f'{cwd}/{_sample_id}.bam_file_meta'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh picard_qc_star_align_2 \
        --cwd "${cwd}" \
        --gtf "${gtf}" \
        --ref-flat "${ref_flat}" \
        --reference-fasta "${reference_fasta}" \
        --java-mem "${java_mem}" \
        --optical-distance ${optical_distance} \
        --input-bam "${_input[0]}" \
        --sample-id "${_sample_id}" \
        --strand "${_strand}" \
        --picard-metrics "${_output["picard_metrics"]}" \
        --picard-rna-metrics "${_output["picard_rna_metrics"]}" \
        --md-bam "${_output["md_bam"]}" \
        --md-metrics "${_output["md_metrics"]}" \
        --bigwig "${_output["bigwig"]}" \
        --output-summary "${_output["output_summary"]}" \
        --var-vcf-file "${varVCFfile}" \
        --zap-raw-bam "${zap_raw_bam}"

[STAR_align_3]
input: group_by = "all"
depends: sos_variable("strand")
output: f'{cwd}/{sample_list:bn}.bam_file_list'
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = 1
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \
        --step star_align_3 \
        --output "${_output}" \
        --input ${" ".join([str(x) for x in _input])} \
        --sample-id ${" ".join([str(x) for x in sample_id])} \
        --strand ${" ".join([str(x) for x in strand])} \
        --var-vcf-file "${varVCFfile}"

Step 5. Post-alignment QC and quantification with RNA-SeQC

Documentation: RNA-SeQC, and the GTEx run_rnaseqc.py. The second QC pass, which also produces the gene- and exon-level quantifications.

[rnaseqc_call_1]
import os
import csv
parameter: bam_list = path
# Reference gene model
parameter: gtf = path
parameter: reference_fasta = path
bam_list = path(f'{bam_list:a}')
gtf = path(f'{gtf:a}')
reference_fasta = path(f'{reference_fasta:a}')
with open(bam_list) as _fh:
    sample_inv = list(csv.DictReader(_fh, delimiter="\t"))
parameter: detection_threshold = 5
parameter: mapping_quality = 255
## Extract strand information if user have specified the strand
strand_list = [r['strand'] for r in sample_inv]
stop_if(not all([x in ["fr", "rf", "unstranded"] for x in strand_list ]), msg = "strand columns should only include ``fr``, ``rf`` or ``unstranded``, please check the bam_list")     
## Extract sample_id
sample_id = [r['sample_id'] for r in sample_inv]
    
## Get the file name for cood_bam data
coord_bam_list_inv = [r['coord_bam_list'] for r in sample_inv]
coord_bam_list = [x if os.path.isabs(x) else f'{cwd}/{x}' for x in coord_bam_list_inv ]
input:  coord_bam_list, group_by = 1, group_with = {"sample_id","strand_list"}
output: f'{cwd}/{_sample_id}.rnaseqc.gene_tpm.gct.gz',
        f'{cwd}/{_sample_id}.rnaseqc.gene_reads.gct.gz',
        f'{cwd}/{_sample_id}.rnaseqc.exon_reads.gct.gz',
        f'{cwd}/{_sample_id}.rnaseqc.metrics.tsv'

task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh rnaseqc_call \
        --cwd "${cwd}" \
        --sample-list "${sample_list}" \
        --data-dir "${data_dir}" \
        --gtf "${gtf}" \
        --reference-fasta "${reference_fasta}" \
        --input-bam "${_input}" \
        --sample-id "${_sample_id}" \
        --strand "${_strand_list}" \
        --detection-threshold "${detection_threshold}" \
        --mapping-quality "${mapping_quality}" \
        --paired-end ${is_paired_end} \
        --numThreads ${numThreads}

The per-sample RNA-SeQC results are merged in the step below.

[rnaseqc_call_2]
parameter: bam_list = path
input: group_by = "all"
output: f'{cwd}/{bam_list:bn}.rnaseqc.gene_tpm.gct.gz',
        f'{cwd}/{bam_list:bn}.rnaseqc.gene_readsCount.gct.gz',
        f'{cwd}/{bam_list:bn}.rnaseqc.exon_readsCount.gct.gz',
        f'{cwd}/{bam_list:bn}.rnaseqc.metrics.tsv'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \
        --step rnaseqc_merge \
        --cwd "${cwd}" \
        --name "${bam_list:bn}" \
        --input ${" ".join([str(x) for x in _input])}

Step 6. Quantify expression with RSEM

Documentation: RSEM, and the GTEx run_RSEM.py. Estimates gene- and isoform-level expression from the transcriptome BAM.

[rsem_call_1]
parameter: bam_list = path
parameter: RSEM_index = path
bam_list = path(f'{bam_list:a}')
RSEM_index = path(f'{RSEM_index:a}')
parameter: max_frag_len = 1000
parameter: estimate_rspd = True

import os
import csv
with open(bam_list) as _fh:
    sample_inv = list(csv.DictReader(_fh, delimiter="\t"))

## Extract strand information if user have specified the strand
strand_list = [r['strand'] for r in sample_inv]
stop_if(not all([x in ["fr", "rf", "unstranded"] for x in strand_list ]), msg = "strand columns should only include ``fr``, ``rf`` or ``unstranded``, please check the bam_list")     
## Extract sample_id
sample_id = [r['sample_id'] for r in sample_inv]
    
## Get the file name for trans_bam_list data
trans_bam_list_inv = [r['trans_bam_list'] for r in sample_inv]
trans_bam_list = [x if os.path.isabs(x) else f'{cwd}/{x}' for x in trans_bam_list_inv]
input: trans_bam_list, group_by = 1, group_with = {"sample_id", "strand_list"} 
output: f'{cwd}/{_sample_id}.rsem.isoforms.results', f'{cwd}/{_sample_id}.rsem.genes.results',f'{cwd}/{_sample_id}.rsem.stat/{_sample_id}.rsem.cnt'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, trunk_size = job_size
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    run_RSEM ${RSEM_index:a} ${_input:a} ${_sample_id} \
        -o ${_output[0]:d} \
        --max_frag_len ${max_frag_len} \
        --estimate_rspd ${'true' if estimate_rspd else 'false'} \
        --paired_end ${"true" if is_paired_end else "false"} \
        --is_stranded ${"true" if _strand_list != "unstranded" else "false"} \
        --threads ${numThreads}

The per-sample RSEM results are merged in the steps below.

[rsem_call_2]
parameter: bam_list = path
input: group_by = "all"
output: f'{cwd}/{bam_list:bn}.rsem_transcripts_expected_count.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_transcripts_tpm.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_transcripts_fpkm.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_transcripts_isopct.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_genes_expected_count.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_genes_tpm.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem_genes_fpkm.txt.gz',
        f'{cwd}/{bam_list:bn}.rsem.aggregated_quality.metrics.tsv'

task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, trunk_size = job_size
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \
        --step rsem_call_2 \
        --cwd "${cwd}" \
        --name "${bam_list:bn}" \
        --input ${" ".join([str(x) for x in _input])}

Step 7. Summarize with MultiQC

MultiQC searches the output directory, recognises the logs and summary files the earlier steps left behind, and writes a single HTML report plus a directory of parsed data. Pointing it at the directory holding all outputs is enough.

[rsem_call_3,rnaseqc_call_3]
output: f'{_input[0]:nnn}.multiqc_report.html'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, trunk_size = job_size
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh multiqc_report \
        --cwd "${cwd}" \
        --input-dir "${_input:d}" \
        --output-report "${_output}" \
        --multiqc-config "${_output:n}.multiqc_config.yml"

[rsem_call_4, rnaseqc_call_4]
# Path to flat reference file, for computing QC metric
output: f'{_input[0]:nn}.picard.aggregated_quality.metrics.tsv'# it will be outputed in a very 'deep' directory with {cwd}/
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, trunk_size = job_size
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.R \
        --step aggregate_picard_qc \
        --input-dir "${_input[-1]:d}" \
        --output "${_output}" \
        --is-paired-end ${is_paired_end} \
        ${"--wasp" if varVCFfile else ""} \
        --numThreads ${numThreads}

Step 8. Filter reads with WASP

Optional. It reads the STAR_align_1 output directly, so no input has to be named.

[filter_reads]
parameter: unique=""
parameter: wasp=""
parameter: mapping_quality = 255
depends: sos_variable('strand')
input: output_from("STAR_align_1"), group_by = 2, group_with = {"sample_id","strand"}
output: cord_bam_wasp_qc = f'{cwd}/{_sample_id}.Aligned.sortedByCoord.out{"_wasp" if wasp else "_nowasp"}{"_qc" if wasp else "_noqc"}.bam',
        trans_bam_wasp_qc = f'{cwd}/{_sample_id}.Aligned.toTranscriptome.out{"_wasp" if wasp else "_nowasp"}{"_qc" if wasp else "_noqc"}.bam'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/molecular_phenotypes/calling/RNA_calling.sh filter_reads \
        --cwd "${cwd}" \
        --input-cord-bam "${_input[0]}" \
        --input-trans-bam "${_input[1]}" \
        --output-cord-bam "${_output["cord_bam_wasp_qc"]}" \
        --output-trans-bam "${_output["trans_bam_wasp_qc"]}" \
        --unique "${unique}" \
        --wasp "${wasp}" \
        --mapping-quality "${mapping_quality}"

Troubleshooting

Step	Substep	Problem	Possible Reason	Solution
rsem_call / rnaseqc_call	-	Target unavailable: .../<name>.Aligned.toTranscriptome.out_*.bam	--bam_list points at a manifest whose bare filenames do not exist under the given --cwd (the step prepends {cwd}/ to each entry)	Use the STAR_align-generated <sample-list-stem>.bam_file_list and set --cwd to the directory that actually holds those BAMs