Préparation de librairie RNA-Seq et quantification

Prérequis

Cette section suppose connus les fichiers FASTQ et les scores Phred présentés en 1.2, ainsi que la démarche de contrôle qualité NGS de 3.1. Les drapeaux SAM et les fichiers BAM sont introduits en 3.3 ; ils serviront ici à vérifier l'orientation des lectures. Le transcriptome de référence et les comptages produits dans cette section seront réutilisés dans la suite du module, notamment pour la normalisation et l'analyse différentielle.

Introduction

Le transcriptome est l'ensemble des ARN présents dans une cellule à un instant donné. Il change avec le type cellulaire, l'environnement, le cycle cellulaire et le traitement appliqué. Contrairement à l'ADN génomique, dont la copie est globalement fixe, l'ARN messager est produit, épissé, transporté et dégradé en permanence. Son abondance peut varier de plusieurs ordres de grandeur entre deux gènes, et un même gène peut produire plusieurs isoformes. Le séquençage de l'ARN (RNA-Seq) transforme cette population en millions de courtes lectures, puis en comptages par gène ou par transcrit. La qualité de cette transformation dépend d'abord de la librairie : une molécule absente, dégradée ou mal orientée avant séquençage ne sera pas reconstituée par l'analyse informatique.

La préparation de librairie enchaîne des choix biologiques et enzymatiques. Il faut extraire un ARN intact, décider si l'on enrichit les ARN messagers polyadénylés ou si l'on retire l'ARN ribosomique (ARNr), fragmenter, copier en ADN complémentaire, marquer le brin, ligaturer des adaptateurs et amplifier. Chaque étape introduit un biais : sélection de taille, conversion incomplète, biais de composition GC, duplicates de PCR, ARNr résiduel. Ces biais ne sont pas de simples détails techniques ; ils déterminent quelles régions du transcriptome sont observables et avec quelle précision.

Sur le plan computationnel, le défi est de passer de lectures à des quantités comparables. Il faut connaître l'orientation de la librairie, distinguer les fragments des lectures, attribuer les lectures ambiguës, mesurer la saturation et la redondance, puis choisir entre comptage par gène et quantification de transcrits. Nous suivons donc la chaîne complète, de l'extraction au comptage, avant d'aborder la normalisation et les tests différentiels.

Théorie

De l'ARN total à l'ADN complémentaire

L'extraction doit inactiver les RNases, enzymes stables qui dégradent l'ARN. Les méthodes au phénol-chloroforme et les colonnes de silice ne conservent pas les mêmes populations : les premières favorisent les longs transcrits, les secondes peuvent retenir les petits ARN. Un traitement à la DNase élimine l'ADN génomique, dont la contamination produirait des lectures introniques et des comptages faussement élevés. Le RIN (RNA Integrity Number) résume l'intégrité électrophorétique ; un RIN élevé ne garantit pas chaque transcrit.

La transcription inverse convertit l'ARN en ADN complémentaire (ADNc). Les amorces oligo-dT s'hybrident sur la queue poly-A ; les amorces aléatoires démarrent partout et réduisent le biais vers l'extrémité 3', au prix de fragments plus courts et d'une composition dépendante du hasard. La transcriptase inverse peut ajouter quelques nucléotides non appariés, comportement exploité par certains protocoles pour conserver le brin. La deuxième chaîne est synthétisée par une ADN polymérase, puis les adaptateurs sont ligaturés.

La fragmentation, enzymatique ou chimique, précède ou suit la transcription inverse. Une sélection de taille retire les inserts trop courts ou trop longs. En séquençage apparié, chaque fragment est lu aux deux extrémités : on distingue les fragments des lectures et on améliore l'alignement sur les jonctions d'épissage. Un insert plus court qu'une lecture produit une traversée d'adaptateur, que fastp détecte et retire.

Poly-A ou déplétion de l'ARNr

La sélection poly-A capture les ARN messagers matures avec des billes oligo-dT. Elle est efficace et peu coûteuse chez les eucaryotes, mais exclut les ARN non polyadénylés, les ARN histone et une partie des ARN longs non codants. Sur un ARN dégradé, elle accentue le biais 3' : les fragments sans queue poly-A sont perdus.

La déplétion retire les ARNr par hybridation de sondes puis digestion ou capture. Elle conserve les ARN non polyadénylés, les pré-ARN messagers et les ARN bactériens ; elle est donc adaptée à Escherichia coli, aux ARN non codants et aux échantillons dégradés. Son revers est un ARNr résiduel variable et un coût plus élevé. La fraction d'ARNr se mesure après séquençage : $f_{\mathrm{ARNr}} = n_{\mathrm{ARNr}} / n_{\mathrm{lectures}}$. Une valeur élevée signale un enrichissement insuffisant ou un ARN de départ dégradé.

