Impact des algorithmes d’analyse de données ARN-seq sur l’estimation de l’expression génique et la prédiction en aval

Ensembles de données de référence FDA SEQC

L’ensemble de données FDA SEQC-benchmark (numéro d’accession Omnibus d’expression génique GSE47792) comprend des données ARN-seq à extrémité appariée générées à l’aide de la plate-forme Illumina HiSeq 2000 avec une longueur de lecture de 100 nucléotides7. Nous avons utilisé un sous-ensemble de données SEQC -benchmark séquencées sur deux sites – Beijing Genomics Institute (BGI) et Mayo Clinic (MAY). Nous avons utilisé quatre échantillons (c.-à-D., A, B, C et D), chacun avec quatre bibliothèques répliquées préparées sur les sites de séquençage. L’échantillon A contient l’ARN de Référence Humain Universel (UHRR), l’échantillon B contient l’ARN de Référence du Cerveau Humain (HBRR), l’échantillon C contient un mélange de A et B (75% A et 25% B) et l’échantillon D contient un mélange de A et B (25% A et 75% B). Nous avons utilisé les données de deux voies d’une seule cellule d’écoulement pour chaque réplique d’échantillon. Le SEQC a également fourni un ensemble de données de référence de PCR quantitative (qPCR) qui comprend 20 801 gènes testés avec PrimePCR (Bio-Rad, Hercules, Californie). Chaque gène PrimePCR a été dosé une fois pour chacun des quatre échantillons (c.-à-D., A, B, C et D). Les ensembles de données et les échantillons de référence de la FDA SEQC sont résumés dans les tableaux supplémentaires S5 et S6.

Ensembles de données sur le neuroblastome et l’adénocarcinome pulmonaire

Nous avons utilisé un ensemble de données sur le neuroblastome de 176 échantillons (un sous-ensemble d’un ensemble de données plus vaste de 498 échantillons; appelé SEQC -neuroblastome dans ce manuscrit) pour évaluer la performance des pipelines d’ARN-seq en termes de prédiction basée sur l’expression génique de l’issue de la maladie. Ces échantillons ont été fournis par l’Hôpital Universitaire pour enfants de Cologne et séquencés au BGI à l’aide de la plateforme Illumina 48. Les 176 échantillons ont été prélevés sur des patients à haut risque qui étaient définis comme étant atteints d’un neuroblastome de stade 4 et d’âge > 18 mois ou présentant des tumeurs amplifiées par MYCN de tout stade ou âge. L’ensemble de données SEQC-neuroblastome a été déposé sur l’Expression génique Omnibus avec le numéro d’accession GSE47792.

Nous avons prédit deux paramètres cliniques : la survie sans événement (EFS), c’est-à-dire la survenue d’événements tels que la progression, la rechute ou le décès, et la survie globale (OS), c’est-à-dire la mort. Pour les deux critères d’évaluation, les patients ont été divisés en deux groupes (c.-à-d. risques élevés par rapport aux risques faibles). Les patients à haut risque ont connu un événement ou sont décédés avant un seuil de survie prédéfini, tandis que les patients à faible risque ont connu un événement ou sont décédés après le seuil, ou leur dernier suivi a dépassé le seuil. Les seuils de survie pour l’EFS et l’OS étaient respectivement de deux et trois ans. Les seuils ont été choisis pour équilibrer le nombre de patients à haut risque et de patients à faible risque. Les détails de l’ensemble de données SEQC-neuroblastome sont fournis dans le tableau supplémentaire S9.

Nous avons également utilisé un jeu de données ARN-seq de l’adénocarcinome pulmonaire de 87 échantillons provenant du dépôt de l’Atlas du génome du Cancer (TCGA). Le critère de prédiction était également la survie, et nous avons utilisé les mêmes critères pour définir les groupes à risque élevé et à faible risque avec le seuil de survie de deux ans. Le seuil de deux ans a été choisi pour équilibrer le nombre de patients à risque élevé et à faible risque. Les détails de l’ensemble de données TCGA-adénocarcinome pulmonaire sont fournis dans le tableau supplémentaire S10.

Filtrage de l’ensemble de données de référence qPCR pour produire un ensemble de gènes de référence

En raison de la variabilité des mesures qPCR et des désaccords entre les plateformes qpcr7, nous avons filtré l’ensemble de données de référence qPCR pour conserver les gènes présentant un comportement « correct”. Nous avons ensuite utilisé ces gènes pour calculer les métriques de référence (c’est-à-dire l’exactitude, la précision, la fiabilité et la reproductibilité). Un tel procédé de filtrage est résumé dans la Fig. L1.

