featureCounts

Prérequis

Cette section suppose les notions de la section 6.1 sur la librairie RNA-seq : ARN messager (ARNm), orientation des transcrits et appariement des lectures. Elle s'appuie sur la section 3.3 pour le format SAM/BAM, les drapeaux et les coordonnées. Le module 6 fournit le fil conducteur transcriptomique ; les comptages produits ici alimenteront les sections de normalisation et de test différentiel.

Introduction

Une expérience de RNA-seq produit des millions de lectures courtes, fragments d'ARN convertis en ADN complémentaire puis séquencés. Après alignement, la question devient concrète : combien de lectures proviennent de chaque gène dans chaque échantillon ? Ce nombre est la matière première de l'analyse différentielle et de l'interprétation fonctionnelle. Pourtant, passer d'un BAM à une matrice de comptages n'est pas un simple comptage. Il faut décider ce qu'est un gène, quels exons lui appartiennent, comment traiter une lecture à cheval sur deux exons, un fragment dont les deux extrémités tombent dans des régions différentes, ou une lecture qui s'aligne aussi bien sur plusieurs loci.

Le défi est double. L'annotation doit être lue correctement : un GTF et un GFF3 ne décrivent pas les mêmes relations, et les identifiants de gènes ne sont pas toujours ceux attendus. Ensuite, chaque lecture doit être confrontée à des milliers de features avec des règles explicites pour les cas ambigus. Un compteur qui additionne les chevauchements exon par exon double les lectures aux jonctions. Un compteur qui ignore le brin peut produire des résultats plausibles mais inversés. Un compteur qui accepte tous les alignements multiples gonfle les familles de gènes.

featureCounts, distribué avec Subread, associe les exons d'un même gène à une méta-feature, puis compte les lectures ou fragments qui la chevauchent. Il produit une matrice de comptages bruts et un fichier de statistiques expliquant les lectures non assignées. Cette section part de l'annotation, construit un comptage à la main, puis reproduit la démarche avec featureCounts sur un BAM synthétique.

Théorie

L'annotation décrit des features

Un GTF (Gene Transfer Format) est un tableau à neuf colonnes séparées par des tabulations : chromosome, source, type de feature, début, fin, score, brin, cadre et attributs. Les coordonnées sont 1-based et inclusives : un exon de 1 à 150 contient les bases 1 et 150. Le brin vaut +, - ou . s'il est inconnu. Les attributs sont des paires clé-valeur terminées par un point-virgule, par exemple gene_id "geneA"; transcript_id "txA";.

Un GFF3 (General Feature Format version 3) utilise aussi neuf colonnes et des coordonnées 1-based inclusives, mais ses attributs suivent la syntaxe clé=valeur et sa hiérarchie repose sur ID et Parent. Un exon n'y porte pas nécessairement gene_id : il pointe vers un transcrit par Parent, lequel pointe vers un gène. featureCounts attend un GTF ou un fichier SAF (simplified annotation format) ; un GFF3 doit être converti ou enrichi avec les attributs attendus avant d'utiliser -g gene_id.

Trois niveaux coexistent. Le gène est l'unité fonctionnelle et la méta-feature qui nous intéresse. Le transcrit est une isoforme produite par épissage ; un gène peut en avoir plusieurs. L'exon est un segment présent dans un transcrit, souvent répété une fois par isoforme. Avec -t exon -g gene_id, featureCounts sélectionne les lignes exon et regroupe celles qui partagent le même gene_id. Les lignes gene et transcript ne sont pas comptées directement.

Le choix de -t et -g fixe le niveau du comptage. -t exon -g gene_id produit des comptages par gène, adaptés à l'analyse différentielle classique. -t exon -g transcript_id produit des comptages par transcrit, plus difficiles à interpréter à cause des isoformes partagées. -t gene compte le corps entier, introns compris : ce n'est généralement pas l'expression mature, mais cela peut servir pour la transcription naissante. Il faut expliciter le niveau avant de comparer des échantillons.

Une lecture, un fragment, une méta-feature

featureCounts distingue la lecture et le fragment. En single-end, chaque lecture est une observation. En paired-end, les deux mates proviennent du même fragment ; les compter séparément revient à compter deux fois le même fragment. La règle de chevauchement est généreuse : un fragment chevauche une feature si au moins une base de l'un de ses mates la chevauche.

Une lecture est assignée lorsqu'elle chevauche exactement une méta-feature compatible. Elle est ambiguë lorsque son alignement chevauche plusieurs gènes, et multi-alignée lorsque la même lecture possède plusieurs alignements possibles. Une lecture ambiguë a un seul alignement qui traverse deux gènes chevauchants ; une lecture multi-alignée a plusieurs alignements, par exemple dans une famille de gènes ou une région dupliquée. featureCounts les traite avec des options différentes : -O pour les chevauchements entre gènes, -M pour les alignements multiples.

Le regroupement des exons est crucial. Si un gène possède les exons 1–150 et 200–300, une lecture de 140 à 210 chevauche les deux. Additionner les comptages exon par exon donnerait deux lectures pour une seule observation. featureCounts construit une méta-feature par gene_id et compte la lecture une seule fois si elle ne chevauche aucun autre gène. C'est l'union des exons, par opposition à la somme des comptages par exon. La colonne Length correspond à la longueur cumulée des exons, mais elle ne corrige pas le comptage : featureCounts renvoie des comptages bruts.

Formalisation et coût

Notons $R_j$ les lectures ou fragments de l'échantillon $j$, $G$ les gènes et $a_{rg}$ l'indicatrice valant 1 si la lecture $r$ est assignée au gène $g$. Le comptage s'écrit

où $w_r = 1$ sans fractionnement et $w_r = 1/K_r$ avec --fraction, $K_r$ étant le nombre de gènes compatibles. Sans -O ni -M, une lecture ambiguë ou multi-alignée a tous ses $a_{rg}$ nuls. Avec ces options, elle peut recevoir plusieurs $a_{rg}=1$ ; sans --fraction, la somme des contributions d'une lecture peut dépasser 1. Le total des comptages d'un échantillon peut donc être supérieur au nombre de lectures assignées.

Brin, protocole et paires

-s contrôle la stratégie de brin. 0 ignore le brin. 1 correspond à un protocole orienté dans le sens du transcrit (forward) : la lecture doit avoir le même brin que le gène. 2 correspond à un protocole inverse (reverse, souvent dUTP) : la lecture doit avoir le brin opposé. Le choix dépend de la chimie de la librairie, pas du BAM. Une stratégie erronée produit des comptages très faibles ou nuls et un enrichissement apparent en transcrits antisens.

Pour les données appariées, -p active le mode paired-end. Depuis Subread 2.0.2, cette option indique seulement que les données sont appariées ; --countReadPairs demande de compter chaque paire une fois. Avec Subread 2.0.6, la commande est donc -p --countReadPairs. Sans cette option, chaque mate est compté séparément et les comptages doublent approximativement. Dans les versions antérieures à 2.0.2, --countReadPairs n'existe pas et -p suffit : un script portable doit vérifier la version.

Deux filtres complètent le mode apparié. -B exige que les deux extrémités soient alignées ; -C exclut les paires dont les mates sont sur des chromosomes ou des brins différents, souvent chimériques. -d et -D encadrent la longueur du fragment, avec des valeurs par défaut de 50 et 600 pb en mode apparié. Un fragment hors intervalle est écarté et apparaît dans les statistiques. Ces filtres changent le dénominateur et doivent être rapportés.

Multi-mapping, chevauchements et longueur

Une lecture multi-alignée provient d'une région répétée ou d'une famille de gènes proches. Par défaut, featureCounts ne la compte pas : elle est signalée comme Unassigned_MultiMapping. -M l'autorise ; si ses alignements touchent plusieurs gènes, -O est nécessaire pour l'assigner à chacun. Sans --fraction, chaque gène reçoit le compte entier ; avec --fraction, le compte est partagé entre les gènes compatibles. Pour une famille très similaire, le fractionnement est souvent plus prudent, mais il produit des comptages non entiers qui doivent être assumés dans l'analyse.

Une lecture peut aussi chevaucher plusieurs gènes avec un seul alignement, notamment dans les génomes compacts ou les annotations imbriquées. Par défaut, featureCounts la déclare ambiguë et ne la compte pas. -O l'assigne à tous les gènes chevauchés, éventuellement avec --fraction. --largestOverlap choisit le gène dont le recouvrement est le plus grand. Aucune option n'est neutre : -O augmente le total, --largestOverlap impose une règle arbitraire, et l'absence d'option jette de l'information. Le fichier .summary mesure l'importance de chaque catégorie.

La longueur intervient à deux niveaux. La longueur du gène, obtenue en cumulant ses exons, influence le nombre de lectures attendu : un gène long capture plus de fragments à expression égale. Les comptages bruts ne corrigent pas cet effet ; TPM et FPKM le font, tandis que DESeq2 modélise la profondeur et la dispersion. La longueur du fragment est contrôlée par -d et -D en paired-end. Il ne faut pas confondre une lecture courte de 50 pb et un fragment de 300 pb.

Exemple travaillé à la main

L'annotation datasets/rnaseq/annotation.gtf contient deux gènes sur chr1. geneA, sur le brin +, possède les exons 1–150 et 200–300 ; geneB, sur le brin -, possède l'exon 400–700. La longueur cumulée de geneA est $150 + 101 = 251$ pb ; celle de geneB est 301 pb. Les lignes transcript et CDS sont ignorées par -t exon.

Quatre lectures suffisent. r1 couvre 10–25 sur le brin + et chevauche le premier exon de geneA. r2 couvre 140–210 sur le brin + : elle chevauche la fin du premier exon et le début du second, mais elle est comptée une seule fois pour geneA. r3 couvre 420–430 sur le brin - et est assignée à geneB. r4 couvre 350–380, entre les deux gènes, et n'est assignée à aucun gène. En mode non strandé, geneA reçoit trois lectures, geneB une, et r4 reste non assignée. En mode -s 1, une lecture + ne peut pas être assignée à geneB ; en mode -s 2, les lectures + ne peuvent pas être assignées à geneA.

Implémentation

Préparer un BAM synthétique

Le jeu datasets/ngs/demo.sam contient un fragment apparié pair1 et une lecture non appariée pair2, tous deux dans geneA. Pour illustrer le comptage de fragments, on conserve les paires correctement appariées avec le drapeau 0x2, puis on convertit en BAM. La commande suppose samtools installé.

mkdir -p results/rnaseq
samtools view -b -f 0x2 datasets/ngs/demo.sam > results/rnaseq/demo_paired.bam
samtools index results/rnaseq/demo_paired.bam

Le filtre -f 0x2 conserve les deux mates de pair1 et élimine pair2. Le BAM contient un seul fragment, avec deux lectures de 12 pb en 10–21 et 40–51, toutes deux dans le premier exon de geneA. Le fragment mesure 42 pb, sous le seuil minimal par défaut de 50 pb ; la commande abaisse donc ce seuil avec -d 20 pour ne pas l'écarter.

La commande exacte

La commande demande quatre threads, active le mode apparié et le comptage des paires, accepte les fragments dès 20 pb, ignore le brin, sélectionne les exons et regroupe par gene_id.

featureCounts \
  -T 4 \
  -p --countReadPairs \
  -d 20 \
  -s 0 \
  -t exon -g gene_id \
  -a datasets/rnaseq/annotation.gtf \
  -o results/rnaseq/counts_demo.txt \
  results/rnaseq/demo_paired.bam

Le fichier de sortie contient une ligne par gène et une colonne par BAM. La sortie illustrative suivante donne les colonnes essentielles ; Length est la longueur cumulée des exons, pas celle du transcrit.

Geneid  Chr   Start  End  Strand  Length  demo_paired.bam
geneA   chr1  1      300  +       251     1
geneB   chr1  400    700  -       301     0

Le fichier .summary ventile les lectures. Extrait illustratif :

Status                  demo_paired.bam
Assigned                1
Unassigned_NoFeatures   0
Unassigned_Ambiguity    0
Unassigned_MultiMapping 0

geneA vaut 1, et non 2, car les deux mates appartiennent au même fragment. Si l'on omet --countReadPairs avec Subread 2.0.6, chaque mate est compté et geneA reçoit 2.

Faire varier le brin

On relance la commande avec -s 1 puis -s 2, en changeant le fichier de sortie. pair1 est sur le brin + et geneA aussi. -s 1 compte le fragment ; -s 2 le rejette, car il faudrait un gène sur le brin opposé.

Ce résultat ne prouve pas que -s 1 convient à une vraie librairie : il dépend du protocole. Sur des données réelles, on estime l'orientation avec infer_experiment.py de RSeQC, ou en comparant les comptages -s 1 et -s 2 à des gènes de contrôle. Un taux d'assignation effondré dans les deux modes doit aussi faire suspecter l'annotation ou le génome.

Plusieurs échantillons

featureCounts accepte plusieurs BAM dans la même commande. La matrice contient alors une colonne par fichier, dans l'ordre fourni.

featureCounts \
  -T 8 \
  -p --countReadPairs \
  -s 1 \
  -t exon -g gene_id \
  -a annotation.gtf \
  -o counts_all.txt \
  sample1.bam sample2.bam sample3.bam

Les noms de colonnes sont les noms de fichiers, pas les noms biologiques. Il faut renommer les colonnes avec un tableau d'échantillons vérifié, jamais selon l'ordre alphabétique. Le fichier .summary doit être conservé : il indique la proportion de lectures assignées, ambiguës, multi-alignées, non alignées ou hors feature.

Vérifier le comptage en Python

Le code suivant reprend les décisions essentielles sur l'annotation synthétique : lecture des attributs, extraction des exons, union par gène et comptage d'un fragment à deux blocs.

import re
 
def parse_attributes(raw):
    attributes = {}
    for chunk in raw.split(";"):
        chunk = chunk.strip()
        if not chunk:
            continue
        match = re.match(r'(\S+)\s+"([^"]*)"', chunk)
        if match:
            attributes[match.group(1)] = match.group(2)
    return attributes
 
def parse_gtf(text):
    features = []
    for line in text.splitlines():
        if not line.strip() or line.startswith("#"):
            continue
        fields = line.split("\t")
        seqname, source, feature, start, end, score, strand, frame, raw = fields
        features.append({
            "seqname": seqname,
            "feature": feature,
            "start": int(start),
            "end": int(end),
            "strand": strand,
            "attributes": parse_attributes(raw),
        })
    return features
 
def assign_fragment(blocks, features, stranded=0):
    genes = set()
    for start, end in blocks:
        for exon in features:
            if exon["feature"] != "exon":
                continue
            if end < exon["start"] or start > exon["end"]:
                continue
            if stranded == 1 and exon["strand"] != "+":
                continue
            genes.add(exon["attributes"]["gene_id"])
    return genes
 
gtf_text = """chr1\tsynthetic\texon\t1\t150\t.\t+\t.\tgene_id "geneA";\nchr1\tsynthetic\texon\t200\t300\t.\t+\t.\tgene_id "geneA";"""
features = parse_gtf(gtf_text)
print(assign_fragment([(140, 210)], features))

Le résultat {'geneA'} confirme qu'un fragment couvrant deux exons du même gène n'est compté qu'une fois. Pour une matrice complète, il faut répéter l'opération par échantillon, puis sommer par gène et par colonne. Les exercices généralisent avec le brin, les fragments, les chevauchements et les lectures multi-alignées.

Étude de cas

Le jeu est synthétique et étiqueté comme tel. L'annotation contient deux gènes, trois exons et des lignes transcript et CDS qui ne doivent pas être comptées. Le BAM est dérivé de datasets/ngs/demo.sam en conservant les paires correctement appariées. La question : combien de fragments sont assignés à chaque gène, et que change la stratégie de brin ?

On inspecte l'annotation pour vérifier identifiants et coordonnées, on convertit le SAM en BAM filtré, puis on lance featureCounts avec -p --countReadPairs -d 20. On répète avec -s 1 et -s 2. Enfin, on lit la matrice et le .summary, et on compare au nombre de fragments du BAM.

Les résultats attendus sont ceux du tableau : geneA vaut 1 en non strandé et en forward, 0 en reverse ; geneB vaut 0 dans les trois cas. Le taux d'assignation est de 100 % pour le fragment conservé, ce qui est attendu sur un jeu construit pour tomber dans un exon. La longueur de geneA est 251 pb, soit 150 + 101, et non 300 pb : les exons sont séparés par un intron de 49 pb. Pour geneB, elle est de 301 pb.

L'interprétation reste prudente. Un comptage de 1 ne mesure pas l'expression de geneA ; il montre que la chaîne d'outils, l'annotation et le mode de brin sont cohérents. Sur un vrai jeu, une différence de quelques lectures entre -s 1 et -s 2 n'est pas nécessairement un signal : il faut examiner la distribution globale des taux d'assignation et l'orientation des lectures.

Le .summary est aussi instructif que la matrice. Si Unassigned_NoFeatures domine, l'annotation ne couvre pas les lectures ou le génome n'est pas le bon. Si Unassigned_Ambiguity domine, les gènes se chevauchent ou les isoformes sont concurrentes. Si Unassigned_MultiMapping domine, les régions répétées posent problème. Ces diagnostics orientent le choix des options ; ils ne se corrigent pas en augmentant arbitrairement -M ou -O.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-6/6.2-featurecounts -q

python -m pytest exercises/module-6/6.2-featurecounts -q

python -m pytest exercises/module-6/6.2-featurecounts/projet -q

Résumé

• featureCounts compte des lectures ou fragments sur des features d'un GTF ou d'un SAF ; un GFF3 doit être converti.
• -t exon -g gene_id regroupe les exons d'un gène en une méta-feature : une lecture à cheval sur deux exons est comptée une fois.
• -s 0, -s 1 et -s 2 correspondent aux modes non strandé, forward et reverse ; le choix se vérifie.
• Avec Subread 2.0.6, -p --countReadPairs compte un fragment apparié une seule fois ; -p seul compte chaque mate.
• -M autorise les lectures multi-alignées, -O les lectures chevauchant plusieurs gènes ; --fraction partage les comptes.
• Les comptages bruts ne sont ni normalisés par la longueur ni par la profondeur ; ils alimentent DESeq2, edgeR ou limma-voom.
• Le fichier .summary est indispensable pour diagnostiquer annotation, brin ou multi-mapping.

Références

• Liao, Y., Smyth, G. K. & Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7), 923–930.
• Liao, Y., Smyth, G. K. & Shi, W. (2019). The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Research, 47(8), e47.
• Subread Users Guide, version 2.0.6. https://subread.sourceforge.net/SubreadUsersGuide.pdf
• Documentation featureCounts. https://subread.sourceforge.net/featureCounts.html
• The Sequence Ontology. GFF3 specification. https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md
• Ensembl. GTF format documentation. https://www.ensembl.org/info/website/upload/gff.html
• Anders, S., Pyl, P. T. & Huber, W. (2015). HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics, 31(2), 166–169.