Le choix dépend de l'organisme et de la question. Pour des transcrits matures humains, la sélection poly-A est standard ; pour les bactéries, les ARN non codants ou les échantillons fixés, la déplétion est préférable. Aucune stratégie n'est universelle : archivez kit, entrée et qualité de l'ARN.

Brin, taille des fragments et duplicates

Une librairie non orientée ne distingue pas le brin sens de l'antisens. Les protocoles orientés marquent un brin, souvent par incorporation de dUTP lors de la deuxième chaîne puis dégradation sélective. featureCounts règle l'orientation avec -s : 0 non orienté, 1 direct, 2 inverse. Un mauvais réglage divise parfois les comptages par deux et crée des différences entre gènes chevauchants.

Les drapeaux SAM aident à contrôler l'orientation. Le bit 0x10 indique le brin inverse ; 0x40 et 0x80 identifient la première et la seconde lecture d'une paire. Un drapeau 99 (0x63) décrit une première lecture forward, un drapeau 147 (0x93) une seconde reverse. Ces bits ne déduisent pas la stratégie de la librairie : il faut les comparer à un gène de référence.

La taille des fragments influence la puissance et les biais. Des inserts courts favorisent les lectures traversantes et l'adaptateur ; des inserts longs peuvent être sous-représentés. Le nombre de cycles PCR doit rester minimal : chaque cycle ajoute des duplicates, issus d'une même molécule amplifiée. Les UMI (unique molecular identifiers), codes aléatoires ajoutés avant amplification, permettent de les reconnaître, mais tous les protocoles n'en ont pas.

Réplicats, profondeur et complexité

Un réplicat technique répète la même librairie ou le même séquençage ; il mesure la variabilité de la machine. Un réplicat biologique répète l'expérience sur un autre individu ou une autre culture ; il mesure la variabilité du système. Les réplicats biologiques sont indispensables pour estimer la dispersion et tester une différence. Trois par condition sont un point de départ, pas une garantie.

La profondeur est le nombre de lectures ou de fragments produits. Elle détermine la détection des transcrits rares et la précision des comptages, mais ne se confond pas avec la complexité de la librairie, le nombre de molécules uniques avant amplification. Une librairie pauvre en molécules uniques peut être séquencée très profondément : les lectures supplémentaires redisent la même information et augmentent les duplicates. La courbe de saturation sépare les deux situations.

On résume la distribution des comptages par l'entropie de Shannon. Pour $G$ gènes de comptages $c_i$ et de probabilités $p_i = c_i / \sum_k c_k$, on pose $H = -\sum_i p_i \log_2 p_i$. L'entropie normalisée $H / \log_2 G$ vaut 1 lorsque les lectures sont également réparties et tend vers 0 lorsqu'un gène domine. C'est un indicateur de diversité, pas une estimation des molécules uniques : il ne distingue pas les duplicates PCR de l'abondance biologique. Il complète le taux de duplication, les UMI et la saturation.

Biais GC, ARNr résiduel et contrôle qualité

La composition GC influence la fragmentation, l'hybridation et l'amplification ; les régions très riches ou très pauvres peuvent être sous-représentées. FastQC compare la distribution GC des lectures à une distribution théorique. Un décalage ou un pic étroit s'interprète avec la composition du transcriptome. fastp signale les séquences surreprésentées, mais aucune correction ne remplace une librairie équilibrée.

L'ARNr résiduel reste un indicateur direct de la préparation : on estime sa fraction après alignement ou on repère des séquences ribosomiques surreprésentées dans FastQC. Une fraction faible ne garantit pas l'absence d'ADN génomique : les lectures introniques et les séquences répétées doivent être examinées. Les adaptateurs non retirés faussent l'alignement et le comptage ; fastp les détecte, mais les paramètres doivent être vérifiés. Un statut FastQC vert peut cacher une faible complexité, une contamination ou un biais de brin : lisez les modules et comparez les échantillons.

Comptage par gène ou quantification de transcrits

featureCounts attribue les lectures aux exons d'une annotation puis les agrège par gène avec -g gene_id. C'est un comptage de lectures par gène, rapide et interprétable, adapté à une analyse différentielle au niveau du gène. Les lectures ambiguës, chevauchant plusieurs gènes, sont traitées selon les options ; le fichier .summary indique les lectures assignées, ambiguës et non comptées. En mode apparié, -p active les paires et --countReadPairs compte les fragments.

