BAMFEATURECOUNT(1)

💖 NAME 💖

bamfeaturecount — she evaluates the alignments your RNA-seq aligner produced 🧬✨

(it's a page from the General Commands Manual, section 1, for the record)

🌸 SYNOPSIS 🌸

bamfeaturecount annotation.gtf mapped.bam [options]

🦋 DESCRIPTION 🦋

okay so bestie, bamfeaturecount takes the alignments that an RNA-seq aligner made (think STAR, like, that whole vibe) and spits out the average coverage it found for every single transcript and every single gene. slay. 💅

To do that it needs two things: an annotation file in GTF format, plus a BAM file stuffed with RNA-seq alignments. and here's the tea 🫖 — that annotation file has to have been run through the filtergtf program first, because that's what guarantees the sorting is compatible with what bamfeaturecount expects. no filtergtf, no vibes, no deal.

The output is giving two flavors of line: one starting with [transcript] and another starting with [gene]. Both of them hand you the info as a tab separated</> set of columns.

Here's what a transcript line actually looks like, fr fr:

[transcript]    chr1:C1orf159   ENST00000472741.5
(0.0618557,5.38144,97)  (0.0496454,11.305,564)  (0.0514372,10.4357,661)
[1116059,1116087)       [1091990,1092103)       [1091045,1091565)

transcript lines always have at least 6 columns, no exceptions. let me break it down because it's a lot:

• column 2 — the reference sequence name (exactly as given in the BAM header and the GTF file) and the name of the gene concerned, glued together with a colon. in the example above that's chr1 and C1orf159.
• column 3 — the transcript identifier (that's transcript_id in the GTF file).
• columns 4, 5 and 6 — each one holds either a triplet of numbers (A,B,C) or the symbol *.
  • column 4 — regions unique to this transcript, i.e. stretches on the genome not shared with any other transcript.
  • column 5 — regions shared by at least one other transcript.
  • column 6 — the info for all regions covered by the transcript (unique and non unique together).
• if a transcript doesn't have any respective intervals (like, if every base is also covered by at least one other transcript, or all the bases are unique to this transcript) then that column just contains the symbol *. otherwise you get a triplet (A,B,C) where:
  • A is the fraction of bases not covered by any alignment (in the example, 0.0618557, so 6% of the bases unique to this transcript are not covered)
  • B is the average coverage (in the example the average sequencing depth on the unique bases for this transcript is 5.38)
  • C is the total number of bases (in the example, 97 bases of this transcript are not shared with any other transcript)
• columns beyond 6 (> 6) — these hold the zero bases intervals of exons for this transcript on the reference sequence.

⚠️ the gene lines? disregard them in the current version of the program. like, they're there, but they're not the moment. 🙈

🔑 the key=value pairs you can pass 🔑

these can all be given as key=value:

• T=<filename> — set the prefix for temporary file names ✨
• verbose=<1> — print some progress report while processing
• threads=<1> — number of threads used for processing. set this to 0 to use all cores detected on the machine 👑
• mapqmin=<255> — minimum mapping quality allowed for alignments considered. by default this is 255 (the value used by STAR to mark uniquely mapped reads).
• mapqmax=<255> — maximum mapping quality allowed for alignments considered. by default this is 255 (the value used by STAR to mark uniquely mapped reads).
• uncoveredthres=<0.1> — maximum fraction of bases allowed to be uncovered in a transcript so the transcript will be reported (i.e. the minimum value allowed for the first number given in column 6 of the transcript output).
• uniqueuncoveredthres=<0.1> — maximum fraction of bases allowed to be uncovered in the unique region of a transcript so the transcript will be reported (i.e. the minimum value allowed for the first number given in column 4 of the transcript output).
• exclude=<SECONDARY> — do not include reads in the output that have any of the given flags set. the flags are given separated by commas. valid flags are:
  • PAIRED — read was paired in sequencing
  • PROPER_PAIR — read has been mapped as part of a proper pair
  • UNMAP — read was not mapped
  • MUNMAP — mate of read was not mapped
  • REVERSE — read was mapped to the reverse strand
  • MREVERSE — mate of read was mapped to the reverse strand
  • READ1 — read was first read of a pair during sequencing
  • READ2 — read was second read of a pair during sequencing
  • SECONDARY — alignment is secondary, i.e. an alternative mapping to the primary alignment in the same file
  • QCFAIL — read as marked as having failed quality control
  • DUP — read is marked as a duplicate of another read in the same file (see bammarkduplicates)
  • SUPPLEMENTARY — read is marked as supplementary alignment
• exportcdna=<0> — instead of feature counting, generate a FastA file containing the CDNA as designated by the GTF annotation file. the second parameter (the BAM file) needs to be specified, but it will not be read. this option requires a reference FastA file suitable for the GTF file to be provided via the reference key.
• reference=<> — name of a reference FastA file. this is required for the exportcdna option.

🎀 AUTHOR 🎀

Written by German Tischler-Höhle. 💜

🐛 REPORTING BUGS 🐛

Report bugs to <germant@miltenyibiotec.de> — bestie will hear you out, fr fr.

©️ COPYRIGHT ©️

Copyright © 2009-2019 German Tischler-Höhle, © 2011-2014 Genome Research Limited. License GPLv3+: GNU GPL version 3 <http://gnu.org/licenses/gpl.html>

This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law.

that's the whole thing, no notes. this page came to you from BIOBAMBAM, dated August 2019 — and bamfeaturecount(1) ate and left no crumbs 🌈⭐💖