QTLtools-rtc(1)

Okay bestie, this is the manual page for QTLtools-rtc(1) 🧬 from the Bioinformatics tools crew, and this is Part 1 of 2 💖

✨ NAME ✨

QTLtools rtc is the Regulatory Trait Concordance score analysis tool 🦋

📜 SYNOPSIS 📜

QTLtools rtc --vcf [in.vcf|in.vcf.gz|in.bcf|in.bed.gz] --bed
quantifications.bed.gz --hotspots hotspots_b37_hg19.bed [--gwas-cis |
--gwas-trans | --mergeQTL-cis | --mergeQTL-trans] variants_external.txt
qtls_in_this_dataset.txt --out output.txt [OPTIONS]

💅 DESCRIPTION 💅

Okay so the RTC algorithm basically checks how likely it is that a GWAS variant and a molQTL share the same functional effect 👑 It does this by measuring how much the significance of the molQTL changes after you correct the molQTL phenotype for the genetic effect of the GWAS variant, then comparing that correction to what happens with every other SNP in the interval. The full method is written up at https://www.nature.com/articles/ng.3981 📖 When you're checking tissue specificity of molQTLs, we use the same method, except the GWAS variant gets treated as a molQTL from a different tissue 🌸

Here's the RTC method. Say you've got a GWAS variant sitting in the same region flanked by recombination hotspots (that's the coldspot 🧊) as a molQTL, with N variants in that coldspot:

1. Correct the phenotype for each variant in the region separately using linear regression. This gives you N pseudo-phenotypes (residuals) ✨
2. Redo the molQTL variant association with all of those pseudo-phenotypes 🔁
3. Sort the resulting p-values in decreasing order, then find the rank of the molQTL-to-GWAS-pseudo-phenotype association among all the molQTL-to-pseudo-phenotype associations 🏆
4. RTC = (N - Rank of GWAS) / N

So the RTC score ends up between 0 and 1, and higher numbers mean a more likely shared functional effect between the GWAS and molQTL variants 💜 Anything greater than or equal to 0.9 counts as a shared functional effect 🔥 If a phenotype has multiple independent molQTLs, the RTC for each one gets assessed after correcting the phenotype with all the other molQTL variants for that phenotype. That correction uses linear regression, and you take the residuals after regressing the phenotype on the other molQTLs 🫶

Now, to turn an RTC score into a probability of sharing, we run two simulations per coldspot region, called H0 and H1 🎀

The H0 scenario is when two variants in a coldspot are tagging different functional effects 🙈 For a coldspot holding colocalized GWAS and molQTL variants, we pick two random hidden causal variants. Then we find two variants (one GWAS, one molQTL) that are linked to those hidden causal variants (default r-squared ≥ 0.5). We generate a pseudo phenotype for the molQTL using the slope and intercept of the observed molQTL, plus randomly distributed residuals of the observed molQTL. After that, we rerun the RTC analysis with this new pseudo-phenotype and the GWAS and molQTL variants 🧚

The H1 scenario is when the two variants are tagging the same functional variant 💖 The setup is exactly the same as H0, except there's only one hidden causal variant, and both the GWAS and molQTL variants are randomly picked from variants linked to that same hidden causal variant 💅

---

⭐ OPTIONS ⭐

• --vcf [in.vcf|in.bcf|in.vcf.gz|in.bed.gz] — 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, like after imputation), that's what gets used as the genotype. If there's only a GT field, that gets used and converted to a dosage instead. REQUIRED. 🔥
• --bed quantifications.bed.gz — Molecular phenotype quantifications in BED format. REQUIRED.
• --out output.txt — The output file. REQUIRED.
• --hotspots recombination_hotspots.bed — Recombination hotspots in BED format. REQUIRED.
• --cov covariates.txt — Covariates you want to correct the phenotype data with 🧁
• --stats-vcf [in.vcf|in.bcf|in.vcf.gz] — Calculate D' and r-squared from this file. It defaults to the --vcf file. You need phased genotypes for D' calculations 🧬
• --stats-vcf-include-samples samples.txt — Samples to include from the --stats-vcf file. One sample ID per line 🌷
• --stats-vcf-exclude-samples samples.txt — Samples to exclude from the --stats-vcf file. One sample ID per line 🌷
• --normal — Rank normal transform the phenotype data so each phenotype comes out normally distributed. RECOMMENDED. ✨
• --conditional — molQTLs contain independent signals, so this runs the conditional analysis 🦋
• --debug — Prints debugging info to stderr. DON'T USE. 💀
• --warnings — Prints all the individual warnings you run into to stdout ⚠️
• --header — Adds a header to the output file when --chunk or --region is provided 🎀
• --individual-Dprime — Calculates D' on an individual variant basis. If you don't provide this, D' won't be calculated after the first unphased genotype is encountered 🙈
• --mem-est — Estimates memory usage and then exits 🧠
• --mem [0|1|2|3] — Decides which calculation results that might get reused stay in memory. 0 = nothing in mem, 1 = only basic, 2 = all in mem but cleaned up once they're unlikely to be reused, 3 = all in mem with no cleaning. DEFAULT=0. RECOMMENDED=2. 💜
• --window integer — Size of the cis window flanking each phenotype's start position. DEFAULT=1000000. RECOMMENDED=1000000. 🌈
• --sample integer — How many simulated RTC values to try to hit for each coldspot, so you can convert RTC into a probability. Each iteration tries to pick a unique combination of variants, so the real number of sample iterations can come in lower than this value because of how many variants are in a region. If you want to run this analysis, please provide at least 100. DEFAULT=0. 🍓
• --max-sample integer — The max number of sample iterations to try to reach --sample before quitting. Provide the actual number, not the multiplier. DEFAULT=--sample * 50. 🫶
• --R2-threshold float — The minimum r-squared required when picking a variant linked to the hidden causal variant(s) during simulations run with --sample. DEFAULT=0.5. ✨
• --D-prime-threshold float — If a pair of variants falls into different coldspots and has a D' greater than this, the RTC calculation gets extended to multiple coldspot regions including both variants. Assumes D' can be calculated. DEFAULT=OFF. NOT RECOMMENDED. 🙅
• --grp-best — Corrects for multiple phenotypes within a phenotype group 👑
• --pheno-col integer — The 1-based phenotype id column number. DEFAULT=1 or 5 when --grp-best 🧾
• --geno-col integer — The 1-based genotype id column number. DEFAULT=8 or 10 when --grp-best 🧾
• --grp-col integer — The 1-based phenotype group id column number. Only matters if --grp-best is in effect. DEFAULT=1 🧾
• --rank-col integer — The 1-based conditional analysis rank column number. Only matters if --conditional is in effect. DEFAULT=12 or 14 when --grp-best 🧾
• --best-col integer — The 1-based phenotype column number. Only matters if --conditional is in effect. DEFAULT=21 or 23 when --grp-best 🧾

