bulk sQTL workflow for eQTLGen: 1. splicing junction

We implement the sQTL workflow in two stages—splicing junction calling and sQTL calling—to ensure robust quality control of junctions.

This is the cookbook for the first stage: splicing junction calling. In concordance with previous consortium efforts (e.g. GTEx, INTERVAL), we use the RegTools and LeafCutter pipeline. To harmonize junctions across all cohorts, filters for junction cluster (total reads) and junction (ratio in the group) are not applied in this workflow to generate per-cohort data. Instead, we will perform a stringent filtering based on the data from all cohorts and generate the harmonized splicing phenotypes for sQTL calling and all other downstream analysis.

Below is the cookbook for per-cohort splicing junction file ({cohort_name}_cluster_perind_numers.counts.gz) preparation.

For this cookbook, we follow the style of eQTLGen Phase II cookbook and aim to require minimal additional user actions.

Eligibility

If you have RNA-seq data on your cohort(s), you will be able to run this analysis. You can either start from FASTQ (recommended) or aligned BAM files.

Prerequisites

If you have run the workflow in eQTLGen Phase II cookbook, you should have all the following prerequisites. However, the current pipelines and configurations are tailored for use with the PBS scheduler by default. If your HPC uses a different scheduler, modify the job scripts to request the appropriate resources accordingly.

• High-Performance Computing (HPC) environment

• Singularity

• Nextflow

Getting Started

Setup

Make analysis root folder for this project, e.g. eQTLGen_splicing_phase1. Please ensure the folder is in data storage disk (not your home directory) to avoid issues.

mkdir eQTLGen_splicing_phase1

Inside eQTLGen_splicing_phase1, download the wrapper of pipelines for analysis:

cd eQTLGen_splicing_phase1
wget -O eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz "https://www.dropbox.com/scl/fi/0o40fueusy4s4mw23emd2/eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz?rlkey=xiuuddvolv9ia2d82r0yazpwx&st=fg905xpv&dl=0"
wget -O eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz.md5 "https://www.dropbox.com/scl/fi/p3hqpuz2e7d737g5g2p2y/eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz.md5?rlkey=gpx33ozddck3hrltgeidefd5z&st=641glw2v&dl=0"

Please check md5sum of the file and extract the content.

# This should print eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz: OK
md5sum -c eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz.md5 
tar -xzvf eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz

You should be able to see two folders: data and singularity_img. Below is the full content (displayed by command tree ./):

