QTLtools-quan(1)

πŸŽ€ QTLtools-quan(1) πŸŽ€

Bioinformatics tools, bestie. Here's the manual for the quan mode πŸ’–

✨ NAME ✨

QTLtools quan is how you quantify gene and exon expression from RNA-seq 🧬


🧚 SYNOPSIS 🧚

QTLtools quan --bam [in.sam|in.bam|in.cram] --gtf gene_annotation.gtf
--out-prefix output [OPTIONS]

πŸ’… DESCRIPTION πŸ’…

okay so this mode is literally about counting how much each gene and exon is expressed 🌈 It uses the reads in your --bam file and checks them against the annotations in your --gtf file. The whole vibe is counting how many reads overlap the exons in that GTF.

Here's the glow-up part ✨ First, all the exons of a gene get turned into meta-exons. Any exons that overlap each other get merged into one big exon that covers all of them. Iconic πŸ‘‘

Now, what counts as a match? Any overlap at all πŸ’œ A read doesn't have to sit fully between an exon's start and end positions to count toward that exon. Just touching it is enough.

Split reads (reads that align to more than one exon) are lowkey a little different πŸ¦‹ Each exon the read overlaps gets a share based on the fraction of that read sitting on that exon. So a split read gives less than one full count to each exon it touches. Not a full count for each one, just a piece.

Reads that align to multiple exons, like overlapping exons from different genes, count toward the quantification of every exon they overlap πŸ“

Paired-end reads are their own thing πŸ™ˆ If your --bam file has paired-end reads and the two mates overlap each other (meaning the insert size is < 0), then each of those reads only contributes less than a single count toward the quantifications. The exception is if you pass --no-merge, which stops overlapping mates from being merged.

Here's a diagram of two genes with overlapping exons, plus one paired-end read where both mates are split reads that overlap each other. It shows how the quantification works πŸ’–

             x             x
            / \           / \
 +---------+   +---------+   +---------+
 | Exon1|1 |   | Exon1|2 |   | Exon1|3 |      Gene1
 +---------+   +---------+   +---------+
                           x
                          / \
            +------------+   +-------------+
            |  Exon2|1   |   |   Exon2|2   |  Gene2
            +------------+   +-------------+
                              x
                             / \
               +------------+   +----+        RNAseq Read Mate1

               |--a-||-b-||c|   |-d--||--e-|
                              x
                             / \
                     +------+   +----------+  RNAseq Read Mate2
  Left Mate1  = ((b * 0.5) + a) / (a + b + d)
  Right Mate1 = (d * 0.5)/(a + b + d)
  Left Mate2  = (b * 0.5)/(b + d + e)
  Right Mate2 = ((d * 0.5) + e)/(b + d + e)

  Exon1|2 = Left Mate1 + Left Mate2
  Exon1|3 = Right Mate1 + Right Mate2
  Exon2|1 = Left Mate1 + Left Mate2
  Exon2|2 = Right Mate2 + Right Mate2

  Gene1 = Exon1|2 + Exon1|3
  Gene2 = Exon2|1 + Exon2|2
⚠️ bestie, listen up πŸ’€ The quan mode in version 1.2 and above is not compatible with quantifications from earlier versions. That's because of bug fixes and small tweaks to how we quantify. DO NOT MIX QUANTIFICATIONS GENERATED BY EARLIER VERSIONS OF QTLTOOLS WITH QUANTIFICATIONS FROM VERSION 1.2 AND ABOVE, because it will create a bias in your dataset 😭

πŸ‘‘ OPTIONS πŸ‘‘

  • --gtf gene_annotation.gtf β€” the gene annotations in GTF format πŸ“œ You can grab them from https://www.gencodegenes.org/. REQUIRED.
  • --bam [in.bam|in.sam|in.cram] β€” your sequence data in BAM/SAM/CRAM format, sorted by chromosome and then position. Each BAM file should be one sample. REQUIRED.
  • --out-prefix output β€” the output prefix 🌸 REQUIRED.
  • --sample sample_name β€” the sample name for the BAM file. If you don't give one, the sample name becomes the BAM file path πŸŽ€
  • --rpkm β€” outputs RPKM values 🫢
  • --tpm β€” outputs TPM values 🫢
  • --xxhash β€” instead of making the unique hash for your options from the GTF file name, it uses the hash of the GTF file itself 🧚
  • --no-hash β€” skips the hash that signifies the options used in the quantification, so it doesn't go in the output file names. NOT RECOMMENDED. πŸ’€
  • --gene-types gene_type ... β€” only quantifies these gene types πŸ’– This needs a gene_type attribute in the GTF. It'll also use transcript_type if that's there.
  • --filter-mapping-quality integer β€” the minimum mapping quality a read or read pair needs before it gets considered. Set this high if you only want uniquely mapped reads πŸ¦‹ DEFAULT=10.
  • --filter-mismatch integer|float β€” the max mismatches allowed in a read. If the value is between 0 and 1, it's treated as a fraction of the read length. This needs an NM attribute in the BAM file ✨ DEFAULT=OFF.
  • --filter-mismatch-total integer|float β€” the max total mismatches allowed across a paired read. If the value is between 0 and 1, it's treated as a fraction of the combined read length. This also needs an NM attribute in the BAM file πŸ“
  • --filter-min-exon integer β€” the minimum exon length to be quantified 🌈 Exons smaller than this won't show up in the exon quantifications, but they still count toward the gene quantifications. DEFAULT=0.