✨ OUTPUT FILE ✨

--out output file 💖 okay so this is where your results land! It's a space separated file, and each row has the columns below. Heads up: columns after the 22nd only show up if you pass --sample 🙈. Also, bestie, please include chunk 0 so you get a header row. Trust me, future you will be so grateful 🫶

• Column 1: other_variant 🦋 — The variant ID from outside this dataset. It could be the GWAS variant or another molQTL.
• Column 2: our_variant 🌸 — The molQTL variant ID that lives inside this dataset.
• Column 3: phenotype — The phenotype ID.
• Column 4: phenotype_group — The phenotype group ID.
• Column 5: other_variant_chr — The external variant's chromosome.
• Column 6: other_variant_start — The external variant's start position.
• Column 7: other_variant_rank 👑 — The rank of the external variant. It only matters if the external variants are part of a conditional analysis.
• Column 8: our_variant_chr — The internal variant's chromosome.
• Column 9: our_variant_start — The internal variant's start position.
• Column 10: our_variant_rank 👑 — The rank of the internal variant. It only matters if the internal variants are part of a conditional analysis.
• Column 11: phenotype_chr — The phenotype's chromosome.
• Column 12: phenotype_start — The start position of the phenotype.
• Column 13: distance_between_variants — The distance between the two variants.
• Column 14: distance_between_other_variant_and_pheno — The distance between the external variant and the phenotype.
• Column 15: other_variant_region_index — The region index of the external variant.
• Column 16: our_variant_region_index — The region index of the internal variant.
• Column 17: region_start — The start position of the region.
• Column 18: region_end — The end position of the region.
• Column 19: variant_count_in_region — How many variants are in the region.
• Column 20: RTC ⭐ — The RTC score.
• Column 21: D' — The D' of the two variants. Only calculated if there are phased genotypes.
• Column 22: r^2 — The r squared of the two variants.
• Column 22: p_value — The p-value of the RTC score.