En partant de l’ensemble initial de 20 801 gènes testés avec PrimePCR, nous avons filtré ces gènes pour ne retenir que les gènes quantifiés non nuls (c’est-à-dire détectés) et avec des valeurs Ct (seuil de cycle) ≤ 35 (35 indique la détection d’une seule molécule dans un échantillon). Le filtrage des données PrimePCR a abouti à 14 014 gènes qui correspondaient également au transcriptome AceView utilisé pour cartographier l’ensemble de données SEQC-benchmark RNA-seq.

Par la suite, nous avons filtré les 14 014 gènes qPCR pour ne retenir que 12 610 gènes présentant l’ordre de titrage correct (TO) et les rapports de mélange attendus (EMR). Les détails de ce processus se trouvent dans la section « Filtrage des gènes qPCR par ordre de titration et rapports de mélange attendus »”

Enfin, comme certaines métriques de référence telles que la précision et la précision sont sensibles aux gènes à expression nulle ou très faible, nous avons également sélectionné des gènes qui étaient exprimés comme non nuls dans toutes les répliques de tous les échantillons de tous les sites de séquençage et de tous les 278 pipelines d’ARN-seq. L’ensemble de référence final ne contient que 10 222 gènes qPCR (appelés « tous les gènes « ) qui ont été utilisés pour calculer les trois métriques de référence pour les pipelines ARN-seq.

D’après l’étude précédente, les gènes à expression plus faible sont plus susceptibles d’être incohérents entre les pipelines49. Ainsi, nous avons également identifié un ensemble de gènes à faible expression dans les 10 222 gènes sur la base de l’expression qPCR moyenne des échantillons A, B, C et D. Les 20% les plus bas des 10 222 gènes (c’est-à-dire 2044 gènes, appelés « gènes à faible expression”) ont également été utilisés pour calculer le même ensemble de mesures de référence pour les pipelines ARN-seq. Cette conception nous a permis d’étudier la capacité des pipelines ARN-seq à estimer l’expression de gènes à faible expression.

Filtrage des gènes qPCR par ordre de titrage et rapports de mélange attendus

Les ensembles de données SEQC-benchmark (RNA-seq et qPCR) ont des propriétés uniques qui permettent d’évaluer l’exactitude de la quantification. Après avoir identifié des gènes qPCR détectables (c’est-à-dire non nuls et Ct ≤ 35) et adaptés à AceView, nous avons utilisé deux métriques (TO et EMR) pour filtrer davantage l’ensemble de données qPCR de référence, ne laissant que des gènes qPCR « corrects”. Les métriques TO et EMR capturent des propriétés de mélange uniques des données, c’est-à-dire

CC= \frac{3}{4}A +\frac{1}{4}B\, \text{ and}\, =\frac{1}{4}A +\frac{3}{4}B.$$

En raison de cette propriété, tous les gènes doivent être exprimés dans l’un des ordres suivants, en fonction de l’expression relative des échantillons A et B :

AA\ge C\ge D\ge B\, \text {or}\, A\le C\le D\le B.$$
{{\stackrel{-}{q}} _{s, \cdot , k} =\frac{1}{N}\sum_{n = 1}^{N}{q}_{s, n, k, }$$

l’ensemble des gènes qPCR qui suivent l’ordre de titrage correct est

Pour un seul ensemble de données qPCR répliquées (par ex., l’ensemble de données PrimePCR que nous avons analysé), la variabilité inhérente d’une seule mesure qPCR peut entraîner certains gènes faux négatifs qui suivent le correct mais ne parviennent pas à être identifiés. À partir de la littérature50, 51, le coefficient de variation pour les mesures qPCR répliquées est généralement de 15% ou plus, nous avons donc utilisé ce nombre pour ajuster la marge permettant de déterminer si un gène suit le bon TO. Mathématiquement, nous avons calculé la plage de plus et moins un écart-type par rapport à chaque mesure qPCR et l’avons utilisée comme marge. Les équations révisées pour \({K}_{TO}\) sont les suivantes:

$${K}_{TO} = {K}_{TO, A\ge B}\cup {K}_{TO, A\le B, } where

où \(a= 1,15, b= 0,85\)

Outre TO, les échantillons doivent en outre présenter un rapport de mélange spécifique. Étant donné que le ratio entre les échantillons A et B est:

$${R}_{A,B}=\frac{A}{B}$$

le DME entre les échantillons C et D

$$EM{R}_{C,D}=\frac{3z\cdot {R}_{A,B}+1}{z\cdot {R}_{A,B}+3}\cdot \frac{z+3}{3z+1}$$
$${R}_{A,B}\in \left\equiv ,$$
$${R}_{C,D}\in \left\equiv \left,\text{ et}$$
$$EM{R}_{C,D}\in \left\equiv ),$$