.
|-- 1_splicing_junction
|   |-- data
|   |   |-- GRCh38.primary_assembly.genome.fa.gz
|   |   |-- gencode.v48.annotation.bed12
|   |   `-- gencode.v48.annotation.gtf
|   |-- singularity_img
|   |   |-- eqtlgen_splicing_junction_strandness.sif
|   |   |-- eqtlgen_splicing_junction_v2.sif
|   |   `-- eqtlgen_splicing_star_v2.sif
|   `-- splicing_junction
|       |-- conf
|       |   |-- Danish_Computerome_profile.config
|       |   |-- JCTF_profile.config
|       |   |-- NSCC_profile.config
|       |   |-- base.config
|       |   |-- junction_processes.config
|       |   |-- local_vm.config
|       |   |-- lsf_per_core.config
|       |   |-- lsf_per_job.config
|       |   |-- pbs.config
|       |   |-- sge.config
|       |   |-- slurm.config
|       |   `-- star_processes.config
|       |-- nextflow.config
|       |-- splicing_junction.sum_stat.nf
|       |-- star_alignment_with_qc.nf
|       |-- submit_1_star_alignment_template.sh
|       |-- submit_2_strandness.sh
|       `-- submit_3_splicing_junction_pipeline_template.sh
|-- eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz
|-- eQTLGen_splicing_1_splicing_junction_wrapper.tar.gz.md5
`-- tools
    `-- nextflow

Test the nextflow executable. You need to load Java >=11 before testing (e.g. module load [Java >=11 module name in your HPC]). We encountered some issues when using Java 21, with various versions of nextflow, it should work, but if it doesn't please try and older version of Java (Java 17 worked perfectly)

cd tools
./nextflow run info
# If this binary works on your machine, you will see:
# N E X T F L O W   ~  version 25.10.2

(Skip this step if nextflow works) If the nextflow does not work, please download and self-install Nextflow executable and test. You might need to load Java >=11 before running self-install (e.g. module load [Java >=11 module name in your HPC]).

wget -qO- https://get.nextflow.io | bash
./nextflow run info

Enter the main work directory, and prepare to run the pipelines.

cd ../1_splicing_junction/splicing_junction

Inside pipeline folder (splicing_junction) is the script template using name format submit_*_template.sh e.g. eQTLGen_splicing_phase1/1_splicing_junction/splicing_junction/submit_1_star_alignment_template.sh. You should adjust this according to your data and environment (e.g. specify the path to your input folder, add required HPC modules), and save it. We specifically marked the sections by the tag: "USER CONFIGURATION REQUIRED".

To ease the process, we have pre-filled some paths in the templates, assuming that you use the recommended folder structure.

You should submit each pipeline from inside pipeline folder.

Data

There are two options as the starting point and a helper script (submit_2_strandness.sh) to infer strand-specificity of RNA-seq library (see details in Strand-specificity)):

• RNA-seq raw FASTQ files (paired-end): Run submit_1_star_alignment_[cohort name].sh, submit_2_strandness.sh (optional), and submit_3_splicing_junction_pipeline_[cohort name].sh.

• Aligned BAM files (see below for requirements of BAM file): Run submit_2_strandness.sh (optional) and submit_3_splicing_junction_pipeline_[cohort name].sh.

By default, we assume the genome assembly used for read alignment is hg38. You need to re-generate BAM files using this workflow if you are using a different version (e.g. hg19).

BAM files generated by STAR (recommended) can be used. Importantly, if your library is unstranded (see Strand-specificity)) if you are not sure), the BAM file must contain the XS tag on spliced reads. This is required for junction strand assignment in this workflow. To manually check if the XS tag is present in your BAM files, you can use samtools view [sample_name].bam | grep XS | head , which will display lines if XS is present.

Configure and run

Step 1: if you start from FASTQ files

Modify the job file submit_1_star_alignment_template.sh in the main work directory. You can search for sections marked by: "USER CONFIGURATION REQUIRED", and at least provide the following:

1. the absolute path of your input FASTQ directory.

2. read length / maximum read length of paired end sequencing read mates (e.g. 150 for pair-end 150).

3. -profile in the Nextflow command (at the bottom of the script) and corresponding config file.

sh submit_1_star_alignment_[cohort name].sh   # or equivalent command for your scheduler

Step 2: (Optional) Infer strandness from BAM files

If you are not sure about what strand parameters to use in the following splicing junction pipeline, run the helper script in Strand-specificity).

Step 3: If you start from aligned BAM files or you have generated BAM files as above

Modify the job file submit_3_splicing_junction_pipeline_template.sh in the main work directory. You can search for sections marked by: "USER CONFIGURATION REQUIRED", and at least provide the following:

1. the absolute path of your input BAM directory.

2. cohort name.

3. strand parameter.

4. -profile in the Nextflow command (at the bottom of the script) and corresponding config file.

sh submit_3_splicing_junction_pipeline_[cohort name].sh   # or equivalent command for your scheduler

Strand-specificity

If you are not sure whether your RNA-seq library is stranded or not, please run the short diagnostic script named submit_2_strandness.sh. This script runs on BAM files, so the first step is to run the alignment pipeline if the starting point is FASTQ files. We use infer_experiment.py from RSeQC.

# Please modify BAM_DIR before you run if you start from BAM files
# Adjust headers based on your scheduler
qsub submit_2_strandness.sh   # or equivalent command for your scheduler

You can check the results in prepare/strandness.txt after the job completes. This usually takes less than 1 minute. For paired-end RNA-seq:

   Usually, strand-specific RNA-seq libraries are prepared using a dUTP-based protocol, resulting in reverse-stranded reads (--strandness RF). In some cases (e.g. older datasets or non–strand-specific protocols), RNA-seq libraries may be unstranded (--strandness XS).
   If the fraction of reads explained by 1+-,1-+,2++,2-- is much larger than that explained by 1++,1--,2+-,2-+ and is close to 1             => RF (reverse-stranded, dUTP-based)
   If 1+-,1-+,2++,2-- is close to 1++,1--,2+-,2-+ and both are close to 0.5             => XS (unstranded)

   Below are two example outputs:

   Example 1: Paired-end, non–strand-specific: 
   This is PairEnd Data
   Fraction of reads failed to determine: 0.0172
   Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903
   Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925

   Example 2: Paired-end, strand-specific (reverse-stranded):
   This is PairEnd Data
   Fraction of reads failed to determine: 0.0072
   Fraction of reads explained by "1++,1--,2+-,2-+": 0.0487
   Fraction of reads explained by "1+-,1-+,2++,2--": 0.9441

If you encounter other scenarios or ambiguous results, please let us know. You can also refer to this table for more details.

Output

The pipeline will generate a output folder with junction files for analysis.

output/
├── 00a_qc_mapq
|   |-- {sample_name}.uniq.qc.txt
|   `-- ... .uniq.qc.txt
├── 01_junc
|   |-- {sample_name}.junc
|   `-- ... .junc
|-- 02_cluster
|   `-- {cohort_name}_cluster_perind_numers.counts.gz
|-- 03_summary_stats
|   `-- {cohort_name}_summary_stats.tsv
|-- bam_alignment
|   |-- 00_fastq_qc
|   |-- 00_star_genome_index
|   |-- aligned_bam
|   `-- star_alignment_qc_summary.tsv
...

