QTLtools-ase(1)

📖 QTLtools-ase(1) · Bioinformatics tools 🧬 · part 1 of 2 💖

---

✨ NAME ✨

okay bestie, the whole tool is QTLtools ase 🧬, and it's the thing that measures allele specific expression (ASE) straight from RNA-seq data 💅✨

---

🧬 SYNOPSIS 🧬

here's how you call it, no cap 💖

QTLtools ase --bam [sample.bam|sample.sam|sample.cram] --vcf
[in.vcf|in.bcf|in.vcf.gz] --sample sample_name_in_vcf --mapq  integer
--out output_file_prefix [OPTIONS]

---

💅 DESCRIPTION 💅

okay so this mode measures ASE from RNA-seq for transcribed heterozygous SNPs 🧬, following the method in <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3918453/>. the vibe in short: for each of the 12 REF/ALT pairs, it works out a reference allele mapping bias separately, as long as there are at least --sites REF/ALT sites with a minimum of --cov-bias reads overlapping 📊. pairs that don't pass that check get their bias calculated from all sites instead 🫶.

the reference allele mapping bias is basically the number of reads across lots of sites that carried the reference allele for that pair, divided by the total number of reads overlapping 💜. then --subsample decides which percentile of the highest covered sites gets subsampled 🌈. this is necessary so the bias isn't estimated mostly from super high coverage sites, but from sites covering the whole genome 🌸.

last bit: each site's ASE p-value is a two-tailed binomial test 🧚. it checks whether the observed reference allele reads are significantly different from random, given the total reads and the probability of seeing a reference allele (which is the reference allele mapping bias for that pair) ✨.

the defaults should slay for most RNAseq experiments 💅👑. but bestie, listen:

NOT advisable 🙈: lowering --baseq below 10, lowering --cov and --cov-bias below 8, or setting --subsample to 1.

NOT recommended 💀: using --keep-bans-for-bias, --keep-discordant-for-bias, --filter-duplicates, --ignore-orientation, or --legacy-options. also refrain from --auto-flip 🙅‍♀️, and instead just create correct VCF/BCF files.

we highly recommend a BCF file over a VCF for the performance boost 🔥. also use --fasta to provide the genome sequence you used, --blacklist to filter out low mappability regions 🚫, and --gtf to annotate the SNPs 🏷️. if you're working with unfiltered imputed genotypes, consider --imp-qual and --geno-prob 🎀. for troubleshooting, --filtered will list why certain variants got filtered out of the analysis 🔍. if you're I/O bound, try --group-by for better performance 🚀. and if your blacklist file has lots of overlapping or contiguous regions, --merge-on-the-fly can decrease memory usage 🧠💖.

---

🎀 OPTIONS 🎀

• --vcf, -v [in.vcf|in.bcf|in.vcf.gz] — genotypes in VCF/BCF format, sorted by chromosome and then position 🧬. it can hold multiple samples. REQUIRED, and a BCF is RECOMMENDED for performance 💖.
• --bam, -b [in.bam|in.sam|in.cram] — sequence data in BAM/SAM/CRAM format, sorted by chromosome and then position 📚. one sample per BAM file. REQUIRED.
• --ind, -i sample_name — the sample name in the VCF that matches your BAM file 🫶. REQUIRED.
• --out, -o output — the output prefix ✨. this generates output.ase and output.ref_bias files. or if you give output.gz or output.bz2, you get output.ase.gz etc. REQUIRED.
• --mapq, -q integer — the minimum mapping quality a read or read pair needs to be considered 🎯. set this to only include uniquely mapped reads. REQUIRED.
• --fasta, -f genome_sequence.fa — the genome sequence in FASTA format 🧬. make sure it's the sequence for the correct genome build. RECOMMENDED.
• --blacklist, -B poor_mappability_regions.bed — poor mappability regions in BED format 🚫. ASE estimates will be unreliable in regions with low mappability. RECOMMENDED.
• --merge-on-the-fly, -t — merges the blacklisted regions on the fly 🧠. this can reduce memory usage if there are many overlapping or contiguous regions in the blacklist.
• --gtf, -g gene_annotation.gtf — gene annotations in GTF format 🏷️. you can grab these from <https://www.gencodegenes.org/>. RECOMMENDED.
• --filtered, -l filename — the file where filtered variants get written 🔍. RECOMMENDED for troubleshooting, especially if you use --suppress-warnings.
• --reg, -r chr:start-end — the genomic region to process 🗺️. e.g. chr4:12334456-16334456, or chr5.
• --fix-chr, -F — tries to match chromosome names to the BAM file by adding or removing chr 🪄. this doesn't apply to --{include,exclude}-positions options. those should use the chromosome names from the VCF.
• --fix-id, -R — converts missing VCF variant IDs to chr_pos_refalt 💅.
• --auto-flip, -x — attempts to fix reference allele mismatches 🙈. it requires a fasta file for the reference sequence, given with --fasta. NOT RECOMMENDED.
• --print-stats, -P — prints out stats for the filtered reads at ASE sites 📊. this will be slower, and if there are sites with more reads than --max-depth, it can potentially give different results.
• --suppress-warnings, -k — suppresses the warnings about individual variants 🤫.
• --illumina13, -j — the base quality is in the Illumina-1.3+ encoding 🧪.
• --group-by, -G integer — groups variants separated by this much into batches 📦. this lets you skip streaming the whole BAM file, and may improve running time 🚀.
• --max-depth, -d integer — pileup max-depth, DEFAULT=1000000 🌊. set it to a more reasonable value if you experience memory or performance issues. it's set high by default since RNA-seq can have very high coverage sites. set to 0 for max.
• --baseq, -Q integer — minimum phred quality for a base to be considered, DEFAULT=10 💖.
• --pvalue, -p float — binomial p-value threshold for ASE output, DEFAULT=1.0 🎯.
• --cov, -c integer — minimum coverage for a genotype to be considered in ASE analysis, DEFAULT=16 📈.
• --cov-bias, -C integer — minimum coverage for a genotype to be considered in reference allele mapping bias analysis, DEFAULT=10 📉.
• --sites, -s integer — minimum number of sites to calculate a reference allele mapping bias from for a specific REF/ALT pair, DEFAULT=200 🧩. pairs with fewer than this many sites will have their bias calculated from all sites instead.
• --imp-qual-id, -I string — the VCF INFO field ID of the imputation score in the VCF, DEFAULT=INFO 🏷️.
• --geno-prob-id, -L string — the VCF FORMAT field ID of the genotype posterior probabilities for RR/RA/AA in the VCF, DEFAULT=GP 🎀.
• --imp-qual, -W float — minimum imputation score for a variant to be considered, DEFAULT=0.0 ✨.
• --geno-prob, -V float — minimum posterior probability for a genotype to be considered, DEFAULT=0.0 🌈.

