QTLtools-cis(1)

🎀 QTLtools-cis(1) · Bioinformatics tools 🎀

✨ NAME ✨

Okay bestie, the whole name is QTLtools cis, which is just cis QTL analysis 💖

🦋 SYNOPSIS 🦋

Here's the vibe for how you call it 💅

QTLtools cis --vcf [in.vcf|in.vcf.gz|in.bcf|in.bed.gz] --bed
quantifications.bed.gz [--nominal float | --permute integer | --mapping
in.txt] --out output.txt [OPTIONS]

💜 DESCRIPTION 💜

Okay so this mode is all about finding cis (proximal) QTLs, the quantitative trait loci that affect a phenotype, and it does it with linear regression 📈. The whole method is written up here: <https://www.nature.com/articles/ncomms15452> 🌈. The flow goes like this: first we regress out the covariates you give us from the phenotype data, then we run the linear regression between those phenotype residuals and the genotype. If you pass --normal and --cov at the same time, the residuals after covariate correction get rank normal transformed ✨. It also has an efficient permutation scheme that handles the differential multiple testing burden of each phenotype, which is so iconic 🔥.

You've got three passes to pick from 🧩. There's a nominal pass (--nominal threshold) that lists all the genotype-phenotype associations below a certain threshold 📋. There's a permutation pass (--permute no_of_permutations) that empirically characterizes the null distribution of associations for each phenotype separately, which then adjusts the nominal p-value of the best association for that phenotype 🎲. And there's a conditional analysis pass (--mapping filename) that discovers multiple proximal QTLs with independent effects on a phenotype 👑.

Multiple molecular phenotypes can belong to higher order biological entities, like exons of genes, so QTLtools cis lets you group phenotypes to maximize the discoveries in those particular cases 🫶. You can either aggregate the phenotypes in a group into a single phenotype via PCA (--grp-pca1) or by taking their mean (--grp-mean). Or you can directly use all the individual phenotypes in an extended permutation scheme that accounts for their number and correlation structure (--grp-best). In their experience, --grp-best outperforms the other options for expression QTLs (eQTLs) 🧬.

The conditional analysis pass is a whole moment 💅. It first uses permutations to derive a nominal p-value threshold per phenotype that varies and reflects the number of independent tests per cis-window. Then it uses a forward-backward stepwise regression to learn the number of independent signals per phenotype, determine the best candidate variant per signal, and assign all significant hits to the independent signal they relate to 🌸.

Real talk: linear regressions assume normally distributed data, so we highly recommend the --normal option. It rank normal transforms the phenotype quantifications so you don't get false positive associations from outliers 🙈 Trust me on this one, periodt 💖.

🎀 OPTIONS 🎀

Okay, here are the options, bestie 🧚 Required ones are flagged, and I kept every default and recommendation.