Send output file to central analysis for joint junction quality control & filtering:

1. Please upload either via your preferred file-sharing service (e.g., Dropbox/Globus) and share with us the link.

• {cohort_name}_cluster_perind_numers.counts.gz (donor-level)

• {cohort_name}_summary_stats.tsv (summary statistics)

2. Please also submit relevant metadata (if available) of your cohort. For example, ancestry, sex and age can be included. This can be a table in either csv/tsv or txt format.

3. Please keep the {cohort_name}_cluster_perind_numers.counts.gz file, which will be used in stage 2 for prepration of phenotype.

Runtime, memory, and storage

Actual resource usage may vary depending on FASTQ/BAM file size, sequencing depth, read length, and filesystem performance.

(from FASTQ) Resource usage estimates below were obtained from test runs on the MAGE LCL subset (N = 327 samples), with a total input FASTQ (.gz) size of 1.4 TB. All benchmarks were performed on a shared HPC filesystem (NFS-based /scratch storage), which is optimized for high-throughput I/O and concurrent access across compute nodes.

Stage	Threads	Runtime	Total CPU-Hours	Peak Memory	Total Disk	Parallelizable
STAR Index (one-time)	16	25 min	6.67	42 GB	75 GB	No
FASTP_QC	4	~5 min/sample	~109	2 GB	~1.4 TB	Yes (327 parallel jobs)
STAR_ALIGN	16	~26 min/sample	~2267.2	42 GB	~7.1 TB	Yes (327 parallel jobs)
Junction Extraction	1	~12 min/sample	~65.4	0.5 GB	~6.5 GB	Yes (327 parallel jobs)
Junction Clustering	1	73 min	1.22	16 GB	4 GB	No (cohort-level)
Summary Stats	1	6 min	0.1	0.5 GB	0.2 GB	No (cohort-level)
TOTAL	–	–	~2449.6 CPU-hours	42 GB	~8.6 TB	Mostly parallel

Troubleshooting

You can check the scheduler-generated stdout/stderr log files (e.g. 1_splicing_junction.o* in PBS) or inspect .nextflow.log.

Please also refer to Section "Running, monitoring and debugging Nextflow pipelines" in eQTLGen Phase II cookbook.

Citation

Nextflow: Di Tommaso, P., Chatzou, M., Floden, E. et al. Nextflow enables reproducible computational workflows. Nat Biotechnol 35, 316–319 (2017). https://doi.org/10.1038/nbt.3820

RSeQC: Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012 Aug 15;28(16):2184-5. doi: 10.1093/bioinformatics/bts356. Epub 2012 Jun 27. PMID: 22743226.

Regtools: Cotto KC, Feng YY, Ramu A, Richters M, Freshour SL, Skidmore ZL, Xia H, McMichael JF, Kunisaki J, Campbell KM, Chen TH, Rozycki EB, Adkins D, Devarakonda S, Sankararaman S, Lin Y, Chapman WC, Maher CA, Arora V, Dunn GP, Uppaluri R, Govindan R, Griffith OL, Griffith M. Integrated analysis of genomic and transcriptomic data for the discovery of splice-associated variants in cancer. Nature Communications. 2023 Mar. pmid: 36949070; doi: 10.1038/s41467-023-37266-6

LeafCutter: Li YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, Pritchard JK. Annotation-free quantification of RNA splicing using LeafCutter. Nat Genet. 2018 Jan;50(1):151-158. doi: 10.1038/s41588-017-0004-9. Epub 2017 Dec 11. PMID: 29229983; PMCID: PMC5742080.

Contact

Please contact e1101920@u.nus.edu if you have any other questions.