📦 OUTPUT FILES ✨

Okay bestie, so QTLtools ase spits out a couple of files and I am about to walk you through every single column 🤓💖

📄 .ase

This is the main ASE results file, so it's the one you'll be staring at 👀 It has the columns below. Heads up: anything after the 23rd column only shows up if you pass --print-stats 🫶

• 1 INDIVIDUAL — the sample ID 🙋‍♀️
• 2 RSID — the SNP ID, pulled straight from the VCF file 🧬
• 3 CHR — the chromosome the SNP lives on
• 4 POS — the SNP's position
• 5 ALLELES — the alleles the SNP carries
• 6 BOTH_ALLELES_SEEN — whether or not both of the SNP's alleles actually got spotted in the RNAseq reads
• 7 MIN_ALLELE_RATIO — the minor allele ratio across the RNAseq reads
• 8 REF_COUNT — how many reference alleles show up in the reads
• 9 NONREF_COUNT — how many alternative alleles show up in the reads
• 10 TOTAL_COUNT — how many reads overlap the SNP
• 11 WEIGHTED_REF_COUNT — REF_COUNT, adjusted for reference mapping bias ⚖️
• 12 WEIGHTED_NONREF_COUNT — NONREF_COUNT, adjusted for reference mapping bias ⚖️
• 13 WRC_MINUS_WNC — WEIGHTED_REF_COUNT minus WEIGHTED_NONREF_COUNT
• 14 ALLELES_SEEN — the alleles observed in the RNAseq reads
• 15 REF_ALLELE — the reference allele
• 16 ALT_ALLELE — the alternative allele
• 17 OTHER_COUNT — how many discordant reads there are, meaning reads that are neither REF nor ALT
• 18 EXPECTED_DISCORDANT — the expected number of discordant alleles. It's calculated by adding up all the base error probabilities at the RNAseq read positions that overlap the SNP 🧮
• 19 DISCORDANT_PVAL — the p-value for the observed number of discordant reads being higher than expected. You can Bonferroni-correct these p-values and drop the significant ones from downstream analyses 🧹
• 20 REF_RATIO — reference allele mapping bias
• 21 PVALUE — the ASE p-value 🔥
• 22 CONCERN — any concerns that could make this SNP unusable get coded here. Check the OUTPUT FILE CODES section below to decode them 🚩
• 23 EXON_INFO — exons overlapping this SNP get listed here if a GTF file is provided. Each exon name is made by joining the gene ID, transcript ID, exon position, and gene name with a colon, and the names are separated by semicolons 🧩
• 24 SECONDARY — the number of secondary alignments that got filtered out
• 25 SUPPLEMENTARY — if --filter-supp is provided, the number of supplementary alignments that got filtered out
• 26 FAIL_MAPQ — the number of alignments filtered out for low mapping quality
• 27 FAILQC — if --keep-failed-qc is not provided, the number of filtered reads that failed QC according to the sequencer
• 28 DUPLICATE — if --filter-duplicates is provided, the number of filtered reads the aligner labeled as duplicates
• 29 MATE_UNMAPPED — if --keep-orphan-reads is not provided, the number of reads filtered because one mate was unmapped
• 30 WRONG_ORIENTATION — if --ignore-orientation is not provided, the number of reads filtered for being in the wrong orientation
• 31 NOT_PROPER_PAIR — if --check-proper-pairing is provided, the number of filtered reads the aligner didn't mark as properly paired
• 32 SKIPPED — the number of alignments filtered out because they had no actual bases overlapping the SNP
• 33 FAIL_BASEQ — the number of filtered reads where the base overlapping the SNP had a base quality below --baseq 📉
• 34 INDEL — if --filter-indel-reads is provided, the number of filtered reads that contained indels
• 35 DEPTH — the number of alignments overlapping the SNP position

