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.IDandfq1are required,fq2is required for paired-end data, andstrandandread_lengthare 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 100strandisrf,fr,unstrandedorstrand_missing; leaving it out letsstrand_detectedinfer it.read_lengthof0means 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 infq1andfq2. 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 bySTAR_align.--RSEM-index-- RSEM reference, required byrsem_call.--gtf-- gene model.STAR_alignandpicard_qctake the transcript-level GTF, whilernaseqc_calltakes 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 PicardCollectRnaSeqMetrics, which reports how bases distribute across UTRs, introns, intergenic regions and coding exons.reference_data_preparationgenerates it in itsRefFlat_generationstep.--bam-list-- the BAM manifest written bySTAR_align, required byrnaseqc_callandrsem_call.
Adaptor trimming, for fastp_trim_adaptor and trimmomatic_trim_adaptor:
--min-len-- reads shorter than this after trimming are discarded;15under fastp,50under Trimmomatic.--window-size-- sliding window for quality trimming, default4.--required-quality-- quality floor within the window, default20.--leading/--trailing-- mean quality below which the leading or trailing window is cut. fastp defaults to20and Trimmomatic to3, 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 instanceTruSeq3-PE.fafrom 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, default100, i.e. read length minus one.--chimSegmentMin-- minimum chimeric segment length, default15;0disables chimeric detection.--mapping-quality-- uniqueness cutoff, default255.--detection-threshold(rnaseqc_call) and--max-frag-len,--estimate-rspd(rsem_call).--optical-distance(picard_qc) -- default100; patterned flowcells such as HiSeq X need2500.--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, defaultcode/script.
Output¶
<sample>_fastqc.htmland<sample>_fastqc.zip-- per-sample FastQC reports.<fastq>.trimmed.fq.gzwith<fastq>.trimmed.fq.htmland.json-- trimmed reads and the fastp report;<sample_list>.trimmed.txtis the manifest of trimmed files. A read whose mate did not survive trimming is discarded by default.trimmomatic_trim_adaptorinstead keeps them, writing paired and unpaired FASTQ separately, four files per paired-end sample.<sample>.Aligned.sortedByCoord.out.bamand<sample>.Aligned.toTranscriptome.out.bam-- STAR genome and transcriptome alignments.filter_readswrites a WASP-filtered..._wasp_qc...bamalongside.<sample_list>.bam_file_list-- the BAM manifest thatrnaseqc_callandrsem_callconsume 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_metricsand the merged<sample_list>.picard.aggregated_quality.metrics.tsv-- Picard duplicate and alignment metrics; the marked BAM flags duplicates with0x400, 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.gzand.exon_readsCount.gct.gz-- RNA-SeQC gene- and exon-level quantifications, merged across samples into<bam_list>.rnaseqc.*.gct.gz, next to<sample>.metrics.tsvand per-transcript coverage statistics.<sample>.rsem.genes.resultsand<sample>.rsem.isoforms.results-- RSEM per-sample estimates, with run statistics in<sample>.rsem.stat/including the.cntread-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.gzand 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 tests/fixtures/rna_calling/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 tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/bamQuality control with fastqc¶
Timing: TBD (on toy dataset)
sos run pipeline/RNA_calling.ipynb fastqc \
--cwd output/rnaseq/fastqc \
--sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastqTrim 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 tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastqTrim 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 tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastqAlign 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 tests/fixtures/rna_calling/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 2Detect 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 tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastqMark 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 tests/fixtures/rna_calling/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.fastaGene-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 tests/fixtures/rna_calling/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_listTranscript-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 tests/fixtures/rna_calling/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_listFilter 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 tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
--data-dir path/to/fastqCommand Interface¶
sos run pipeline/RNA_calling.ipynb -husage: 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 |
- DeLuca, D. S., Levin, J. Z., Sivachenko, A., Fennell, T., Nazaire, M.-D., Williams, C., Reich, M., Winckler, W., & Getz, G. (2012). RNA-SeQC: RNA-seq metrics for quality control and process optimization. Bioinformatics, 28(11), 1530–1532. 10.1093/bioinformatics/bts196
- Signal, B., & Kahlke, T. (2022). how_are_we_stranded_here: quick determination of RNA-Seq strandedness. BMC Bioinformatics, 23(1). 10.1186/s12859-022-04572-7