• --vcf [in.vcf|in.bcf|in.vcf.gz|in.bed.gz] — Your genotypes in VCF/BCF format, or another molecular phenotype in BED format. If a variant's genotype FORMAT has a DS field (the dosage calculated from genotype probabilities, e.g. after imputation), that's what gets used as the genotype. If there's only the GT field in the genotype FORMAT, then that's used and it gets converted to a dosage. REQUIRED.
• --bed quantifications.bed.gz — Your molecular phenotype quantifications in BED format. REQUIRED.
• --out output.txt — The output file. REQUIRED.
• --cov covariates.txt — The covariates you want to correct the phenotype data with 🧼
• --normal — Rank normal transforms the phenotype data so each phenotype is normally distributed. RECOMMENDED.
• --std-err — Calculates and outputs the standard error of the beta (slope) 📏
• --window integer — The size of the cis window flanking each phenotype's start position. DEFAULT=1000000. RECOMMENDED=1000000.
• --nominal float — Calculates the nominal p-value for the genotype-phenotype associations and prints out the ones that pass the threshold you give. Give 1.0 to print everything. Mutually exclusive with --permute and --mapping.
• --permute integer — Adjusts the best nominal p-value for this phenotype, accounting for the number of variants and the linkage disequilibrium in its cis-window. We recommend at least 1000 permutations for the final analysis, and in most cases you'll see diminishing returns when going over 5000. But if you're doing exploratory stuff like figuring out which/how many covariates to include, you can go as low as 100. Mutually exclusive with --nominal and --mapping. RECOMMENDED=1000.
• --mapping thresholds_filename — The conditional analysis 🧩 First you need to run a permutation analysis, then generate a thresholds file using the runFDR_cis.R script in the script directory. Mutually exclusive with --nominal and --permute.
• --grp-best — Corrects for multiple phenotypes within a phenotype group. Mutually exclusive with --grp-pca1 and --grp-mean.
• --grp-pca1 — Runs PCA on the phenotypes within a phenotype group and uses PC1 for association testing. Mutually exclusive with --grp-best and --grp-mean.
• --grp-mean — Averages the phenotypes within a group and uses the results for association testing. Mutually exclusive with --grp-best and --grp-pca1.
• --chunk integer1 integer2 — For parallelization 🚀 Divides the data into integer2 number of chunks and processes chunk number integer1. Chunk 0 will print a header. Mutually exclusive with --region and --region-pair. The minimum number of chunks has to be at least the same number of chromosomes in the --bed file.
• --region chr:start-end — The genomic region to be processed 🗺️ E.g. chr4:12334456-16334456, or chr5. Mutually exclusive with --chunk and --region-pair.
• --region-pair chr:start-end chr:start-end — The genomic region for genotypes, followed by the region for phenotypes, to be processed. Mutually exclusive with --chunk and --region.

The --permute output file is where all the juicy results live 💖 It's space separated, and bestie, which columns show up depends on the options you pick 🎀 We really recommend including chunk 0 so you get a header and nobody gets confused 🫶

• Column 1: phe_id | grp_id — the phenotype ID, or the phenotype group ID if one of the grouping options is provided 🌸
• Column 2: phe_chr — the phenotype's chromosome 🧚
• Column 3: phe_from — where the phenotype starts 🦋
• Column 4: phe_to — where the phenotype ends ✨
• Column 5: phe_strd — the phenotype's strand 🎀
  • Column 5.1: phe_id | ve_by_pc1 | n_phe_in_grp — only printed with --group-best | --group-pca1 | --group-mean. It's the phenotype ID for --group-best, the variance explained by PC1 for --group-pca1, or the number of phenotypes in the group for --group-mean 💅
  • Column 5.2: n_phe_in_grp — only printed with --group-pca1 | --group-mean. The number of phenotypes in the group 🌈
• Column 6: n_var_in_cis — the number of variants in the cis window for this phenotype 🍓
• Column 7: dist_phe_var — the distance between the variant and the phenotype start positions 💜
• Column 8: var_id — the ID of the most significant variant 👑
• Column 9: var_chr — the most significant variant's chromosome 🔥
• Column 10: var_from — the start position of the most significant variant ⭐
• Column 11: var_to — the end position of the most significant variant ✨
• Column 12: dof1 — the number of degrees of freedom used to compute the p-values 🦋
• Column 13: dof2 — the estimated number of degrees of freedom used in the beta approximation p-value calculations 🌸
• Column 14: bml1 — the first shape parameter of the fitted beta distribution (alpha). These should be close to 1 💖
• Column 15: bml2 — the second shape parameter of the fitted beta distribution (beta). This one corresponds to the effective number of independent tests in the region 🧚
• Column 16: nom_pval — the nominal p-value of the association between the most significant variant and the phenotype 🎀
• Column 17: r_squared — the r squared of the linear regression 💅
• Column 18: slope — the beta (slope) of the linear regression 🌈
  • Column 18.1: slope_se — the standard error of the beta. Only printed if --std-err is provided 🫶
• Column 19: adj_emp_pval — the adjusted empirical p-value from permutations. This is the adjusted p-value not using the beta approximation, and it's literally calculated as:

(number of p-values observed during permutations that were smaller than or equal to the nominal p-value + 1) / (number of permutations + 1)