La quantification de transcrits poursuit un autre objectif. Salmon ou kallisto estiment l'abondance des isoformes à partir des lectures et d'un modèle d'inférence ; ils produisent des comptages estimés ou des TPM. L'approche est plus fine mais plus fragile : les isoformes partagent des exons et l'inférence dépend du modèle. Un comptage par gène ne permet pas de conclure sur l'usage des exons ni sur une isoforme minoritaire.

Le comptage brut n'est pas une normalisation : il dépend de la profondeur, de la longueur des gènes et de la composition. DESeq2 attend des comptages bruts et modélise ces effets ; lui donner des TPM ou des FPKM est une erreur.

Ce qui change par rapport au séquençage génomique

Le séquençage génomique mesure un génome de copie globalement constante. Le RNA-Seq mesure une population dynamique, avec des abondances variables, des transcrits épissés et des isoformes. Les lectures chevauchent des jonctions d'épissage, ce qui exige STAR, HISAT2 ou un pseudo-alignement. L'annotation devient centrale : un gène absent n'est pas compté.

Le RNA-Seq est aussi plus sensible aux contaminations (ARNr, ADN, adaptateurs, duplicates, GC), la couverture n'est pas uniforme, et l'orientation comme la taille des fragments sont des paramètres importants.

Exemple travaillé à la main

Le fichier datasets/rnaseq/counts.tsv contient quatre gènes et quatre échantillons synthétiques. Pour S1, les comptages 120, 80, 15 et 500 totalisent $N_1 = 715$. Les fractions valent $120/715 \approx 0.168$, $80/715 \approx 0.112$, $15/715 \approx 0.021$ et $500/715 \approx 0.699$. Pour S3, le total est 766 et les fractions sont 0.052, 0.287, 0.021 et 0.640. Alpha paraît divisé par environ 3,2, beta multiplié par environ 2,6, et delta reste majoritaire avec une part qui baisse. Les comptages bruts de delta (500 contre 490) semblent stables, alors que sa part relative change : la composition compte autant que le total.

Pour l'entropie, prenons deux gènes de comptages 75 et 25. Les probabilités sont 0,75 et 0,25 ; $H = -0.75 \log_2 0.75 - 0.25 \log_2 0.25 \approx 0.811$ bit. Comme $\log_2 2 = 1$, la complexité normalisée vaut environ 0,811. Avec quatre gènes équilibrés, $H = 2$ bits et la complexité vaut 1. Avec les comptages agrégés (335, 585, 63, 2005), elle tombe à 0,6591 : la distribution est concentrée.

Implémentation

Le pipeline minimal comporte quatre étapes : contrôle qualité brut avec FastQC, nettoyage avec fastp, alignement avec un aligneur épissé, puis comptage avec featureCounts. Les commandes suivantes sont illustratives et supposent les outils installés ; options conformes au frontmatter.

fastqc --threads 4 --outdir qc/ raw/S1_R1.fastq.gz raw/S1_R2.fastq.gz
 
fastp --in1 raw/S1_R1.fastq.gz --in2 raw/S1_R2.fastq.gz \
  --out1 trim/S1_R1.fastq.gz --out2 trim/S1_R2.fastq.gz \
  --detect_adapter_for_pe --qualified_quality_phred 20 \
  --length_required 36 --json qc/S1.fastp.json --html qc/S1.fastp.html
 