🧪 .ref_bias

This one details the reference allele mapping bias results. It has these columns:

• 1 INDIVIDUAL — the sample ID
• 2 ALLELES — the reference/alternative allele pair
• 3 REF_ALL — total reference alleles observed across all sites
• 4 NONREF_ALL — total alternative alleles observed across all sites
• 5 SITES — number of sites for this REF/ALT pair that pass the thresholds
• 6 SUBSAMPLED_SITES — number of sites that got subsampled down to SUBSAMPLED_TO because their coverage was too high 🫧
• 7 SUBSAMPLED_TO — the coverage that SUBSAMPLED_SITES were subsampled to
• 8 PERC — reference allele mapping bias
• 9 SOURCE — whether the reference allele mapping bias was calculated from allele-specific sites or from all sites

🚫 --filtered filename

This file is the tea on every variant that got left out of the analysis and why 🍵 The first column is the coded omission reason, then comes the variant ID. The codes are all explained in the next section.

---

🗝️ OUTPUT FILE CODES

🚩 .ase file codes for the CONCERN column

Real talk: you probably want to exclude SNPs with RM, NRA, MDTA, or MDTR concerns from your analyses, because these likely have wrong genotypes 😭 SNPs with other concerns may be OK 🤷‍♀️

• RM — the reference allele in the VCF mismatches the reference sequence. Requires --fasta
• DP — multiple variants observed at the same position in the VCF file
• NRA — no reference or alternative allele observed in the RNAseq reads
• MDTA — more discordant (not reference or alternative) alleles than alternative alleles in the RNAseq reads
• MDTR — more discordant (not reference or alternative) alleles than reference alleles in the RNAseq reads
• BANS — both alleles of the SNP were not observed in the RNAseq reads
• LMAR — among the RNAseq reads, the minor allele ratio was less than 2%
• PD — pileup depth (--max-depth) was potentially exceeded around this SNP. This may stop some overlapping reads from being counted, so we recommend rerunning with a higher --max-depth 📈

🧹 --filtered file codes

• VMA — multi-allelic variant
• VU — variant position or ID is excluded by the user
• VCNIB — variant chromosome is not in the BAM file
• VB — variant is in a blacklisted region
• VI — variant is an indel
• VMRA — variant is missing either the reference or the alternative allele
• VNIF — variant is not in the reference genome. Only if --fasta was provided
• VS — reference and alternative alleles were swapped (not excluded). Only if --auto-flip and --fasta were provided 🔄
• VF — reference and alternative alleles' strand was flipped (not excluded). Only if --auto-flip and --fasta were provided
• VWR — reference allele does not match the reference sequence. Only if --auto-flip and --fasta were provided
• VBI — variant failed the imputation quality filter
• VMGT — variant with a missing GT field, or variant is not diploid
• VMG — missing genotype
• VH — homozygous variant
• VBG — variant failed the genotype probability filter
• VD — duplicate variant position
• VMI — variant with a missing ID (not excluded). It gets renamed if --fix-id is provided 🏷️
• VDK — heterozygous variant has another variant at the same position (not excluded)
• BC — variant didn't have enough coverage for reference allele mapping bias calculations
• BBANS — both alleles were not observed for reference allele mapping bias calculations
• BMDTRA — variant had more discordant alleles than reference or alternative alleles, so it was excluded from reference allele mapping bias calculations
• AC — variant didn't have enough coverage for ASE calculations
• ABANS — both alleles were not observed for ASE calculations

---

🎀 EXAMPLES

ASE analysis of a sample mapped with STAR 🌟:

QTLtools ase --vcf multi_sample.bcf --bam sample1.bam --ind sample1 --mapq 255 --out sample1 --filtered sample1.filtered.gz --gtf gencode.v19.annotation.gtf.gz --blacklist poor_mappability_regions.bed --fasta hg19.fa --fix-chr --fix-id

🔗 SEE ALSO

QTLtools(1) 💅

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

🐛 BUGS

Found a bug? Submit it here: https://github.com/qtltools/qtltools 🪲

📚 CITATION

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

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

---

And that's the whole page, bestie! ✨ This one is from QTLtools v1.3, dated 06 May 2020, and it's the QTLtools-ase(1) page. 🦋 Go forth and run some ASE analysis, you absolutely slay 💅🌸