nf-rna-wasp-allele-count

A Nextflow DSL2 pipeline to perform allele-specific expression analysis from RNA-seq data using a WASP-aware alignment strategy.

Note

This pipeline adopts the genotype-aware alignment strategy used in:

Taylor DJ, Chhetri SB, Tassia MG, et al. Sources of gene expression variation in a globally diverse human cohort. Nature 632, 122–130 (2024). https://doi.org/10.1038/s41586-024-07708-2

This workflow uses STAR alignment with WASP correction to prevent reference-mapping bias at heterozygous sites. While the original study focused on regulatory and splicing variation, this pipeline extends the same alignment principles to explicit allele-specific read counting from RNA-seq data.

Features

• Genotype-aware RNA-seq alignment with STAR + WASP correction.
• Optional allele-specific read counting (REF/ALT) from RNA-seq BAMs, performed only if a regions VCF is provided and exists.
• Reports total depth and variant allele fraction.

Requirements

• Nextflow ≥ 22.10.0
• Singularity / Apptainer or Docker (or compatible container runtime)
• Phased VCFs for each sample (e.g., 1000G high-coverage release)
• Reference genome FASTA + GTF matching alignment build

Quick Start

1. Prepare your sample sheet (CSV)

samplesheet.csv

sample,fastq_1,fastq_2
sampleID,sampleID_1.fastq.gz,sampleID_2.fastq.gz

2. Download reference data

Download the high-coverage phased WGS VCFs an store all the files in a directory:

https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased

Expected output:

1000G_2504_high_coverage_20201028_3202_phased_vcfs/
├── CCDG_14151_B01_GRM_WGS_2020-08-05_chr*.filtered.shapeit2-duohmm-phased.vcf.gz
├── CCDG_14151_B01_GRM_WGS_2020-08-05_chr*.filtered.shapeit2-duohmm-phased.vcf.gz.tbi
├── CCDG_14151_B01_GRM_WGS_2020-08-05_chrX.filtered.eagle2-phased.v2.vcf.gz
└── CCDG_14151_B01_GRM_WGS_2020-08-05_chrX.filtered.eagle2-phased.v2.vcf.gz.tbi

3. Reference genome and annotation

mkdir reference
cd reference
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gtf.gz
gunzip gencode.v38.annotation.gtf.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

Expected output:

reference/
├── GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
└── gencode.v38.annotation.gtf

4. Define SNPs of interest (VCF) (optional)

If you want to perform allele-specific counting, provide a VCF file with the SNPs of interest. If this file is not provided or does not exist, the allele counting step will be skipped.

query_sites.vcf

##fileformat=VCFv4.2
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
chr1    123 .   C   T   .   .  .

5. Set pipeline parameters

params {
    sample_sheet = samplesheet.csv
    phased_vcf_dir = 1000G_2504_high_coverage_20201028_3202_phased_vcfs
    reference_fa = GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
    gencode_gtf = gencode.v38.annotation.gtf
    regions_vcf = query_sites.vcf
    outdir = results/
}

6. Run the pipeline

nextflow run main.nf -profile docker

Workflow Overview

The pipeline consists of the following main steps:

• Adapter removal (ADAPTER_TRIM)

  • Used trimgalore to remove adapter sequences from the reads. Output: trimmed reads and fastqc report.

• Genome index generation (STAR_GENOME_INDEX)

  • Builds a STAR genome index from GRCh38 and GENCODE v38.
  • Executed once and reused across all samples.
  • Output: STAR_INDEX/

• Extraction of phased variants (SUBSET_1KGP_VCF)

  • Subsets 1000 Genomes phased WGS VCFs to a single sample.
  • Extracts biallelic heterozygous SNPs only.
  • Output: sampleID.1KGP.snps.het.vcf

• Alignment with WASP correction (STAR_ALIGNMENT_WASP)

  • Aligns reads using STAR with genotype-aware variant input
  • Applies WASP correction to mitigate reference mapping bias
  • Uses per-sample phased variants via --varVCFfile.
  • Output: sampleID.Aligned.sortedByCoord.out.bam (coordinate sorted BAM file with WASP tags)

• Allele-specific read counting (ALLELE_COUNT) (optional)

  • This step is only executed if a regions VCF is provided as input.
  • Uses bcftools mpileup restricted to SNPs of interest.
  • Reports:
    • Total depth (DP)
    • REF and ALT read counts (AD)
    • Computes variant allele fraction (VAF = ALT / (REF + ALT)).
  • Output: sampleID.allele_counts.with_qual_and_vaf.tsv.gz

• Summary (REPORT)

  • Generates a report using multiqc.
  • Output: multiqc_report.html

Customization

Resource requirements, containers, and execution profiles can be adjusted in:

• nextflow.config
• conf/base.config

Components

Tools:

Component	Version
STAR	2.7.11b
SAMTOOLS	1.21
BCFTOOLS	1.23
TABIX	1.11
trim-galore	2.1.0
multiqc	1.34

Container Image:

• ghcr.io/maurya-anand/nf-rna-wasp-allele-count