et, enfin, détermine un ensemble de gènes qui satisfait aux EMR critère comme suit:

$${K}_{EMR}= \left\{k|\left({{R}_{C, D}^{Inférieur}\le{EMR}_{C, D}^{Supérieur}|}_{{k, R}_{C, D}\ge EM{R}_{C, D}}\right)\vee\left({{R}_{C, D}^{Supérieur}\ge {EMR}_{C, D}^{ Lower}|}_{{k, R}_{C, D}\le EM{R}_{C, D}}\right)\right\}$$

Pipelines d’analyse de données RNA-seq — cartographie, quantification et normalisation

Nous avons étudié 278 pipelines RNA-seq comprenant treize algorithmes de mappage de séquences18,19,20,21,22,23,24,25,26,27,28,29, trois catégories d’algorithmes de quantification d’expression31,32, 33 et sept méthodes de normalisation d’expression. Les tableaux supplémentaires S2 à S4 résument toutes les options envisagées pour chaque composant de pipeline (mappage de séquence, quantification d’expression et normalisation d’expression). Les treize algorithmes de mappage étudiés sont Bowtie18, Bowtie219, BWA20, GSNAP21, Magic22 (un nouveau pipeline développé par NCBI pour le projet SEQC : ftp://ftp.ncbi.nlm.nih.gov/repository/acedb/Software/Magic ), MapSplice23, Novoalign (un paquet commercialisé développé par Novocraft : https://www.novocraft.com/products/novoalign/), OSA24, RUM25, STAR26, Subread27, TopHat28 et WHAM29. Certains utilisent un mappage non épissé des lectures sur le transcriptome, et d’autres effectuent un mappage épissé sur le génome. Magic utilise les deux en parallèle et compare la qualité de chaque alignement pour garder le meilleur sur plusieurs cibles. Les algorithmes de mappage peuvent signaler uniquement un mappage unique ou autoriser plusieurs emplacements de mappage par lecture. Les algorithmes de quantification comprennent des méthodes simples basées sur le comptage (c.-à-d. HTSeq31) et des méthodes probabilistes basées sur la distribution de Poisson appliquées aux données de cartographie génomiques (c.-à-d. Cufflinks32) ou transcriptomiques (c.-à-d. RSEM33). La Magie, le RHUM et la Sous-lecture (c.-à-d., featureCounts52) les pipelines comprennent des méthodes de quantification intégrées qui entrent dans la catégorie des méthodes simples basées sur le comptage. Les méthodes de normalisation comprennent des méthodes de mise à l’échelle simples (fragments par million de fragments cartographiés, fragments par kilobase de longueur de gène par million de fragments cartographiés, quartile médian et supérieur), des méthodes de mise à l’échelle robustes (expression logarithmique relative et moyenne rognée des valeurs m) et des méthodes intégrées dans des pipelines spécifiques (indice d’expression magique).

Mappage de séquences

Nous avons mappé des séquences à chaque référence par étapes successives en utilisant des algorithmes de mappage non épissés ou épissés. Le mappage non épissé fait référence à des algorithmes qui alignent des séquences de lecture entières (par exemple, Bowtie2, BWA et Novoalign), tandis que le mappage épissé fait référence à des algorithmes qui divisent les lectures en segments pour tenir compte de longs espaces ou introns dans une lecture (par exemple, TopHat et MapSplice). Dans la première étape de la cartographie non épissée, nous avons tenté de mapper toutes les séquences d’extrémité appariées à la référence ERCC/MT/ARNr (i.e., L’ARN externe contrôle les séquences du Consortium, le génome mitochondrial et les séquences d’ARN ribosomiques). Toutes les paires de lecture non mappées ont ensuite été mappées sur le transcriptome AceView. Enfin, toutes les paires de lecture qui ne correspondaient pas aux références ERCC/MT/ARNr ou AceView ont été mappées à la référence du génome humain. Les coordonnées de cartographie transcriptomique ont ensuite été traduites en coordonnées de cartographie génomique et fusionnées avec les résultats de cartographie du génome humain pour produire les résultats finaux (fig. S21, panneau de gauche). Nous avons utilisé Bowtie2 comme mappeur pour la première étape de tous les pipelines de mappage épissés (Fig. S21, panneau de droite). Les algorithmes de mappage épissés ont soit directement mappé les lectures sur le génome humain (par exemple, MapSplice et GSNAP), soit mappé des lectures entières non épissées sur le transcriptome, puis fusionné ces résultats de mappage avec les résultats de mappage épissés des lectures restantes sur le génome humain (par exemple, TopHat et OSA). Le tableau supplémentaire S2 résume tous les outils de cartographie étudiés dans cette étude.

Bowtie2, GSNAP, Novoalign, TopHat et WHAM permettent de contrôler le nombre de mappages signalés par paire de lectures. Par défaut, ces algorithmes signalent généralement un seul meilleur emplacement de mappage par paire de lecture. Cependant, certains algorithmes de quantification peuvent utiliser des informations sur plusieurs emplacements de cartographie ambigus pour améliorer l’estimation de l’expression génique. Ainsi, en plus des rapports à un seul coup, nous avons également généré des résultats de mappage qui ont rapporté jusqu’à 200 coups par lecture (coups multiples). Nous avons également inclus le pipeline de cartographie Bowtie avec des paramètres de cartographie spécifiques pour la quantification avec RSEM, comme décrit dans la section 33 suivante.

Les options de ligne de commande pour tous les outils d’alignement des séquences sont détaillées dans la Note supplémentaire 1.

Quantification de l’expression génique

L’étape de quantification comprenait trois catégories de quantificateurs: les quantificateurs basés sur le comptage (c’est-à-dire HTSeq et les quantificateurs intégrés pour les pipelines Magic, RUM et Subread), les quantificateurs basés sur des modèles de probabilité pour la cartographie génomique (c’est-à-dire, Boutons de manchette), et des quantificateurs basés sur des modèles de probabilité pour la cartographie transcriptomique (c.-à-d., RSEM). Les principales caractéristiques de ces quantificateurs sont résumées dans le tableau supplémentaire S3. Cufflinks est un quantificateur basé sur un modèle de Poisson qui estime les probabilités d’affectation en lecture sur la base des informations d’alignement32. Il est capable à la fois d’assembler des transcriptions et de quantifier des expressions de gènes ou de transcriptions. Dans cette étude, nous avons désactivé la fonction d’assemblage et fourni le fichier GTF d’annotation du génome comme référence de quantification. HTSeq est un quantificateur naïf basé sur le comptage qui attribue des lectures mappées à genes31. HTSeq est capable de quantifier l’expression des gènes, mais pas l’expression de la transcription. RSEM est également un quantificateur basé sur un modèle de Poisson dont le concept est similaire à Cufflinks33. Les informations provenant des lectures à plusieurs coups sont importantes pour les boutons de manchette et les RSEM. Ces algorithmes utilisent des informations de lecture multi-coups pour estimer plus précisément l’expression d’un gène ou d’un transcrit.

Les résultats de cartographie des pipelines d’alignement n’étaient pas toujours compatibles avec les trois catégories de quantificateurs. Les boutons de manchette nécessitent que les résultats d’alignement soient triés par coordonnées d’alignement et que les lectures multi-coups soient marquées avec la balise ‘NH’ dans le champ d’attribut du fichier SAM. HTSeq exige que les résultats d’alignement soient triés par noms de lecture et que la balise ’NH’ soit absente du fichier SAM. RSEM ne quantifie que la cartographie transcriptomique, c’est-à-dire les lectures mappées et rapportées en coordonnées transcriptomiques. De plus, RSEM ne gère que les alignements non écartés. Ainsi, un filtrage est nécessaire pour supprimer les alignements espacés. En raison de ces exigences, nous avons pré-traité tous les résultats d’alignement avant la quantification. En résumé, vingt pipelines d’alignement, y compris des pipelines épissés, non épissés, à un seul coup et à plusieurs coups, convenaient à la quantification basée sur le comptage. Seize pipelines d’alignement convenaient aux boutons de manchette et seulement dix convenaient au RSEM. RSEM est spécialement conçu pour bien fonctionner avec Bowtie. Ainsi, nous avons également inclus ce pipeline de cartographie et de quantification intégré.

Les options de ligne de commande pour tous les outils de quantification sont détaillées dans la Note supplémentaire 1.

Normalisation de l’expression génique

La normalisation des données ARN-seq permet la comparaison entre échantillons. Généralement, les méthodes de normalisation corrigent la taille de la bibliothèque (c’est-à-dire le nombre total de lectures dans un échantillon), qui est la principale source de variance entre échantillons. Nous avons étudié sept méthodes de normalisation – fragments par million de fragments cartographiés (FPM), fragments par kilobase de longueur de gène par million de fragments cartographiés (FPKM), médiane (Med.), quartile supérieur (UQ), expression log relative (RLE), moyenne découpée des valeurs M (TMM) et indice d’expression (EIndex, qui est spécifique au pipeline magique) (voir Tableau supplémentaire S4). Nous décrivons chacune de ces méthodes de normalisation sur la base de la description mathématique suivante de l’ensemble de données SEQC-benchmark.

$$\overline{x}_{s,\cdot, k} = \frac{1}{N}\mathop\sum\limits_{n= 1}^{N}x_{s, n,k}we

nous avons défini l’ensemble des gènes présents comme étant

et l’ensemble final des gènes présents est

KK_{p} =K_{p , BGI}\cap K_ {p, MAI}.$$

Nous avons utilisé le même ensemble de gens présents pour toutes les méthodes de normalisation pour un pipeline RNA-seq.

Le nombre total de gènes présents pour un échantillon donné s et n répliqués est

$$x_{s, n}=\mathop\sum\limits_{{k\in K_{p}}} x_{s, n,k}, ,

et le nombre total moyen de gènes présents pour toutes les données d’un seul site est

$$\bar{x}=\frac {1}{4}\frac {1}{N}\mathop\sum\limits_{s}\mathop\sum\limits_{{n= 1}}^{N}x_{{s,n}}.$$

Ainsi, nous avons défini l’expression normalisée FPM pour chaque échantillon s, réplication n et gène k comme

$$y_{s, n, k}^{FPM}= \frac{{x_{s, n, k}\cdot\overline{x}}}{{x_{s, n}}}.expression

L’expression normalisée du quartile médian et supérieur pour chaque échantillon s, la réplication n et le gène k sont alors définies comme

$$y_{s, n,k}^{Med}=\frac{{x_{s, n,k}\cdot\tilde{x}}}{{\tilde{x}_{s, n}}}{\text{ et }}y_{s, n,k}^{ UQ} = \frac {{x_{s, n, k}\cdot\chapeau {x}}} {{\chapeau{x}_{s, n}}}.{}$$

Pour la normalisation FPKM, nous avons défini la longueur d’un gène k comme \(\ell_{k}\), qui est la longueur de l’union de tous les exons liés au gène tels que définis par le transcriptome d’AceView. La formulation originale de FPKM utilisait arbitrairement des facteurs d’échelle de 1 × 103 pour la longueur du gène et de 1 × 106 pour le nombre total de fragments cartographiés. Afin de maintenir une plage dynamique comparable entre toutes les méthodes de normalisation, nous avons plutôt mis à l’échelle la longueur moyenne des gènes et le nombre total moyen de tous les gènes présents. La longueur moyenne de tous les gènes présents est

$$\overline{\ell}=\frac{1}{{\left|{K_{p}}\right/|}\mathop\sum\limits_{{k\in K_{p}}}\ell_{k}.$$

Ainsi, l’expression normalisée FPKM redimensionnée pour chaque échantillon s, réplication n et gène k est

$$y_{s, n, k}^{FPKM}= \frac{{x_{s, n, k}\cdot\overline{\ell}\cdot\overline{x}}}{{x_{s, n}\cdot\ell_{k}}}.p

Les méthodes de normalisation TMM et RLE sont similaires à la normalisation FPM, mais introduisent un facteur d’échelle supplémentaire pour ajuster la taille de la bibliothèque. Nous avons utilisé le paquet edgeR dans R pour estimer un facteur d’échelle pour chaque réplique d’échantillon36,53. La méthode TMM sélectionne une bibliothèque de référence à partir d’un pool de bibliothèques de répliques d’échantillons, puis calcule des rapports d’expression logarithmique génique (valeurs M) et des valeurs d’expression logarithmique moyenne génique (valeurs A) entre la bibliothèque cible et la bibliothèque de référence. Les nombres extrêmes dans les valeurs M et les valeurs A sont rognés, et le facteur d’échelle de la bibliothèque cible est la moyenne pondérée des valeurs M restantes. La méthode RLE détermine un facteur de mise à l’échelle en définissant d’abord la bibliothèque médiane comme la moyenne géométrique génique des réplicats d’échantillon35. Le rapport médian de chaque bibliothèque cible à la bibliothèque médiane est pris comme facteur d’échelle. L’expression normalisée TMM et RLE pour chaque échantillon s, répliquée n et gène k est alors définie comme suit :

où \(\hat{f}_{s,n}^{TMM}\) et \(\hat{f}_{s,n}^{RLE}\) sont le facteur d’échelle pour l’échantillon s, répliquée n.

Les métriques de performance du pipeline RNA-seq

Les métriques de référence pour les pipelines RNA-seq sont résumées dans le tableau supplémentaire S7.

Précision mesurée en écart par rapport aux références qPCR

{{\stackrel{-}{y}}_{s,\cdot, k} = \frac{1}{N}\sum_{n =1}^{N}{y}_{s,n,k}$$

Compte tenu des échantillons A et B, l’écart absolu du rapport log de l’expression basée sur l’ARN seq par rapport à l’expression basée sur qPCR l’expression d’un gène k est

$$\Delta_{\frac{A}{B}, k} = \left|\log_2\left(\frac{\bar{x}_{A, ., k }} {\barre {x} _ {B, ., k }}\ droite) – \log_2\gauche (\frac{\bar{q}_{A, ., k }} {\bar {q} _ {B, ., k}}\right)\right |,

et la métrique de précision finale a été définie comme la médiane de tous les \({\Delta}_{{\frac{A}{B}, k}}\), \(k = 1\ldots K\).

Précision mesurée en tant que variation de l’expression des gènes dans les bibliothèques de réplication

Nous avons calculé le coefficient de variation (CoV) pour chaque gène et chaque échantillon dans quatre bibliothèques de réplication comme suit:

CoVCoV_{s, k} = \frac{{sd\left({x_{s,\cdot, k}}\right)}} {{\overline{x}_{s,\cdot, k}}},,

Fiabilité mesurée comme corrélation intra-échantillon de l’expression génique

La fiabilité d’un système de mesure peut être évaluée par le coefficient de corrélation intraclasse (ICC) 54,55. ICC est applicable aux mesures qui peuvent être organisées en groupes, et il décrit comment les mesures du même groupe sont similaires les unes aux autres. La définition moderne de l’ICC emprunte le cadre d’analyse de la variance (ANOVA), ou plus précisément d’ANOVA à effets aléatoires 55. Le type d’ANOVA dépend du plan expérimental et suit généralement la définition de l’article de Shrout publié en 197955. ICC(1,1) et ICC(1, k) sont basés sur le modèle à effets aléatoires unidirectionnels et s’appliquent au cas où chaque groupe est évalué par un ensemble différent de k évaluateurs sélectionnés au hasard parmi une population plus importante de évaluateurs. ICC(2,1) et ICC(2,k) sont basés sur le modèle des effets aléatoires bidirectionnels et s’appliquent au cas où un échantillon aléatoire de k évaluateurs est présélectionné dans une population plus grande et que chaque évaluateur évalue chaque groupe exactement une fois (c.-à-d., chaque évaluateur évalue n groupes au total). ICC(3,1) et ICC(3,k) sont basés sur le modèle des effets mixtes bidirectionnels et s’appliquent au cas où chaque groupe est évalué par chacun des mêmes évaluateurs k, qui sont les seuls évaluateurs de la population. Le deuxième paramètre en ICC (,) indique si l’ICC doit mesurer la fiabilité d’une seule mesure ou la moyenne de k mesures.

Pour l’ensemble de données de référence SEQC avec des bibliothèques de répliques pour chaque échantillon, ICC(1,1) ou ICC(1,k) correspondait à notre objectif car, pour un gène g, l’expression génique des bibliothèques de répliques pour différents échantillons (ou différents groupes dans le contexte précédent) n’a pas été évaluée exactement dans les mêmes conditions (ou évaluée par les mêmes évaluateurs dans le contexte précédent). Nous avons choisi d’utiliser ICC(1, k) car des bibliothèques de réplication sont disponibles pour la plupart des expériences. Mathématiquement, un modèle d’effets aléatoires unidirectionnels peut être formulé comme suit :

YY_{ij} =\mu +\alpha_{j} +{\varepsilon_{ij}}, ,
$$ICC\left({1, k}\right) =\frac{BMS-WMS}{{BMS}},We

Nous avons calculé ICC pour chaque gène k, \ (k = 1 \ ldots K \), puis a utilisé la médiane de tous les ICC comme mesure finale de fiabilité.

Nous avons également étudié d’autres métriques potentielles, telles que la reproductibilité, qui est définie comme la corrélation de Spearman entre deux bibliothèques répliquées d’un même échantillon (Note supplémentaire 2). La corrélation de Spearman variait de 0.993 à 0,996 (Fig. S8) en utilisant AllGenes. Nous avons écarté la mesure de reproductibilité en raison de la plage dynamique relativement petite.

Évaluation de l’utilité des métriques de référence pour la sélection des pipelines RNA-Seq

Nous avons classé les pipelines RNA-seq en fonction du rang moyen des trois métriques de référence (c’est-à-dire la précision, la précision et la fiabilité). Nous avons ensuite évalué l’utilité des mesures de référence en examinant si les pipelines de bons et de mauvais résultats identifiés sur la base des mesures de référence étaient informatifs pour déduire la performance de la prédiction basée sur l’expression génique de l’issue de la maladie et la signification statistique de la stratification des patients pour tous les paramètres cliniques (c.-à-d. les paramètres SEQC- neuroblastome EFS et OS et le paramètre de survie TCGA-poumon-adénocarcinome).

Tout d’abord, pour les 278 pipelines RNA-seq représentatifs appliqués à l’ensemble de données SEQC-benchmark, nous avons calculé le rang moyen en utilisant un sous-ensemble de mesures de référence comme indicateur de performance final pour chaque pipeline. Au total, nous avions 6 métriques (3 métriques de référence × 2 ensembles de gènes) et nous avons étudié 12 sous—ensembles (4 × 3) des 6 métriques en utilisant les critères suivants:

  1. (1)

    Quatre combinaisons des trois métriques de référence avec au moins deux dans un sous-ensemble – une combinaison avec les trois métriques de référence, trois combinaisons avec deux des trois métriques de référence.

  2. (2)

    Trois sous-ensembles formés par des métriques dérivées de tous les gènes, celles dérivées de gènes à faible expression, ou une combinaison des deux.

Deuxièmement, pour chacun des 278 pipelines RNA-seq représentatifs (156 pour le critère de survie TCGA-adénocarcinome pulmonaire), nous avons calculé l’ASC et la MCC de validation croisée imbriquées, comme décrit dans la section ”Méthode  » ” Modélisation prédictive du neuroblastome et de l’adénocarcinome pulmonaire », ce qui donne 834 (468 pour le critère de survie TCGA-adénocarcinome pulmonaire) valeurs de l’ASC et de la MCC pour chaque critère clinique ( c.-à-d., 278 pipelines × 3 classificateurs, ou 156 pipelines × 3 classificateurs) (Tableaux supplémentaires S11, S12). Nous avons également modélisé les fonctions de survie en utilisant l’analyse de Kaplan-Meier pour chaque pipeline, comme décrit dans la section ”Méthode » ”Analyse de survie de Kaplan–Meier ». Pour chaque pipeline d’ARN-seq, nous avons résumé la performance de la prédiction basée sur l’expression génique de l’issue de la maladie en utilisant à la fois l’ASC moyenne et la MCC entre les classificateurs et le taux de réussite de la stratification des patients (c.-à-d., séparation statistiquement significative de deux courbes de Kaplan-Meier) à travers toutes les itérations et classificateurs dans le cadre de validation croisée imbriquée.

Enfin, nous avons identifié les 10 % de pipelines les plus performants et les 10 % les moins performants sur la base du classement moyen d’un sous-ensemble des trois indicateurs de référence. La performance de prédiction correspondante (i.e., AUC et MCC) des pipelines performants ont été testés par rapport à ceux des pipelines peu performants en utilisant le test unilatéral de somme de rang de Wilcoxon avec l’hypothèse nulle que la médiane du premier groupe n’était pas plus grande que celle du dernier groupe.

Modélisation prédictive du neuroblastome et de l’adénocarcinome pulmonaire

Nous avons évalué la performance de 278 pipelines ARN-seq en termes de prise de décision basée sur l’expression génique à l’aide de l’ensemble de données SEQC-neuroblastoma48. L’ensemble de données SEQC-neuroblastome et les paramètres cliniques associés sont résumés dans le tableau supplémentaire S9. Les pipelines RNA-seq ont été évalués en termes de prédiction des résultats des patients atteints de neuroblastome pour deux paramètres cliniques à l’aide d’une validation croisée imbriquée (Fig. supplémentaire. S13) 56,57. Nous avons également évalué de la même manière la performance de 156 pipelines ARN-seq appliqués à l’ensemble de données TCGA-poumon-adénocarcinome pour prédire l’issue de la maladie. L’ensemble de données TCGA-adénocarcinome pulmonaire et le paramètre clinique associé sont résumés dans le tableau supplémentaire S10.

La validation croisée imbriquée implique la formation et la mise à l’essai d’un modèle de prédiction optimal. Ceci est accompli en utilisant l’optimisation triple ou la validation croisée interne, appliquée au sous-ensemble d’entraînement à partir de la validation croisée externe quintuple. Une fois que les paramètres finaux du modèle de prédiction optimal (c’est-à-dire les hyperparamètres du classificateur et la taille de la caractéristique) sont identifiés, le modèle final est entraîné en utilisant l’ensemble du sous-ensemble d’apprentissage, puis testé en utilisant le pli restant de la validation croisée externe quintuplée. Ce processus a été répété pendant dix itérations. Nous avons effectué une validation croisée imbriquée séparément pour chacun des trois classificateurs (c.-à-d., amplification adaptative, régression logistique et machines vectorielles de support) et a utilisé la méthode de sélection de fonctionnalités de redondance minimale et de pertinence maximale (mRMR) pour choisir des tailles de fonctionnalités optimales dans la plage de 5 à 40 avec une taille de pas de 558.

Analyse de survie de Kaplan–Meier

Pour chaque pipeline et classificateur RNA-seq (c’est-à-dire 278 pipelines × 3 classificateurs pour les paramètres SEQC-neuroblastome et 156 pipelines × 3 classificateurs pour le paramètre de survie TCGA–poumon-adénocarcinome), nous avons modélisé les fonctions de survie de Kaplan-Meier en fonction des étiquettes prédites de chaque échantillon. Nous avons ensuite utilisé le test de log-rank à deux queues pour déterminer si les courbes de survie estimées pour chaque groupe de patients prédits étaient statistiquement différentes.

Analyse de la variance et calcul de la contribution de chaque facteur de pipeline ARN-seq à la variance globale du pipeline

Nous avons utilisé l’analyse de la variance (ANOVA) pour déterminer si chaque facteur de pipeline ARN-seq contribue de manière significative à la variance de chacune des trois mesures de référence (c.-à-d. précision, précision et fiabilité) ainsi qu’à la variance des performances de prédiction (c.-à-d. ASC et MCC). Pour chacune des trois métriques de référence, nous avons utilisé un modèle linéaire (fonction R « lm”) pour ajuster les données des 278 pipelines en utilisant la métrique comme variable dépendante et les facteurs de pipeline RNA-seq comme variables catégorielles indépendantes. Nous avons considéré les facteurs suivants comme des variables catégorielles indépendantes : algorithme de mappage, stratégie de mappage (c’est-à-dire épissé par rapport à non épissé), rapport de mappage (c’est-à-dire à un seul coup par rapport à plusieurs coups), algorithme de quantification et algorithme de normalisation. Nous avons inclus tous les facteurs et leurs interactions bidirectionnelles dans le modèle linéaire. Pour chacun des paramètres de prédiction, nous avons appliqué la même technique pour ajuster les données des 278 pipelines en utilisant l’ASC moyenne ou la MCC comme variable dépendante et le même ensemble de facteurs de pipeline RNA-seq comme variables catégorielles indépendantes. Nous avons ensuite effectué l’ANOVA sur le modèle linéaire (fonction R « anova »). ANOVA calcule une ”somme de carrés » (c’est-à-dire une variance) attribuée à chaque facteur ou interaction et utilise un test F pour déterminer si la variance est statistiquement significative. Nous avons calculé le pourcentage que chaque facteur ou interaction contribue à la variance totale en calculant le rapport de la « somme des carrés” pour chaque facteur à la somme totale des carrés.

Analyse de régression

Nous avons étudié la relation entre les profils d’alignement ou les caractéristiques de distribution de l’expression génique et les métriques de référence. Les profils d’alignement comprenaient le nombre total de fragments mappés, le nombre total de lectures couvrant la région intronique, le nombre total de lectures avec insertions ou suppressions, le nombre total de lectures parfaitement appariées, le nombre total de lectures avec au plus une non-concordance et le nombre de non-concordances par lecture mappée. Chaque algorithme d’alignement était représenté par les statistiques moyennes sur 2 sites de séquençage, 4 échantillons, 4 bibliothèques de réplication et 2 voies. En utilisant le package ”MASSE » en R, nous avons adopté l’approche de pondération de Huber pour adapter des modèles de régression linéaire robustes entre une variable dépendante (performance métrique de référence) et une variable explicative (profil d’alignement). L’estimation M avec l’approche de pondération de Huber est une méthode de régression robuste en présence de valeurs aberrantes. Les caractéristiques de distribution de l’expression génique comprenaient le quartile inférieur, le quartile médian, le quartile supérieur, le maximum, la plage interquartile, l’écart-type, l’asymétrie, le kurtose et l’entropie d’une distribution d’expression génique. Nous avons utilisé la même méthode d’estimation M avec pondération de Huber pour adapter des modèles de régression linéaire robustes, puis nous avons rapporté l’erreur-type résiduelle pour chaque modèle.

Avertissement

Les points de vue présentés dans cet article ne reflètent pas nécessairement l’opinion ou la politique actuelle ou future de la Food and Drug Administration des États-Unis. Toute mention de produits commerciaux est faite à titre de clarification et non comme une approbation.

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *