Pipeline RNA-Seq intégré

Prérequis

Cette section assemble des briques déjà vues : le comptage et la normalisation de 6.2 et 6.3, le modèle différentiel de 6.4 et l'enrichissement de 6.5. Elle suppose le cadrage du projet de la section 7.1 et le pipeline de variants de 3.5, présentés dans le module 7. Les identifiants de gènes, les formules de CPM et les sorties DESeq2 ne sont pas redéfinis : ils sont réutilisés, contrôlés et reliés.

Introduction

Une expérience RNA-Seq compare des populations de transcrits entre conditions. La question biologique est souvent simple — quels gènes changent, dans quel sens, et que disent ces changements du processus étudié ? — mais la réponse dépend d'une chaîne de traitements. Une librairie mal contrôlée, un comptage sur le mauvais brin, une normalisation inadaptée ou un contraste mal défini peuvent produire une liste de gènes convaincante et fausse. Le pipeline intégré consiste précisément à rendre chaque maillon visible, vérifiable et reproductible.

L'enjeu computationnel est double. D'abord, il faut faire dialoguer des formats et des outils : lectures brutes, alignements, annotation, matrice de comptages, métadonnées, résultats différentiels, termes fonctionnels. Ensuite, il faut distinguer ce qui relève de la description, de l'estimation et de l'interprétation. Un log2 fold change n'est pas une preuve mécanistique ; une p-value ajustée n'est pas une probabilité que le gène soit vrai ; un terme GO enrichi n'est pas une voie causale.

Cette section suit un enchaînement complet : contrôle qualité, quantification, normalisation, analyse DESeq2, puis enrichissement. Le fil conducteur est un jeu synthétique de quatre gènes et quatre échantillons, assez petit pour être vérifié à la main et trop petit pour soutenir une conclusion biologique forte. Nous verrons aussi comment ce pipeline s'articule avec l'appel de variants et comment un rapport reproductible conserve les versions, les paramètres et les limites.

Théorie

De la lecture au comptage

Le contrôle qualité intervient avant toute comparaison. FastQC résume la qualité par base, la composition, les adaptateurs, la duplication et la contamination ; fastp ou Trimmomatic retirent les adaptateurs et les bases de mauvaise qualité ; MultiQC rassemble les rapports. Trois signaux doivent alerter : une chute de qualité en fin de lecture, une proportion élevée de duplicats compatible avec une faible complexité, et une contamination par l'ARN ribosomal. Cette dernière est particulièrement trompeuse en RNA-Seq : elle consomme la profondeur sans informer sur les transcrits d'intérêt.

Vient ensuite le choix de la quantification. Un aligneur classique — STAR ou HISAT2 — place les lectures sur le génome, puis featureCounts compte les lectures qui tombent dans les exons d'un gène. Une méthode de quantification sélective — Salmon ou kallisto — travaille directement sur les transcrits et estime les abondances. Les deux familles peuvent alimenter une analyse différentielle, mais elles ne produisent pas les mêmes objets : featureCounts renvoie des comptages entiers par gène, Salmon renvoie des estimations par transcrit qu'il faut agréger.

Le comptage par gène exige une annotation cohérente avec la version du génome et une stratégie de brin correcte. Pour des lectures appariées, featureCounts se lance avec -p ; le paramètre -s vaut 0 pour une banque non strandée, 1 pour le brin direct et 2 pour le brin inverse. Se tromper de valeur ne produit pas toujours une erreur : le taux d'assignation chute, ou les comptages se déplacent vers des gènes chevauchants. Vérifiez le taux d'assignation et comparez les trois valeurs sur un sous-échantillon.

La matrice de comptages est un objet de contrôle autant qu'un résultat. On vérifie que les identifiants de gènes sont uniques, que les colonnes d'échantillons correspondent exactement aux métadonnées, que les valeurs sont des entiers positifs ou nuls, que les tailles de bibliothèque ne sont pas nulles et que les lignes vides ou dupliquées sont explicites. On documente aussi la version du GTF, la longueur des gènes si l'on calcule des TPM, et les éventuelles colonnes de batch. Une matrice propre mais mal annotée reste inutilisable.

Normaliser sans déformer

La normalisation la plus simple est le comptage par million (CPM) :

Elle corrige la profondeur de séquençage, pas la longueur du gène ni la composition. Les TPM et les FPKM corrigent en plus la longueur, ce qui les rend utiles pour comparer des gènes entre eux, mais pas pour alimenter un modèle de comptage. DESeq2 attend les comptages bruts et estime des facteurs de taille par la médiane des rapports entre chaque échantillon et une pseudo-référence. Cette méthode est plus robuste qu'une simple somme lorsqu'un petit nombre de gènes très exprimés domine la librairie.

Prenons geneA dans le jeu synthétique. Les tailles de bibliothèque valent 715, 733, 766 et 774. Les CPM de geneA valent environ 167 832 et 177 353 chez les témoins, puis 52 219 et 58 140 chez les traités. La moyenne des CPM passe d'environ 172 593 à 55 179, soit un rapport d'environ 0,32 et un log2 fold change descriptif de $-1{,}64$. Ce calcul est reproductible à la main ; il ne remplace pas le test différentiel, car il ignore la variance entre réplicats.

Le choix de la normalisation doit être cohérent avec la visualisation. Pour une heatmap de gènes, une transformation stabilisante de la variance (VST) ou un logarithme des CPM sont adaptés ; pour un test DESeq2, on garde les comptages bruts.

Modèle, design et contrastes

DESeq2 modélise le comptage du gène $i$ dans l'échantillon $j$ par une binomiale négative :

Le paramètre $s_j$ est le facteur de taille, $\alpha_i$ la dispersion du gène et $x_j$ la variable de condition. La dispersion capture la surdispersion : la variance vaut $\mu_{ij} + \alpha_i \mu_{ij}^2$, supérieure à celle d'une loi de Poisson. DESeq2 estime d'abord une dispersion par gène, puis la rapproche d'une tendance moyenne lorsque le nombre de réplicats est faible. C'est ce compromis qui rend l'analyse possible avec deux ou trois réplicats, au prix d'une incertitude qu'il faut reconnaître.

Le design est la partie la plus facile à mal écrire. Une formule ~ condition compare les conditions en tenant compte d'un niveau de référence ; ~ batch + condition ajuste un effet de lot ; ~ condition + condition:batch ajoute une interaction. Le niveau de référence détermine le signe des contrastes : avec control comme référence, un coefficient positif signifie une hausse chez les traités. Un design confondu — par exemple tous les témoins dans un lot et tous les traités dans un autre — ne se rattrape pas par le modèle : l'effet biologique et l'effet technique sont indiscernables.

Un contraste se définit explicitement, par exemple contrast = c("condition", "treated", "control"). Il ne faut pas confondre le coefficient affiché par resultsNames() et l'interprétation biologique : un coefficient d'interaction ne se lit pas comme un effet principal.

Tester, corriger, mesurer l'effet

Le test de Wald compare l'estimation du coefficient à son erreur standard : $z = \hat{\beta} / \mathrm{SE}(\hat{\beta})$. La p-value associée teste l'hypothèse nulle $\beta = 0$ pour un gène donné. Tester des milliers de gènes produit mécaniquement des faux positifs : au seuil de 0,05, environ 5 % des gènes réellement nuls seraient déclarés significatifs. La procédure de Benjamini-Hochberg contrôle le taux de fausses découvertes :

où $V$ est le nombre de faux positifs et $R$ le nombre de rejets. La colonne padj de DESeq2 est une p-value ajustée par cette procédure, après filtrage indépendant. Elle ne signifie pas « probabilité que le gène soit faux » : elle quantifie une erreur attendue dans l'ensemble des gènes déclarés significatifs.

La taille d'effet est aussi importante que la significativité. Un gène très significatif avec un log2 fold change de 0,2 peut être biologiquement négligeable ; un gène à fort effet mais à variance élevée peut manquer de puissance. DESeq2 propose lfcShrink pour stabiliser les estimations, en particulier aux faibles comptages. Le paramètre lfcThreshold permet de tester un effet minimal plutôt que zéro, ce qui change l'hypothèse nulle et doit être annoncé. Enfin, un gène marqueur demande un effet reproductible, une expression suffisante et, si possible, une validation indépendante.

Enrichir sans surestimer

L'enrichissement fonctionnel demande si une liste de gènes contient plus de gènes associés à un terme GO ou à une voie KEGG que ce que l'univers de référence produirait au hasard. Le test hypergéométrique compare le chevauchement observé à l'attendu ; clusterProfiler fournit enrichGO, enrichKEGG et enricher pour un catalogue personnalisé. Trois choix déterminent le résultat : la liste d'intérêt, l'univers (tous les gènes testés ou tous les gènes annotés ?), et la source d'annotation.

Un enrichissement n'est pas une expérience. Il ne démontre ni causalité ni spécificité : un terme très général, comme « processus métabolique », peut être enrichi parce que la liste est petite ou parce que l'annotation est biaisée. Les termes redondants doivent être regroupés ou lus ensemble. Avec quatre gènes différentiels, l'univers est minuscule, les termes sont peu nombreux et le FDR est instable : mieux vaut présenter le résultat comme une piste que comme une conclusion.

Implémentation

Le script R suivant enchaîne lecture, contrôle, filtrage, DESeq2, diagnostics et export. Il est écrit pour le jeu synthétique datasets/rnaseq/counts.tsv, avec deux témoins et deux traités. Les versions cibles sont R 4.4+, Bioconductor 3.19+, DESeq2 1.38+ et featureCounts 2.0.6. Le bloc d'enrichissement est volontairement conditionnel : les identifiants synthétiques ne sont pas des identifiants Entrez, et un échec de mapping ne doit pas interrompre le reste du rapport.

# pipeline_rnaseq.R — analyse différentielle du jeu synthétique counts.tsv
suppressPackageStartupMessages({
  library(DESeq2)
})
 
counts_path <- file.path("datasets", "rnaseq", "counts.tsv")
out_dir <- file.path("results", "rnaseq")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
 
raw <- read.delim(counts_path, check.names = FALSE, stringsAsFactors = FALSE)
metadata <- data.frame(
  sample = c("S1", "S2", "S3", "S4"),
  condition = factor(
    c("control", "control", "treated", "treated"),
    levels = c("control", "treated")
  ),
  stringsAsFactors = FALSE
)
rownames(metadata) <- metadata$sample
 
stopifnot(all(metadata$sample %in% colnames(raw)))
count_mat <- as.matrix(raw[, metadata$sample, drop = FALSE])
rownames(count_mat) <- raw$gene_id
storage.mode(count_mat) <- "integer"
stopifnot(!anyNA(count_mat), all(count_mat >= 0))
stopifnot(!anyDuplicated(rownames(count_mat)))
stopifnot(identical(colnames(count_mat), rownames(metadata)))
 
dds <- DESeqDataSetFromMatrix(
  countData = count_mat,
  colData = metadata,
  design = ~ condition
)
dds <- dds[rowSums(counts(dds) >= 10) >= 2, ]
dds$condition <- relevel(dds$condition, ref = "control")
dds <- DESeq(dds)
 
res <- results(
  dds,
  contrast = c("condition", "treated", "control"),
  alpha = 0.05
)
coef_name <- grep("^condition_.*_vs_control$", resultsNames(dds), value = TRUE)
stopifnot(length(coef_name) == 1L)
res_shrunk <- lfcShrink(dds, coef = coef_name, type = "normal")
 
summary(res)
sum(res$padj < 0.05, na.rm = TRUE)
 
write.csv(
  as.data.frame(res_shrunk),
  file = file.path(out_dir, "deseq2_treated_vs_control.csv"),
  row.names = TRUE
)
 
png(file.path(out_dir, "ma_plot.png"), width = 900, height = 700)
plotMA(res, ylim = c(-3, 3), main = "MA plot — treated vs control")
dev.off()
 
# vst peut échouer sur un jeu minuscule ; normTransform sert alors de repli.
vsd <- tryCatch(
  vst(dds, blind = FALSE),
  error = function(e) normTransform(dds)
)
png(file.path(out_dir, "pca.png"), width = 900, height = 700)
print(plotPCA(vsd, intgroup = "condition"))
dev.off()
 
if (requireNamespace("pheatmap", quietly = TRUE)) {
  top <- head(order(res$padj, na.last = TRUE), 20)
  top <- top[!is.na(res$padj[top])]
  if (length(top) > 1) {
    png(file.path(out_dir, "heatmap_top.png"), width = 1000, height = 800)
    pheatmap::pheatmap(
      assay(vsd)[top, , drop = FALSE],
      scale = "row",
      annotation_col = as.data.frame(
        colData(dds)[, "condition", drop = FALSE]
      ),
      main = "Top gènes différentiels"
    )
    dev.off()
  }
}
 
# Enrichissement : à exécuter avec de vrais identifiants mappés.
if (requireNamespace("clusterProfiler", quietly = TRUE) &&
    requireNamespace("org.EcK12.eg.db", quietly = TRUE)) {
  suppressPackageStartupMessages({
    library(clusterProfiler)
    library(org.EcK12.eg.db)
  })
  sig <- rownames(res)[
    !is.na(res$padj) & res$padj < 0.05 &
      abs(res$log2FoldChange) >= 1
  ]
  entrez <- tryCatch(
    suppressWarnings(
      mapIds(
        org.EcK12.eg.db,
        keys = sig,
        keytype = "SYMBOL",
        column = "ENTREZID"
      )
    ),
    error = function(e) setNames(rep(NA_character_, length(sig)), sig)
  )
  entrez <- unname(entrez[!is.na(entrez)])
  if (length(entrez) > 0) {
    ego <- enrichGO(
      gene = entrez,
      OrgDb = org.EcK12.eg.db,
      keyType = "ENTREZID",
      ont = "BP",
      pAdjustMethod = "BH",
      pvalueCutoff = 0.05,
      qvalueCutoff = 0.2,
      readable = TRUE
    )
    write.csv(
      as.data.frame(ego),
      file = file.path(out_dir, "enrichGO_BP.csv"),
      row.names = FALSE
    )
    png(file.path(out_dir, "enrichGO_dotplot.png"), width = 1000, height = 700)
    print(dotplot(ego, showCategory = 15))
    dev.off()
  } else {
    message("Aucun identifiant mappé : enrichissement ignoré.")
  }
} else {
  message("clusterProfiler ou org.EcK12.eg.db absent : enrichissement ignoré.")
}
 
sessionInfo()

Les sorties R ne sont pas recopiées : avec quatre gènes, une dispersion stable n'a pas de sens biologique. Le format attendu de summary(res) est le suivant ; les valeurs sont omises.

# Sortie illustrative — non exécutée
# out of ... with nonzero total read count
# adjusted p-value < 0.05
# LFC > 0 (up)     : ...
# LFC < 0 (down)   : ...
# outliers [1]     : ...
# low counts [2]   : ...

Le point critique est le filtrage : rowSums(counts(dds) >= 10) >= 2 retire les gènes qui n'ont pas au moins deux comptages suffisants. Ce filtre précède le modèle et ne regarde pas les p-values ; il diffère du filtrage indépendant que DESeq2 applique ensuite. Un autre point est le contraste : coef_name est retrouvé par expression régulière plutôt que codé en dur, car le nom du coefficient dépend des niveaux du facteur. Enfin, lfcShrink avec type = "normal" ne demande pas de package supplémentaire ; apeglm ou ashr sont des alternatives qui modifient la forme de la shrinkage et doivent être documentées.

L'enchaînement peut être décrit dans un Snakefile. L'extrait suivant suppose des BAM triés dans results/align/ et une annotation GTF ; il montre les règles, pas la totalité du projet.

SAMPLES = ["S1", "S2", "S3", "S4"]
 
rule all:
    input:
        "results/rnaseq/deseq2_treated_vs_control.csv",
        "reports/rnaseq.html"
 
rule featurecounts:
    input:
        bams = expand("results/align/{sample}.bam", sample=SAMPLES),
        gtf = "datasets/rnaseq/annotation.gtf"
    output:
        counts = "results/counts/featurecounts.tsv",
        summary = "results/counts/featurecounts.tsv.summary"
    params:
        strandedness = 2
    threads: 4
    shell:
        "featureCounts -T {threads} -p -s {params.strandedness} "
        "-a {input.gtf} -o {output.counts} {input.bams}"
 
rule counts_matrix:
    input:
        "results/counts/featurecounts.tsv"
    output:
        "results/counts/matrix.tsv"
    shell:
        "cut -f1,7- {input} > {output}"
 
rule deseq2:
    input:
        counts = "results/counts/matrix.tsv",
        metadata = "datasets/rnaseq/metadata.tsv"
    output:
        "results/rnaseq/deseq2_treated_vs_control.csv"
    script:
        "scripts/deseq2.R"
 
rule report:
    input:
        "results/rnaseq/deseq2_treated_vs_control.csv"
    output:
        "reports/rnaseq.html"
    shell:
        "quarto render reports/rnaseq.qmd --output rnaseq.html"

La règle counts_matrix mérite un commentaire : la sortie de featureCounts contient cinq colonnes d'annotation avant les colonnes d'échantillons, et son en-tête reprend les chemins des BAM. Le cut garde Geneid et les comptages, mais il faut renommer les colonnes dans R avec les vrais noms d'échantillons. La valeur strandedness = 2 est un exemple : elle doit être vérifiée sur la banque, pas recopiée.

Le rapport Quarto peut recevoir le contraste en paramètre, ce qui évite de dupliquer le document pour chaque comparaison.

---
title: "Rapport RNA-Seq"
params:
  contrast: "treated_vs_control"
format:
  html:
    code-fold: true
    toc: true
execute:
  freeze: auto
---

Le rendu s'effectue avec quarto render reports/rnaseq.qmd --execute-params params.yml. Le paramètre doit être lu dans le code R et affiché dans le rapport, pas seulement dans le titre. Avec freeze: auto, Quarto réutilise les résultats ; cela ne dispense pas de conserver versions et données.

Étude de cas

Le jeu contient quatre gènes et quatre échantillons. Les échantillons S1 et S2 sont des témoins ; S3 et S4 sont traités. Les comptages sont les suivants.

La question est double : le pipeline retrouve-t-il les directions attendues, et que peut-on légitimement conclure ? Les tailles de bibliothèque valent 715, 733, 766 et 774 ; aucune n'est vide, mais geneD domine avec environ 500 comptages, soit près des deux tiers de chaque librairie. Cette composition justifie de regarder les CPM et de ne pas comparer les comptages bruts entre échantillons.

Le filtrage min_count = 10, min_samples = 2 conserve les quatre gènes : geneC atteint 14 à 18 comptages dans les quatre échantillons. Un seuil plus strict, min_count = 100, retirerait geneC. Ce choix doit être annoncé : il ne signifie pas que geneC est biologiquement nul, seulement qu'il est peu informatif pour ce modèle. Les moyennes de CPM donnent geneA en baisse d'un facteur d'environ 3 et geneB en hausse d'un facteur d'environ 2,6. geneC et geneD restent proches de leur niveau témoin. Ces calculs sont descriptifs et reproductibles ; ils ne sont pas des p-values.

Avec deux réplicats par condition, DESeq2 peut estimer une dispersion, mais l'incertitude est grande. Le script fournit les figures et le tableau attendus ; les valeurs de padj ne doivent pas être présentées comme une preuve. Le contraste treated contre control devrait retrouver geneA négatif et geneB positif, avec des intervalles larges. geneC, très peu compté, peut être filtré indépendamment et recevoir un padj manquant. geneD, très exprimé, peut être significatif pour un effet faible : c'est le cas typique où la taille d'effet doit être lue avant la p-value.

L'interprétation reste prudente. geneA et geneB sont des candidats marqueurs dans ce jeu synthétique ; ils ne prouvent rien sur un mécanisme. Aucun enrichissement GO n'est interprétable avec quatre gènes : l'univers est trop petit, les termes seraient choisis après observation et le FDR serait instable. Pour une vraie expérience, on vérifierait la qualité des librairies, la cohérence des réplicats par ACP, l'effet de lot, puis on confirmerait les candidats par qPCR ou par un jeu indépendant.

Le croisement avec les variants exige les mêmes identifiants et la même annotation. On intersecte les gènes différentiels avec les variants annotés de 3.5. Un variant dans le promoteur de geneA ou dans la région codante de geneB est une hypothèse, pas une explication : il faut vérifier la profondeur, la qualité de l'appel, le type de variant et son effet attendu. L'expression allélique spécifique peut renforcer l'hypothèse, mais l'appel de variants sur RNA-Seq reste sensible à l'épissage, aux duplications et à l'expression différentielle des allèles.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-7/7.3-pipeline-rnaseq -q

python -m pytest exercises/module-7/7.3-pipeline-rnaseq -q

python -m pytest exercises/module-7/7.3-pipeline-rnaseq/projet -q

Résumé

• Le pipeline RNA-Seq enchaîne contrôle qualité, quantification, normalisation, modèle différentiel et enrichissement ; chaque étape possède ses contrôles et ses formats.
• Une matrice de comptages se vérifie avant l'analyse : identifiants uniques, colonnes cohérentes, entiers positifs, tailles de bibliothèque non nulles, annotation versionnée.
• DESeq2 attend des comptages bruts et estime des facteurs de taille ; les CPM, TPM et transformations logarithmiques servent aux descriptions et aux figures.
• Le design et le contraste doivent être écrits et vérifiés avant le modèle ; un facteur confondu avec le lot ne se corrige pas après coup.
• Le FDR, la taille d'effet et les comptages normalisés se lisent ensemble ; padj n'est pas une probabilité que le gène soit faux.
• Un gène marqueur est un candidat robuste, pas le gène le plus significatif d'une liste ; il demande un effet reproductible et, si possible, une validation indépendante.
• L'enrichissement sur une liste courte reste descriptif ; l'univers, l'annotation et la redondance des termes limitent l'interprétation.
• Le croisement avec les variants exige les mêmes identifiants et la même annotation ; il produit des hypothèses, pas une causalité.

Références

• Love, M. I., Huber, W. et Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550.
• Liao, Y., Smyth, G. K. et Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7), 923–930.
• Yu, G., Wang, L.-G., Han, Y. et He, Q.-Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5), 284–287.
• Robinson, M. D. et Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology, 11(3), R25.
• Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L. et Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods, 5(7), 621–628.
• Mölder, F. et al. (2021). Sustainable data analysis with Snakemake. F1000Research, 10, 33.
• Documentation Bioconductor. https://bioconductor.org/packages/DESeq2/