The most significant p-value you could possibly get is 1 / (number of permutations + 1) 😭

• Column 20: adj_beta_pval — the adjusted empirical p-value given by the fitted beta distribution 👑

We strongly recommend using this adjusted p-value in any downstream analysis. Seriously, bestie, use this one 💀

🗺️ --mapping output file

okay so this is the mapping file 🗺️ it's space-separated, and some columns only show up depending on your options 💅 we recommend including chunk 0 so you get a header and nobody gets confused 🎀

• Column 1: phe_id | grp_id: the phenotype ID, or the phenotype group ID if you used one of the grouping options 🧚
• Column 2: phe_chr: the phenotype's chromosome 🌸
• Column 3: phe_from: where the phenotype starts 🦋
• Column 4: phe_to: where the phenotype ends 💜
• Column 5: phe_strd: the phenotype's strand 🌈
• Column 5.1: phe_id | ve_by_pc1 | n_phe_in_grp: only shows up with --group-best | --group-pca1 | --group-mean. It's the phenotype ID, the variance explained by PC1, or the number of phenotypes in the group, for --group-best, --group-pca1, and --group-mean respectively ✨
• Column 5.2: n_phe_in_grp: only shows up with --group-pca1 | --group-mean. It's the number of phenotypes in the phenotype group 🫶
• Column 6: n_var_in_cis: how many variants are in the cis window for this phenotype 🍓
• Column 7: dist_phe_var: the distance between the variant and the phenotype start positions 📏
• Column 8: var_id: the ID of the most significant variant 👑
• Column 9: var_chr: the most significant variant's chromosome 🌸
• Column 10: var_from: the start position of the most significant variant 🦋
• Column 11: var_to: the end position of the most significant variant 💜
• Column 12: rank: the rank of the association. This tells you whether the variant got mapped to the best signal (rank=0), the second best (rank=1), and so on 🔥 Because of that, the highest rank for a given phenotype tells you how many independent signals there are (e.g. rank=2 means 3 independent signals) 🧠
• Column 13: fwd_pval: the nominal forward p-value of the association between the most significant variant and the phenotype 💅
• Column 14: fwd_r_squared: the r squared of the forward linear regression 📈
• Column 15: fwd_slope: the beta (slope) of the forward linear regression 📈
• Column 15.1: fwd_slope_se: the standard error of the forward beta. Only printed if --std-err is provided 🎀
• Column 16: fwd_best_hit: whether or not this variant was the forward most significant variant 👑
• Column 17: fwd_sig: whether this variant was significant. Currently all variants are significant, so this one is redundant 💀
• Column 18: bwd_pval: the nominal backward p-value of the association between the most significant variant and the phenotype 💅
• Column 19: bwd_r_squared: the r squared of the backward linear regression 📈
• Column 20: bwd_slope: the beta (slope) of the backward linear regression 📈
• Column 20.1: bwd_slope_se: the standard error of the backward beta. Only printed if --std-err is provided 🎀
• Column 21: bwd_best_hit: whether or not this variant was the backward most significant variant 👑
• Column 22: bwd_sig: whether this variant was significant. Currently all variants are significant, so this one is redundant 💀

---

✨ NOMINAL ANALYSIS EXAMPLES ✨

• Nominal pass, correcting for technical covariates, rank normal transforming the phenotype, printing out associations with a p-value <= 0.01 on chromosome 22 between 17000000 and 18000000 bp, and excluding some samples (see QTLtools (1)) 🌸:

QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --nominal 0.01 --region chr22:17000000-18000000 --normal --out nominals.txt --exclude-samples sample_names_to_exclude.txt

• Nominal pass with parallelization, correcting for technical covariates, rank normal transforming the phenotype, and printing out associations with a p-value <= 0.01. To make parallelization on a compute cluster easier, we built an option to run the analysis in chunks of molecular phenotypes 🧚 For instance, to run the analysis on chunk 12 when splitting the example data set into 20 chunks, run:

QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --nominal 0.01 --chunk 12 20 --normal --out nominals_12_20.txt