✨ OUTPUT FILES ✨

Okay so bestie, unless you pass --no-hash, every output file gets a hash baked in that matches the exact combo of options you used πŸ”. Why does this matter? Because you do NOT want to merge quantifications from samples that were quantified differently. That would quietly bias your whole dataset and it is literally the worst πŸ’€ So the hash keeps everyone in their lane 🫢

.gene.count.bed .exon.count.bed .gene.rpkm.bed .exon.rpkm.bed .gene.tpm.bed .exon.tpm.bed 🌸

These are the main quantification results files, and they all have the same seven columns ✨

  • Column 1 β€” chr β€” the phenotype's chromosome πŸ—ΊοΈ
  • Column 2 β€” start β€” the phenotype's start position (0-based) πŸŽ€
  • Column 3 β€” end β€” the phenotype's end position (1-based) πŸ¦‹
  • Column 4 β€” gene|exon β€” the gene or exon ID 🧚
  • Column 5 β€” info|geneID β€” info about the gene, or the gene ID of the exon. Gene info is separated by semicolons, and the keys are L= (gene length), T= (gene type), R= (gene positions), and N= (gene name) πŸ’…
  • Column 6 β€” strand β€” the phenotype's strand 🌈
  • Column 7 β€” sample_name β€” the sample name from the BAM file πŸ‘‘

.stats πŸ“Š

This one details all the stats for the quantification. Here are the rows, in order, and yes there are seventeen of them πŸ’€

  • Row 1 β€” filtered_secondary_alignments_(does_not_count_towards_total_reads) β€” the number of secondary alignments. These do NOT count towards the total reads, so don't panic πŸ™ˆ
  • Row 2 β€” total_reads β€” the number of reads in the BAM file πŸ“š
  • Row 3 β€” filtered_unmapped β€” the number of unmapped reads πŸ‘»
  • Row 4 β€” filtered_failqc β€” the number of reads carrying the failed QC tag ❌
  • Row 5 β€” filtered_duplicate β€” the number of duplicate reads πŸ‘―
  • Row 6 β€” filtered_mapQ_less_than_X β€” the number of reads below the mapping quality threshold X πŸ“‰
  • Row 7 β€” filtered_notpaired β€” the number of pairs that were either not in the correct orientation or not properly paired πŸ™…β€β™€οΈ
  • Row 8 β€” filtered_mismatches_greater_than_X_Y β€” the number of reads that failed the mismatches-per-read filter (X) and the mismatches-total filter (Y) 🚫
  • Row 9 β€” filtered_unmatched_mate_pairs β€” the number of reads where a paired read had a missing mate πŸ’”
  • Row 10 β€” total_good β€” the number of reads that passed ALL the filters πŸ†
  • Row 11 β€” total_exonic β€” the number of reads that aligned to exons AND passed all filters 🧬
  • Row 12 β€” total_exonic_multi_counting β€” the number of reads that aligned to exons when we count reads that align to multiple exons multiple times πŸ”
  • Row 13 β€” total_merged_reads β€” the number of reads where the mate pairs were overlapping and got merged 🫠
  • Row 14 β€” total_exonic_multi_counting_after_merge_(used_for_rpkm) β€” the number of reads that aligned to exons when we merge overlapping mate pairs. This is the one used for RPKM πŸ”₯
  • Row 15 β€” good_over_total β€” the number of good reads over the total number of reads βž—
  • Row 16 β€” exonic_over_total β€” the number of exonic reads over the total number of reads βž—
  • Row 17 β€” exonic_over_good β€” the number of exonic reads over the number of good reads βž—

✨ EXAMPLE ✨

Bestie, here's a look at a real vibe. This one quantifies a sample mapped with GEM, outputs TPM and RPKM values, and takes the hash of the GTF file 🧁

QTLtools quan --bam HG00381.chr22.bam --gtf
gencode.v19.annotation.chr22.gtf.gz --out-prefix HG00381 --sample
HG00381 --rpkm --tpm --xxhash --filter-mismatch-total 8 --filter-
mapping-quality 150

✨ SEE ALSO ✨

Go check out QTLtools(1) for more ✨

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

✨ BUGS ✨

Found something broken? Slide into the DMs of the devs at <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) and Olivier Delaneau (olivier.delaneau@gmail.com) πŸ‘‘πŸ«Ά


Okay bestie that's the whole tea β˜•πŸ’– QTLtools-v1.3, last updated 06 May 2020, QTLtools-quan(1) forever. Stay iconic πŸ¦‹βœ¨