Yes, 22 shows up twice in the numbering. The page really said that 💀

• Column 23: unique_picks_H0 — The number of unique combinations of variants in the H0 simulations.
• Column 24: unique_picks_H1 — The number of unique combinations of variants in the H1 simulations.
• Column 25: rtc_bin_start — The lower bound of the RTC bin, based on the observed RTC score.
• Column 26: rtc_bin_end — The upper bound of the RTC bin, based on the observed RTC score.
• Column 27: rtc_bin_H0_proportion — The proportion of H0 simulated values that sit between rtc_bin_start and rtc_bin_end.
• Column 28: rtc_bin_H1_proportion — The proportion of H1 simulated values that sit between rtc_bin_start and rtc_bin_end.
• Column 29: median_r^2 — The median r-squared in the region.
• Column 30: median_H0 — The median RTC score in the H0 simulations.
• Column 31: median_H1 — The median RTC score in the H1 simulations.
• Column 32: H0 — The RTC scores observed in the H0 simulations.
• Column 33: H1 — The RTC scores observed in the H1 simulations.

---

🧚 EXAMPLES 🧚

💅 Example 1: the basic glow-up

Run RTC with GWAS variants and cis-eQTLs. This one corrects for technical covariates and rank normal transforms the phenotype ✨

QTLtools rtc --vcf genotypes.chr22.vcf.gz --bed
genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz
--hotspot hotspots_b37_hg19.bed --gwas-cis GWAS.b37.txt
permutations_all.significant.txt --normal --out rtc_results.txt

🌈 Example 2: the full main character run

Here we get RTC with GWAS variants and cis-eQTLs, plus simulations. It corrects for technical covariates, rank normal transforms the phenotype, runs a conditional analysis, and keeps the data in memory 🧠 To make parallelizing on a compute cluster easier, there's an option to run the analysis in chunks of molecular phenotypes. For example, to run chunk 12 when the example data set is split into 20 chunks, run:

QTLtools rtc --vcf genotypes.chr22.vcf.gz --bed
genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz
--hotspot hotspots_b37_hg19.bed --gwas-cis GWAS.b37.txt
conditional_all.significant.txt --normal --conditional --mem 2 --chunk
12 20 --sample 200 --out rtc_results_12_20.txt

🦋 Example 3: the whole squad on a cluster

Want to submit the whole analysis as 20 jobs on a compute cluster? Just run this. Heads up that qsub needs to be swapped for whatever job submission system you use (bsub, psub, etc...) 🫶

for j in $(seq 0 20); do
    echo "QTLtools rtc --vcf genotypes.chr22.vcf.gz --bed
    genes.50percent.chr22.bed.gz --cov genes.covariates.pc50.txt.gz
    --hotspot hotspots_b37_hg19.bed --gwas-cis GWAS.b37.txt
    conditional_all.significant.txt --normal --conditional --mem 2
    --chunk $j 20 --sample 200 --out rtc_results_$j_20.txt" | qsub
done

---

🔗 SEE ALSO 🔗

QTLtools(1) 🌸 Go check out the parent page, it's iconic.

QTLtools website: https://qtltools.github.io/qtltools 💜

🐛 BUGS 🐛

Found something weird? Please submit bugs to https://github.com/qtltools/qtltools 🙈

📚 CITATION 📚

Ongen H, Brown AA, Delaneau O, et al. Estimating the causal tissues for complex traits and diseases. Nat Genet. 2017;49(12):1676-1683. doi:10.1038/ng.3981 https://doi.org/10.1038/ng.3981 🎀

👑 AUTHORS 👑

Halit Ongen (halitongen@gmail.com), Olivier Delaneau (olivier.delaneau@gmail.com) 💖

---

Bestie, that's the whole glow-up! 🌈 This page is from QTLtools v1.3, last updated 06 May 2020, and it's still giving 💅✨