• If you want to submit the whole analysis as 20 jobs on a compute cluster, just run this (qsub needs to be swapped out for the job submission system you use [bsub, psub, etc...]) 🔥:

for j in $(seq 0 20); do
    echo "QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --nominal 0.01 --chunk $j 20 --normal --out nominals_$j\_20.txt" | qsub
done

---

✨ PERMUTATION ANALYSIS EXAMPLES ✨

• Permutation pass, correcting for technical covariates, rank normal transforming the phenotype, and running 1000 permutations with a specific random seed on chromosome 22 between 17000000 and 18000000 bp 🎀:

QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --permute 1000 --region chr22:17000000-18000000 --normal --seed 1354145 --out permutation.txt

• Permutation pass with parallelization, correcting for technical covariates, rank normal transforming the phenotype, and running 5000 permutations. To make parallelization on a compute cluster easier, we built an option to run the analysis in chunks of molecular phenotypes 🧚 For instance, to run the analysis on chunk 12 when splitting the example data set into 20 chunks, run:

QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --permute 5000 --chunk 12 20 --normal --out permutations_12_20.txt

• If you want to submit the whole analysis as 20 jobs on a compute cluster, just run this (qsub needs to be swapped out for the job submission system you use [bsub, psub, etc...]) 🔥:

for j in $(seq 0 20); do
    echo "QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --permute 5000 --chunk $j 20 --normal --out permutations_$j\_20.txt" | qsub
done

• When your phenotypes in the BED file are grouped, you can run a permutation pass at the phenotype group level to discover group-level QTLs 👑:

QTLtools cis --vcf genotypes.chr22.vcf.gz --bed exons.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --permute 1000 --normal --grp-best --out permutation.group.txt

---

✨ CONDITIONAL ANALYSIS EXAMPLE ✨

Conditional analysis to discover independent signals 🦋

1. First we need to run a permutation analysis (see previous section), then calculate the nominal p-value threshold for each gene. Here an FDR of 5% is given as an example 💅:

cat permutations_*.txt | gzip -c > permutations_all.txt.gz
Rscript ./script/qtltools_runFDR_cis.R permutations_all.txt.gz 0.05 permutations_all

2. Now you can proceed with the actual conditional analysis. Here we split into 20 chunks, and when they're all complete, concatenate the results 🫶:

for j in $(seq 0 20); do
    echo "QTLtools cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz --mapping permutations_all.thresholds.txt --chunk $j 20 --normal --out conditional_$j\_20.txt" | qsub
done
cat conditional_*.txt > conditional_all.txt

3. If you're interested in the most significant variants per independent signal, you can filter the results using the backward p-value 👑:

awk '{ if ($21 == 1) print $0 }' conditional_all.txt > conditional_top_variants.txt

---

✨ SEE ALSO ✨

• QTLtools(1) 🌈
• QTLtools website: https://qtltools.github.io/qtltools 💖

✨ BUGS ✨

• Heads up, bestie 😭 Versions up to and including 1.2 have a bug when reading missing genotypes in VCF/BCF files. It hits variants with a DS field in their genotype's FORMAT and a missing genotype (DS field is .) in one of the samples. When that happens, genotypes for all the samples get set to missing, which effectively removes that variant from the analyses. Not it 💀

Found a bug that isn't on this list? Report it at https://github.com/qtltools/qtltools 🫶

---

🎀 CITATION 🎀

If you cite QTLtools, this is the one 📚✨

Delaneau, O., Ongen, H., Brown, A. et al. A complete tool set for
molecular QTL discovery and analysis. Nat Commun 8, 15452 (2017).
<https://doi.org/10.1038/ncomms15452>

---

👑 AUTHORS 👑

The girlies behind this tool are Olivier Delaneau (olivier.delaneau@gmail.com) and Halit Ongen (halitongen@gmail.com) 💜💅

That's the whole page, bestie! 🦋 It came from QTLtools-v1.3, dated 06 May 2020, and lives under QTLtools-cis(1) 🌈✨