sqtlseeker2-nf

A pipeline for splicing quantitative trait loci (sQTL) mapping.

The pipeline performs the following analysis steps:

Index the genotype file
Preprocess the transcript expression data
Test for association between splicing ratios and genetic variants in cis (nominal pass)
Obtain an empirical P-value for each phenotype (permutation pass, optional)
Control for multiple testing

For details on each step, please read sQTLseekeR2 documentation.

The pipeline uses Nextflow as the execution backend. Please check Nextflow documentation for more information.

Requirements

Unix-like operating system (Linux, MacOS, etc.)
Java 8 or later
Docker (v1.10.0 or later) or Singularity (v2.5.0 or later)

Quickstart (~2 min)

Install Nextflow:
```
curl -fsSL get.nextflow.io | bash
```
Make a test run:
```
./nextflow run guigolab/sqtlseeker2-nf -with-docker
```
Note: set -with-singularity to use Singularity instead of Docker.

Important: Since release 22.12.0-edge, DSL1 is not further supported in Nextflow. Until sqtlseeker2-nf is migrated to DSL2, the pipeline should be run using an older Nextflow release. This can be done using NXF_VER before Nextflow commands, e.g. NXF_VER=22.04.0 ./nextflow run guigolab/sqtlseeker2-nf -with-docker.

Pipeline usage

Launching the pipeline with the --help parameter shows the help message:

nextflow run sqtlseeker2-nf --help

N E X T F L O W  ~  version 0.27.2
Launching `sqtlseeker2.nf` [admiring_lichterman] - revision: 28c86caf1c

sqtlseeker2-nf ~ A pipeline for splicing QTL mapping
----------------------------------------------------
Run sQTLseekeR2 on a set of data.

Usage: 
    sqtlseeker2-nf [options]

Options:
--genotype GENOTYPE_FILE    the genotype file
--trexp EXPRESSION_FILE     the transcript expression file
--metadata METADATA_FILE    the metadata file
--genes GENES_FILE          the gene location file
--dir DIRECTORY             the output directory
--mode MODE                 the run mode: nominal or permuted (default: nominal)
--win WINDOW                the cis window in bp (default: 5000)
--covariates COVARIATES     include covariates in the model (default: false)
--fdr FDR                   false discovery rate level (default: 0.05)
--min_md MIN_MD             minimum effect size reported (default: 0.05)
--svqtl SVQTLS              report svQTLs (default: false)

Additional parameters for mode = nominal:
--ld LD                     threshold for LD-based variant clustering (default: 0, no clustering)
--kn KN                     number of genes per batch in nominal pass (default: 10)

Additional parameters for mode = permuted:
--kp KP                     number of genes per batch in permuted pass (default: 10)
--max_perm MAX_PERM         maximum number of permutations (default: 1000)

Input files and format

sqtlseeker2-nf takes as input files the following:

Genotype file. Contains the genotype of each sample, coded as follows: 0 for REF/REF, 1 for REF/ALT, 2 for ALT/ALT, -1 for missing value. The first four columns should be: chr, start, end and snpId. This file needs to be sorted by coordinate.
Transcript expression file. Contains the expression of each transcript in each sample (e.g. read counts, RPKM, TPM). It is not recommended to use transformed (log, quantile, or any non-linear transformation) expression. Columns trId and geneId, corresponding to the transcript and gene IDs, are required.
Metadata file. Contains the covariate information for each sample. In addition, it defines the groups or conditions for which sQTL mapping will be performed. The first columns should be: indId, sampleId, group, followed by the covariates. This file defines which samples will be tested.
Gene location file. Contains the location of each gene. Columns chr, start, end and geneId are required. This file defines which genes will be tested.

Example data is available for the test run.

Pipeline results

sQTL mapping results are saved into the folder specified with the --dir parameter. By default it is the result folder within the current working directory.

Output files are organinzed into subfolders corresponding to the different groups specified in the metadata file:

result
└── groups
    ├── group1                            
    │   ├── all-tests.nominal.tsv          
    │   ├── all-tests.permuted.tsv         
    │   ├── sqtls-${level}fdr.nominal.tsv      
    │   └── sqtls-${level}fdr.permuted.tsv     
    ├── group2
   ...

Note: if only a nominal pass was run, files *.permuted.tsv will not be present.

Output files contain the following information:

all-tests.nominal.tsv

geneId: gene name
snpId: variant name
F: test statistic
nb.groups: number of genotype groups
md: maximum difference in relative expression between genotype groups (sQTL effect size)
tr.first/tr.second: the transcript IDs of the two transcripts that change the most, in opposite directions
info: number of individuals in each genotype group, including missing values (-1,0,1,2)
pv: nominal P-value

if --svqtl true

F.svQTL: svQTL test statistic
nb.perms.svQTL: number of permutations for svQTL test
pv.svQTL: svQTL nominal P-value

if --ld ${r2}

LD: other variants in linkage disequilibrium with snpId above a given r² threshold > 0

sqtls-${level}fdr.nominal.tsv (in addition to the previous)

fdr: false discovery rate (computed across all nominal tests)
fdr.svQTL: svQTL FDR

all-tests.permuted.tsv

geneId: gene name
variants.cis: number of variants tested in cis
LD: median linkage disequilibrium in the region (r²)
best.snp: ID of the top variant
best.nominal.pv: P-value of the top variant
shape1: first parameter value of the fitted beta distribution
shape2: second parameter value of the fitted beta distribution (effective number of independent tests in the region)
nb.perm: number of permutations
pv.emp.perm: empirical P-value, computed based on permutations
pv.emp.beta: empirical P-value, computed based on the fitted beta distribution
runtime: run time in minutes

sqtls-${level}fdr.nominal.tsv (in addition to the previous)

fdr: false discovery rate (computed across empirical P-values)
p_tn: gene-level threshold for nominal P-values

Cite sqtlseeker2-nf

If you find sqtlseeker2-nf useful in your research please cite the related publication:

Garrido-Martín, D., Borsari, B., Calvo, M., Reverter, F., Guigó, R. Identification and analysis of splicing quantitative trait loci across multiple tissues in the human genome. Nat Commun 12, 727 (2021). https://doi.org/10.1038/s41467-020-20578-2