featureCounts -T 4 -p --countReadPairs -s 2 -t exon -g gene_id \
  -a annotation/gencode.gtf -o counts/all_samples.txt aligned/*.bam

FastQC produit un rapport HTML. fastp retire les adaptateurs, filtre les bases de faible qualité et écrit un rapport JSON lisible par machine. featureCounts exige une annotation au format GTF ou GFF et un BAM trié par coordonnées. L'option -s 2 doit être adaptée à la librairie ; -p --countReadPairs compte les fragments. Le fichier .summary donne les lectures assignées, ambiguës et non comptées.

Le code Python suivant reproduit les métriques essentielles sur de petites données. Il est autonome et ne lit aucun fichier : les comptages sont ceux du fichier counts.tsv, agrégés par gène. La sortie a été vérifiée le 2026-11-01 avec Python 3.13.

Aucune sortie. Exécutez le code pour voir le résultat.

La sortie vérifiée est :

0.75
0.6591
('first', '+') ('second', '-')
40.0

Trois points méritent attention. D'abord, rrna_fraction compte une lecture si elle est une sous-chaîne exacte d'un ARNr : c'est un test de dépistage, pas un alignement. Ensuite, la complexité normalisée ne remplace pas une estimation de duplicates : elle décrit la forme de la distribution des comptages. Enfin, strand ne devine pas la stratégie de la librairie ; il lit les bits du drapeau SAM. Pour régler featureCounts, il faut confronter cette information au kit et à un contrôle sur des gènes connus.

Étude de cas

Jeu public. Le paquet Bioconductor airway contient un jeu RNA-Seq de cellules musculaires lisses des voies aériennes humaines. Il comprend huit bibliothèques appariées : quatre lignées cellulaires, chacune observée dans une condition traitée par la dexaméthasone et dans une condition non traitée. Les métadonnées cell et dex décrivent l'appariement et le traitement ; la matrice counts(airway) fournit des comptages par gène. Le paquet ne contient pas les FASTQ bruts, mais il sert de référence réelle pour la normalisation et l'analyse différentielle des sections suivantes. Nous n'inventons ici aucune statistique : les chiffres cités pour l'étude de cas proviennent du fichier synthétique local.

Jeu synthétique local. datasets/rnaseq/counts.tsv contient quatre gènes et quatre échantillons. Les totaux sont 715, 733, 766 et 774 lectures. Le gène delta domine avec des fractions de 0,699, 0,696, 0,640 et 0,652 ; alpha passe de 0,168 à 0,052 entre S1 et S3, tandis que beta passe de 0,112 à 0,287. gamma reste autour de 0,02. La complexité normalisée par échantillon varie de 0,628 à 0,637, et celle des comptages agrégés vaut 0,6591. Ces valeurs sont des sorties vérifiées du code ci-dessus et du fichier local.

La question est de savoir si les comptages bruts sont directement comparables. Non : les profondeurs diffèrent (715 à 774), la composition change (alpha et beta s'échangent presque leur part relative, delta reste majoritaire), et le plan n'est pas un simple ensemble de réplicats. La démarche vérifie les métadonnées et l'unité de comptage, inspecte totaux et fractions, mesure la complexité, contrôle la qualité technique si les FASTQ sont disponibles, puis choisit une normalisation et un modèle adaptés.

Pour airway, les FASTQ ne sont pas fournis : on travaille sur les comptages et les métadonnées. L'appariement impose un modèle avec effet de lignée, et le faible nombre de réplicats limite la puissance. La normalisation doit tenir compte de la composition. Le comptage par gène convient à une question au niveau du gène ; une question sur les isoformes demanderait une quantification de transcrits et une validation indépendante. Quatre gènes et quatre échantillons ne suffisent pas pour conclure à une différence biologique : le jeu synthétique ne possède pas de réplicats explicites. Il montre pourquoi les comptages bruts ne sont pas des proportions et pourquoi une matrice doit être accompagnée de métadonnées, d'une annotation et d'un rapport de qualité.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-6/6.1-preparation-librairie-rnaseq -q

python -m pytest exercises/module-6/6.1-preparation-librairie-rnaseq -q

python -m pytest exercises/module-6/6.1-preparation-librairie-rnaseq/projet -q

Résumé

• La préparation conditionne l'observabilité : extraction, enrichissement, fragmentation, ADNc, brin et amplification introduisent des biais durables.
• La sélection poly-A convient aux ARN matures eucaryotes ; la déplétion de l'ARNr est nécessaire pour les bactéries et les ARN non codants.
• Le brin se vérifie sur le kit et les drapeaux SAM ; featureCounts attend -s 0, -s 1 ou -s 2.
• Les réplicats biologiques mesurent la variabilité du système ; les techniques ne la remplacent pas.
• La profondeur ne corrige pas une librairie pauvre en molécules uniques ; duplicates, saturation et complexité sont distincts.
• ARNr résiduel, biais GC et adaptateurs se contrôlent avec FastQC et fastp, mais aucun statut global ne remplace la lecture des modules.
• featureCounts produit des comptages par gène ; la quantification de transcrits est une inférence différente, adaptée aux isoformes.
• Les comptages bruts ne sont pas comparables : la normalisation et le modèle différentiel arrivent ensuite.

Références

• 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.
• Robinson, M. D. et Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology, 11(3), R25.
• 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.
• Chen, S., Zhou, Y., Chen, Y. et Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890.
• Documentation Bioconductor du paquet airway. https://bioconductor.org/packages/airway/
• Documentation FastQC. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• Documentation Subread et featureCounts. https://subread.sourceforge.net/