Informations

Quelle est la meilleure méthode pour vérifier l'expression différentielle des gènes dans un grand ensemble de données ? Je veux identifier les gènes du vieillissement exprimés de manière différentielle


Mes données sont des données RNAseq et je veux savoir quels gènes (que je recherche) sont exprimés de manière différentielle avec l'âge. J'ai essayé le test de Kruskal-Wallis et les méthodes de changement de pli. Avez-vous de meilleures suggestions?


Un pipeline standard pour DGE serait Salmon pour la (pseudo)cartographie + DESeq2 pour l'analyse statistique.

Salmon fait partie d'un ensemble de logiciels de cartographie modernes, rapides et précis. Nécessite la définition d'un transcriptome.

DESeq2 est un package R (bioconducteur) mature qui peut gérer des conceptions raisonnablement complexes, y compris des séries chronologiques, mais pas des modèles mixtes. Les vignettes sont assez complètes. Il peut gérer des données provenant d'un logiciel de pseudo-mapping ainsi que des données générées à l'aide de méthodes plus traditionnelles, par ex. Star + FeatureCounts.

L'utilisation d'un package statistique RNAseq prédéfini est une bien meilleure idée que d'essayer de le faire, par ex. un GLM Poisson/Neg.bin manuellement pour chaque gène (même si c'est finalement ce que fait DESeq2).

(d'accord que cette question convient mieux à SE Bioinf)


Fond

Les analyses statistiques à grande échelle sont devenues la marque de fabrique de la recherche biologique post-génomique en raison des progrès des tests à haut débit et de l'intégration de grandes bases de données biologiques. Un problème associé est l'estimation simultanée des valeurs p pour un grand nombre de tests d'hypothèses. Dans de nombreuses applications, une hypothèse paramétrique dans la distribution nulle telle que la normalité peut être déraisonnable, et les valeurs p basées sur le rééchantillonnage sont la procédure préférée pour établir la signification statistique. L'utilisation de procédures basées sur le rééchantillonnage pour des tests multiples est gourmande en calculs et nécessite généralement un grand nombre de rééchantillonnages.

Résultats

Nous présentons une nouvelle approche pour affecter plus efficacement les rééchantillonnages (tels que les échantillons bootstrap ou les permutations) dans un cadre de tests multiples non paramétrique. Nous avons formulé une approche d'inspiration bayésienne de ce problème et conçu un algorithme qui adapte l'affectation des rééchantillonnages de manière itérative avec un espace et un temps d'exécution négligeables. Dans deux études expérimentales, un ensemble de données de puces à ADN pour le cancer du sein et un ensemble de données d'étude d'association à l'échelle du génome pour la maladie de Parkinson, nous avons démontré que notre procédure d'attribution différentielle est nettement plus précise que l'attribution de rééchantillonnage uniforme traditionnelle.

Conclusion

Nos expériences démontrent que l'utilisation d'une stratégie d'allocation plus sophistiquée peut améliorer notre inférence pour les tests d'hypothèses sans augmentation drastique de la quantité de calculs sur les données aléatoires. De plus, nous gagnons davantage en efficacité lorsque le nombre de tests est important. Le code R de notre algorithme et la méthode de raccourci sont disponibles sur http://people.pcbi.upenn.edu/


Données expérimentales

Les données utilisées dans ce flux de travail sont une expérience RNA-Seq de cellules musculaires lisses des voies respiratoires traitées avec de la dexaméthasone, un stéroïde glucocorticoïde synthétique ayant des effets anti-inflammatoires. Les glucocorticoïdes sont utilisés, par exemple, chez les patients asthmatiques pour prévenir ou réduire l'inflammation des voies respiratoires. Dans l'expérience, quatre lignées cellulaires primaires de muscles lisses des voies respiratoires humaines ont été traitées avec 1 micromolaire de dexaméthasone pendant 18 heures. Pour chacune des quatre lignées cellulaires, nous avons un échantillon traité et un échantillon non traité. La référence de l'expérience est :

Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q &ldquoRNA-Seq Le profilage du transcriptome identifie CRISPLD2 comme un gène sensible aux glucocorticoïdes qui module la fonction des cytokines dans les cellules musculaires lisses des voies respiratoires.&rdquo PLoS One. 2014 juin 139(6) : e99625. PMID : 24926665. GEO : GSE52778.


Matériaux et méthodes

Données de séquençage d'ARN

Les données anonymisées de niveau 3 pour les échantillons de cancer de la prostate et tous les échantillons non malins disponibles de ces patients atteints de cancer de la prostate ont été téléchargées à partir du portail de données The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov). Le niveau 3 décrit les données qui ont été traitées et agrégées pour donner des signaux d'expression génique pour un échantillon. Pour chaque échantillon, les données contiennent le nombre d'expressions pour un maximum de 20 531 transcrits d'ARN codant et non codant, ainsi que des informations cliniques telles que l'âge, le stade, le score de Gleason, le niveau de PSA et la race/l'origine ethnique. Avant l'analyse, les échantillons tumoraux et non malins ont été tirés au hasard pour obtenir un pool de 225 échantillons appariés selon l'âge et le stade (tableau S1). Un total de 173 échantillons de cancer de la prostate et 52 échantillons non malins provenant de 204 patients uniques ont été analysés. Les informations cliniques du patient sont présentées dans le tableau 1.

Expression génétique différentielle

L'environnement de programmation R (version 3.1.2) [61] a été utilisé pour traiter les données brutes, effectuer des calculs statistiques et effectuer une analyse d'expression différentielle. Après appariement de l'âge et du stade, 393 transcrits ont été supprimés car ils manquaient d'expression dans les 225 échantillons composant l'ensemble de données. Les comptes d'ARN pour les 20 138 transcrits restants ont été arrondis au nombre entier le plus proche et compilés dans une matrice pour créer l'ensemble de données. L'amplitude des changements d'expression par rapport aux échantillons non malins a également été calculée en prenant le logarithme en base 2 du rapport d'expression moyenne tumeur/non malin. Pour les gènes sans expression dans les échantillons tumoraux ou non malins, le log2 les changements de pli ont été ajustés en ajoutant un à chaque moyenne, puis en calculant le rapport. Tous les journaux2 les valeurs citées sont des valeurs après de tels ajustements. Des changements de pli négatifs indiquaient une régulation à la baisse dans les échantillons de tumeur alors que des valeurs positives indiquaient une régulation à la hausse. Le package R DESeq2 (version 1.6.3) [62] a été utilisé pour identifier les DEG dans les données d'ARN des patients TCGA. Le calcul a été effectué sur le cluster de calcul haute performance de la Florida State University. DESeq2 a renvoyé une valeur P déterminée par les statistiques de Wald et une valeur P ajustée (valeur Q) pour corriger les tests de comparaisons multiples à l'aide de la méthode Benjamini-Hochberg pour déterminer le taux de fausses découvertes (FDR). Les DEG ont été définis comme des gènes différents avec un FDR inférieur à 1 % (Q < 0,01).

Pour évaluer l'importance des DEG identifiés, des analyses ont été menées pour rechercher les voies surreprésentées, l'enrichissement de l'ensemble de gènes et l'impact des voies de signalisation. Dans un premier temps, des éléments surreprésentés ont été identifiés parmi les DEG. Le système de classification de l'analyse des protéines à travers les relations évolutives (PANTHER) et des outils d'analyse ont été utilisés pour catégoriser les DEG par classe de protéines PANTHER, fonction moléculaire de l'ontologie génique (GO) et processus biologique GO afin de déterminer si l'une de ces classes ou termes GO était surreprésenté [ 63]. Le test de surreprésentation PANTHER (version 20150430) a été utilisé pour rechercher les données dans la base de données PANTHER (version PANTHER 10.0 publiée le 2015-05-15) et la base de données GO (publiée le 2015-05-09) afin d'identifier les classes de protéines ou les annotations GO surreprésentées. dans nos données par rapport à un génome humain de référence. Les valeurs p ont été ajustées à l'aide d'une correction de Bonferroni.

Analyse de la voie

L'analyse d'enrichissement du jeu de gènes (GSEA) [64] a été utilisée pour identifier des groupes de gènes enrichis dans la tumeur ou dans l'état non malin. L'outil d'analyse GSEA (version 2.2.0) a été téléchargé à partir du site Web du Broad Institute (http://www.broadinstitute.org/gsea/index.jsp). Des ensembles de gènes sélectionnés des voies BioCarta et Reactome ont été téléchargés à partir de la base de données de signatures moléculaires du Broad Institute. Un ensemble de gènes supplémentaires a été construit à partir des voies de l'Encyclopédie des gènes et des génomes de Kyoto (KEGG) [65]. Les voies les moins pertinentes pour le cancer de la prostate ont été exclues. Les voies KEGG incluses dans l'analyse sont répertoriées dans les informations complémentaires (tableau S2). L'intégralité de la matrice de comptage d'expression d'ARN a été chargée dans l'application GSEA sans limiter l'entrée aux seuls DEG. Les petits (< 5 gènes) et les grands (> 500 gènes) ensembles de gènes ont été exclus de l'analyse.

L'analyse d'impact des voies de signalisation (SPIA) a été utilisée pour évaluer l'importance des voies enrichies en termes d'impact et de capacité à activer ou inhiber une voie [66]. L'analyse SPIA a été réalisée à l'aide du package R « SPIA » (version 2.18.0) [67]. Entrez les identifiants, connectez-vous2 les changements de pli et les valeurs Q pour tous les gènes ont été compilés. Le seuil d'expression différentielle utilisé dans l'algorithme SPIA était basé sur la valeur Q ajustée par FDR. L'analyse a été effectuée en utilisant la même liste personnalisée de voies que celle utilisée dans la GSEA (tableau S2) et des versions mises à jour de ces voies ont été téléchargées avant d'exécuter l'analyse (consulté le 29/07/2015).


Résultats

Méta-analyse de l'expression différentielle

Nous avons évalué les niveaux globaux d'expression génique dans les ensembles de données en attribuant des valeurs de classement à chaque gène en fonction de sa valeur d'expression moyenne dans un ensemble de données. Alors que nous observons des variations entre les études, il apparaît toujours un schéma clair de gènes qui sont constamment fortement ou faiblement exprimés dans le cerveau (données non présentées) soutenant la faisabilité de comparer les études les unes aux autres.

Nous avons utilisé la régression linéaire pour évaluer les changements d'expression génique par rapport à quatre facteurs (âge, sexe, pH cérébral et PMI), au sein de chaque ensemble de données. Nous avons pris en compte les deux directions de changement (régulation à la hausse et à la baisse) pour chaque facteur, créant jusqu'à huit scénarios différents pour chaque ensemble de données. Bien que l'objectif de cet article soit de rendre compte des résultats de la méta-analyse à travers des ensembles de données, nous avons brièvement résumé ici les résultats d'études individuelles. Le nombre de gènes dont la preuve change avec chacun des facteurs (q < 0,01) est indiqué dans le tableau 2 . Sans surprise, les ensembles de données avec une taille d'échantillon plus petite ont montré moins de changements statistiquement significatifs associés aux facteurs. Dans l'ensemble, les facteurs associés à l'expression différentielle la plus robuste étaient l'âge et le pH du cerveau. Les gènes qui changent dans les ensembles de données individuels mais pas à partir de la méta-analyse peuvent être trouvés dans le tableau supplémentaire 7 à http://chibi.ubc.ca/

Tableau 2

Gènes significatifs (qπ.01) identifiés dans chaque ensemble de données individuel

Base de donnéesÂgeSexepHPMI
Vers le basEn hautFemelleHommeVers le basEn hautVers le basEn haut
Stanley Chen00270000
Stanley Kato00110703000
<"type":"entrez-geo","attrs":<"text":"GSE1572","term_id":"1572">> GSE15721626416n / An / A00
<"type":"entrez-geo","attrs":<"text":"GSE2164","term_id":"2164">> GSE216400120000
<"type":"entrez-geo","attrs":<"text":"GSE8919","term_id":"8919">> GSE8919153133619n / An / A2278
<"type":"entrez-geo","attrs":<"text":"GSE11512","term_id":"11512">> GSE11512003131000
<"type":"entrez-geo","attrs":<"text":"GSE11882","term_id":"11882">> GSE118824280214n / An / A10
<"type":"entrez-geo","attrs":<"text":"GSE13162","term_id":"13162">> GSE131620000n / An / A00
<"type":"entrez-geo","attrs":<"text":"GSE5390","term_id":"5390">> GSE539011n / An / An / An / A00
McLeanPFC001500n / An / A
<> GSE379000113n / An / An / An / A
Signature ‘Union’689198103370402278
Chevauchement des méta-signatures4151026191910146

Les valeurs Q ont été calculées à partir des valeurs p de régression pour chaque facteur dans chaque ensemble de données. Le nombre de gènes rapportés ici est significatif à q < 0,01. La signature ‘union’ représente l'union de gènes uniques identifiés pour chaque facteur dans tous les ensembles de données. Le chevauchement de méta-signature indique le nombre de gènes de signature d'union chevauchant la méta-signature correspondante.

Pour examiner les changements dans l'expression des gènes qui étaient cohérents dans tous les ensembles de données, ou étayés par des preuves provenant de plusieurs ensembles de données, nous avons mis en œuvre une approche de méta-analyse croisée (voir Méthodes). Le résultat de cette analyse est de huit méta-signatures (à la hausse ou à la baisse pour chacun des quatre facteurs). Les dix principaux gènes de chaque méta-signature se trouvent dans le tableau 7 . Pour examiner les résultats, nous avons d'abord extrait les valeurs p correspondantes de chaque ensemble de données individuel et les avons visualisées sous forme de tracés (lissés) dans l'ordre déterminé par l'ordre des méta-signatures (Figure 1). Nous observons que les gènes qui ont de bonnes valeurs méta-q ont tendance à avoir de bonnes valeurs p dans plusieurs, mais pas nécessairement toutes les études. Nous démontrons également cela en utilisant une visualisation de carte thermique alternative (Figure supplémentaire 1). Des résultats plus détaillés sont tracés pour certains exemples de gènes sur la figure 2 . Ces graphiques montrent que les valeurs p pour un gène donné varient d'un ensemble de données à l'autre, et la méta-analyse peut identifier des gènes qui ne présentent que des effets faibles ou non significatifs dans certains ensembles de données. Par exemple, sur la figure 2, pour les gènes d'âge Gfap et Rgs4, nous constatons de faibles changements dans le niveau d'expression (respectivement à la hausse et à la baisse), qui ne sont pas significatifs après plusieurs corrections de tests dans ces études. D'autre part, nous trouvons également des gènes qui montrent des effets significatifs dans la plupart sinon toutes les études (par exemple, Xist, Figure 2). Dans la figure 3, nous avons rassemblé les 50 principaux gènes régulés à la baisse avec l'âge et tracé les niveaux d'expression dans chaque ensemble de données avec des échantillons classés par âge croissant. Pour la plupart des études, nous observons un gradient à travers l'ensemble de données à mesure que l'expression diminue de niveaux élevés à faibles, illustrant que la méta-analyse récupère de nombreux gènes qui montrent des tendances assez cohérentes entre les ensembles de données.

Pour chaque ensemble de données utilisé, les valeurs p des gènes ont été tracées par rapport à la valeur meta-q correspondante et un ajustement de loess a été calculé pour générer une courbe lisse entre les points. Le fait que la plupart des ensembles de données montrent une augmentation des valeurs p corrélées aux valeurs méta-q indique la contribution de signaux d'intensité variable aux méta-signatures. Une vue alternative des données utilisant des cartes thermiques est disponible dans le supplément. Les courbes déformées pour le sexe sont dues aux effets puissants d'un petit nombre de gènes avec de très petites valeurs méta-q (notez la différence d'échelle de D par rapport à A𠄼).

Pour les gènes sélectionnés de chacune des méta-signatures, nous avons tracé les valeurs p de régression logarithmique de chaque ensemble de données. Les cercles vides représentent les ensembles de données pour lesquels le gène s'est avéré significatif après plusieurs corrections de tests (q < 0,01). La ligne pointillée indique un niveau de signification de la valeur p par étude de 0,05 à titre de référence.

Les 50 principaux gènes régulés à la baisse par l'âge ont été sélectionnés en fonction du classement de la valeur q de la méta-analyse. Pour chaque gène, les données correspondantes de chaque étude ont été extraites et converties en une carte thermique. Les valeurs d'expression ont été normalisées entre les échantillons de chaque ensemble de données et classées par âge. L'âge est tracé en haut de chaque carte thermique. Les valeurs de lumière dans la carte thermique indiquent une expression plus élevée. Les barres grises indiquent les valeurs manquantes. Tous les ensembles de données sont approximativement à la même échelle horizontale, sauf le dernier, qui est compressé pour tenir sur la page.

Tableau 7

Âge régulé à la baisseÂge régulé à la hausse
OLFM1olfactomédine 1NEBLnébuleuse
KCNF1canal potassique voltage-dépendant, sous-famille F, membre 1MED12médiateur complexe sous-unité 12
RGS4régulateur de la signalisation des protéines G 4BCL2LLC à cellules B/lymphome 2
PPP3CBprotéine phosphatase 3 (anciennement 2B), sous-unité catalytique, isoforme bêtaGMPRguanosine monophosphate réductase
ADCY2adénylate cyclase 2 (cerveau)GFAPprotéine acide fibrillaire gliale
SVOPHomologue de la protéine apparentée à SV2 (rat)ADH1Balcool déshydrogénase 1B (classe I), polypeptide bêta
EFNB3éphrine-B3WWOXdomaine WW contenant de l'oxydoréductase
ATP2B2ATPase, transportant Ca++, membrane plasmique 2PLEC1plectine 1, protéine de liaison au filament intermédiaire 500kDa
HPCAhippocalcineVCANversique
CALB1calbindine 1, 28kDaAHCYL1S-adénosylhomocystéine hydrolase-like 1
pH régulé à la baissepH régulé à la hausse
FGF2facteur de croissance des fibroblastes 2 (basique)SLC1A1transporteur de soluté famille 1 (transporteur neuronal/épithélial de glutamate de haute affinité, système Xag), membre 1
AHCYL1S-adénosylhomocystéine hydrolase-like 1GRANDlike-glycosyltransférase
DTNAdystrobrevine, alphaC17orf81cadre de lecture ouvert du chromosome 17 81
MAPKAPK2protéine kinase activée par un mitogène protéine kinase 2 activéeHISPPD2Adomaine histidine acide phosphatase contenant 2A
TJP1protéine de jonction serrée 1 (zona occludens 1)PRKCDprotéine kinase C, delta
S100A13S100 protéine de liaison au calcium A13DLG3disques, grand homologue 3 (Drosophila)
RBBP6protéine de liaison au rétinoblastome 6KCNAB1canal potassique voltage-dépendant, sous-famille liée à l'agitateur, membre bêta 1
GNG12protéine de liaison aux nucléotides de guanine (protéine G), gamma 12SLC8A1famille de porteurs de solutés 8 (échangeur sodium/calcium), membre 1
ANP32Efamille des phosphoprotéines nucléaires acides (riches en leucine) 32, membre EGABRA5Récepteur de l'acide gamma-aminobutyrique (GABA) A, alpha 5
BAALCcerveau et leucémie aiguë, cytoplasmiqueRIT2Ras-like sans CAAX 2
Femelle régulée à la hausseMâle régulé à la hausse
XISTTranscrit spécifique de X (inactif) (non codant pour une protéine)JARID1Djumonji, domaine interactif riche en AT 1D
HDHD1Adomaine hydrolase de type haloacide déshalogénase contenant 1AUSP9Ypeptidase 9 spécifique à l'ubiquitine, liée à l'Y (de type facettes graisseuses, drosophile)
UTXrépétition tétratricopeptide transcrite de manière ubiquitaire, chromosome XEIF1AYfacteur d'initiation de la traduction eucaryote 1A, lié à l'Y
JARID1Cjumonji, domaine interactif riche AT 1CCYorf15Bcadre de lecture ouvert 15B du chromosome Y
TSIXARN antisens XIST (non codant pour les protéines)DDX3YPolypeptide 3 de la boîte DEAD (Asp-Glu-Ala-Asp), lié à l'Y
USP9Xpeptidase 9 spécifique de l'ubiquitine, liée à l'XUTYgène répété de tétratricopeptide transcrit de manière ubiquitaire, lié à l'Y
LOC554203domaine alanyl-ARNt synthétase contenant 1 pseudogèneRPS4Y1protéine ribosomique S4, liée à l'Y 1
STSstéroïde sulfatase (microsomal), isozyme STTTY15transcription spécifique aux testicules, liée à l'Y 15
ZFXprotéine à doigt de zinc, liée à l'XCYorf15Acadre de lecture ouvert du chromosome Y 15A
PNPLA4domaine phospholipase de type patatine contenant 4TMSB4Ythymosine bêta 4, liée à l'Y
PMI régulé à la baissePMI régulé à la hausse
BRD8bromodomaine contenant 8GOSR2membre du complexe récepteur golgi SNAP 2
RBM5Protéine de motif de liaison à l'ARN 5CYB5Bcytochrome b5 type B (membrane mitochondriale externe)
PUM2homologue pumilio 2 (Drosophila)SNF8SNF8, sous-unité du complexe ESCRT-II, homologue (S. cerevisiae)
ARHGEF7Facteur d'échange de nucléotides Rho guanine (GEF) 7GRLF1récepteur des glucocorticoïdes facteur de liaison à l'ADN 1
NCOA3coactivateur des récepteurs nucléaires 3C6orf1chromosome 6 cadre de lecture ouvert 1
DAPK1protéine kinase 1 associée à la mortMGMTO-6-méthylguanine-ADN méthyltransférase
SORL1récepteur lié à la sortieline, L (classe DLR) A contenant des répétitionsMAXIMUMFacteur X associé à MYC
KIAA0841KIAA0841ST3GAL2ST3 bêta-galactoside alpha-2,3-sialyltransférase 2
CAMK2Gprotéine kinase II gamma dépendante du calcium/calmodulinePTGER3récepteur 3 de la prostaglandine E (sous-type EP3)
DIP2CDIP2 disco-interacting protein 2 homologue C (Drosophila)EXT1exostoses (multiples) 1

Pour chaque méta-signature, nous avons répertorié les dix principaux gènes, classés par valeur méta-q (tous à q < 0,001). Pour chaque gène, nous avons répertorié le symbole du gène et un nom de gène. Les listes complètes de méta-signatures pour chaque facteur se trouvent dans le tableau supplémentaire 6 http://chibi.ubc.ca/

Alors que nous avons présenté les résultats d'une analyse qui a traité chaque facteur indépendamment (régression linéaire), nous avons également effectué une méta-analyse qui modélise les facteurs simultanément dans une analyse de covariance, donnant des résultats très similaires (disponible sur http://chibi.ubc .Californie/

ministère/). Nous avons également tenté de modéliser les interactions entre les facteurs, mais cela a été difficile en raison du manque de puissance.

Évaluation du méta-profil

Nous avons testé la robustesse de nos méta-profils en utilisant une procédure jackknife. Cela impliquait de supprimer séquentiellement un ensemble de données, d'effectuer la méta-analyse sur les ensembles de données restants, puis de sélectionner les gènes à un méta-q < 0,01. Cette procédure a été répétée pour chaque ensemble de données à tour de rôle, et les gènes trouvés dans tous les cycles ont été conservés comme signature 𠇌ore”. Chacune des signatures de base englobe plus de la moitié des gènes trouvés dans les méta-profils correspondants, à l'exception des méta-profils de pH. Les signatures 𠇜ore” se trouvent dans le tableau supplémentaire 6 à http://chibi.ubc.ca/

Les études que nous avons sélectionnées pour la méta-analyse n'étaient en général pas conçues pour tester les effets de l'âge, du sexe, du pH ou du PMI. ). Cependant, comme les échantillons de cerveau humain sont difficiles à obtenir, une variabilité suffisante est autorisée dans les études pour nous permettre d'effectuer une méta-analyse significative. Nous nous demandions toujours dans quelle mesure nos résultats concorderaient avec des études plus ciblées. Par conséquent, dans le but de valider notre approche, nous avons identifié des listes de gènes indépendantes de la littérature pour l'âge, le sexe et le pH cérébral (nous n'avons pas pu trouver d'ensemble de validation complet pour le PMI). Des détails sur les ensembles de validation peuvent être trouvés dans le tableau supplémentaire 2. Chaque liste de gènes de validation a ensuite été séparée en deux groupes en fonction de la direction du changement, pour correspondre aux méta-signatures obtenues à partir de notre méta-analyse. De toute évidence, aucune de ces listes de gènes de validation ne peut être considérée comme de véritables étalons-or, mais aide à mettre nos résultats dans le contexte des découvertes précédentes.

Pour quantifier le pouvoir prédictif de nos méta-signatures d'analyse par rapport aux ensembles de validation correspondants, nous avons d'abord effectué une analyse standard des caractéristiques de fonctionnement du récepteur (ROC). Le score rapporté pour chaque profil est l'aire sous la courbe ROC (AUC), une valeur comprise entre 0 et 1, où 1,0 correspond parfaitement à la liste externe et 0,5 refléterait un ordre aléatoire. Les valeurs AUC pour chacun des profils sont rapportées dans le tableau 3 . Des parcelles ROC individuelles peuvent être trouvées dans la figure 2 supplémentaire. Nous avons également testé l'effet de l'utilisation d'un seuil statistique spécifique pour sélectionner des gènes à partir des méta-profils, en collectant des gènes à deux niveaux de signification -q < 0,001). Le chevauchement avec l'ensemble de validation était significatif (pπ.001, test exact de Fisher Tableau 3 ) pour toutes les signatures. Nous avons également trouvé un chevauchement comparable entre chacune des signatures de base et les ensembles de validation.

Tableau 3

Comparaison des profils de méta-signature avec des ensembles de gènes de validation

Nombre de gènes de profil (q < 0,001)Chevauchement avec l'ensemble de validation (q < 0,001)Nombre de gènes de profil (q < 0,01)Chevauchement avec l'ensemble de validation (q < 0,01)Score ASC
Gènes d'âge
Régulé à la hausse404401241690.74
Régulé à la baisse113411322471360.76
Gènes du pH
Régulé à la hausse25113681310.91
Régulé à la baisse2155510181220.86
Gènes sexuels
Homme1473870.90
Femelle1911281n / A
Gènes PMI
Régulé à la hausse75n / A691n / An / A
Réglementé à la baisse4n / A49n / An / A

Les profils de pH du cerveau ont donné des scores AUC élevés de 0,91 et 0,86 (gènes régulés à la hausse et à la baisse, respectivement), lorsqu'ils ont été validés par rapport aux gènes sensibles au pH DLPFC tirés de Vawter et al. (2006). De plus, nous avons obtenu des scores raisonnablement élevés en utilisant une liste de gènes de pH indépendante plus petite obtenue de Mexal et al. (2006), malgré une différence dans la région cérébrale utilisée entre l'étude de validation et la méta-analyse (données non présentées). Étant donné que le pH lui-même varie probablement d'une région du cerveau à l'autre (Mexal et al., 2006), nos résultats sont cohérents avec l'hypothèse selon laquelle les changements liés au pH dans l'expression des gènes sont similaires dans toutes les régions du cerveau. Les profils d'âge, en revanche, présentaient des scores légèrement inférieurs à ceux obtenus pour le pH cérébral. Erraji-Benchekroun et al. ont utilisé des échantillons des zones Brodmann 9 (PFC dorsolatéral, BA9) et 47 (PFC orbitofrontal, BA47) de chaque sujet, pour évaluer les différences d'expression selon l'âge, montrant des changements comparables dans les deux régions du cerveau (Erraji-Benchekroun et al., 2005). En tant que tel, notre ensemble de validation se composait de gènes montrant des changements d'expression de l'âge collectivement dans les deux régions cérébrales néocorticales BA9 et BA47. Alors que bon nombre de ces gènes apparaissent en haut de nos listes de classement, certains sont également dispersés dans notre classement. La méta-signature de genre de notre méta-analyse a obtenu un score élevé lorsqu'elle a été validée avec un ensemble de gènes de Galfalvy et al. (2003), avec la plupart des gènes de validation en haut du classement.

Chevauchement avec les gènes candidats de la schizophrénie

Nous avons comparé des gènes significatifs (méta-q < 0,001) de chacune de nos méta-signatures avec des gènes connus pour être associés à la schizophrénie. Nous avons extrait une liste de 34 gènes candidats de la schizophrénie fournie dans une revue de littérature complète (Colantuoni et al., 2008), et avons recherché cette liste de gènes dans chacune de nos méta-signatures. Nous avons constaté que 12 de ces gènes s'identifiaient avec au moins une de nos méta-signatures, bien que la majorité des chevauchements aient été observés parmi les méta-profils d'âge et de pH (tableau 4). Le chevauchement des gènes avec chacune des méta-signatures d'âge était significatif à p < 0,01.

Tableau 4

Analyse du gène candidat de la schizophrénie

Gènes de la schizophrénie identifiés dans les méta-profils
ÂgeVers le basOpcml, Pldn, Nrg1, Rgs4, Bdnf, Dlg4, Gad67
En hautNtrk2, Ppp1r1b, Erbb3
pHVers le basNtrk2, Slc1a2
En hautRgs4, Gad67
PMIEn hautNrg1

Voies biologiques affectées

Pour dériver une interprétation biologique de haut niveau de nos méta-signatures, nous avons effectué une analyse d'enrichissement de l'ontologie génétique (GO) (Ashburner et al., 2000) en utilisant ErmineJ (Lee et al., 2005). Nous avons extrait les catégories ‘top’ GO pour chacun des profils et une comparaison entre elles a révélé divers ‘processus biologiques’ uniques à chaque méta-signature, en plus d'un certain nombre de catégories partagées. Dans la figure 4, nous avons affiché les dix premières catégories pour chaque méta-signature d'âge représentant le nombre de gènes représentés dans chaque catégorie GO et la valeur p associée (corrigée pour plusieurs tests). Les données ORA générées pour les huit méta-signatures se trouvent dans le tableau supplémentaire 5 à l'adresse http://chibi.ubc.ca/

Pour les méta-signatures d'âge, nous avons affiché les 10 principaux termes GO identifiés à l'aide d'une analyse de surreprésentation GO. L'axe des y principal affiche le nombre de gènes de méta-signature qui entrent dans la catégorie donnée « processus biologique » ; tandis que l'axe secondaire affiche la valeur p associée. Les termes GO étaient réduits au terme parent si le parent et l'enfant figuraient tous les deux dans les dix premiers.

Pour les gènes dont l'expression augmente avec la progression de l'âge, les principales catégories GO comprenaient ceux impliqués dans la croissance et la prolifération cellulaires, et l'interaction cellule-cellule, conformément aux études précédentes (Erraji-Benchekroun et al., 2005, Hong et al., 2008) . D'autres processus comprenaient la voie de signalisation des récepteurs de l'insuline, englobant un certain nombre de gènes impliqués dans la longévité et le vieillissement (Bartke, 2008). Les gènes régulés à la baisse par l'âge présentaient un enrichissement de l'activité synaptique et/ou des récepteurs avec des catégories GO telles que la reconnaissance des neurones, le transport et la sécrétion des neurotransmetteurs et la régulation des niveaux de neurotransmetteurs. Cette découverte est concordante avec les études de vieillissement existantes chez la souris et l'homme (Erraji-Benchekroun et al., 2005, Oh et al., 2009). De même, nous avons trouvé un enrichissement des gènes impliqués dans la sécrétion de neurotransmetteurs dans la méta-signature régulée à la hausse du pH, en plus des gènes impliqués dans le métabolisme et un éventail différent de voies (c'est-à-dire la signalisation des protéines G). Les méta-signatures féminines et masculines ont identifié un enrichissement similaire de termes incluant des processus pertinents tels que la détermination du sexe. Cependant, les gènes associés à ces termes partagés étaient assez différents entre les deux méta-signatures. L'analyse d'enrichissement fonctionnel de nos méta-signatures, bien que loin de fournir des preuves cellulaires solides, a tout de même fourni une indication utile des processus biologiques altérés par chaque facteur et une meilleure compréhension au niveau moléculaire.

Corrélation entre facteurs et méta-profils

Parce que nous avons analysé chaque facteur indépendamment, nous avons souhaité vérifier si les facteurs étaient corrélés entre eux dans les 415 échantillons (tableau 5). L'âge et le PMI affichaient la corrélation la plus élevée de 0,35, ce qui correspond à une corrélation positive rapportée dans (Galfalvy et al., 2003). L'âge et le pH du cerveau ont affiché une légère corrélation négative de 𢄠.2. Une étude plus approfondie de ces deux facteurs a révélé que la catégorisation des valeurs d'âge en groupes ‘young’ (< 50 ans) et ‘old’ (≥ 50 ans) résultait en un pH moyen plus faible. dans le groupe &# x02018old&# x02019 contre le &# x02018young&# x02019. C'était la tendance générale au sein de chaque ensemble de données (Figure supplémentaire 3).

Tableau 5

Classement des corrélations entre les facteurs

ÂgeSexepHPMI
Âge
Sexe𢄠.12 **
pH𢄠.2 * 0.17
PMI0.35 *** 𢄠.2 *** −.06

Les corrélations de rang de Spearman ont été calculées en utilisant les informations caractéristiques de l'échantillon pour les sujets individuels.

En raison de ces corrélations, nous nous attendions à ce que les gènes individuels dans certains méta-profils puissent se chevaucher avec d'autres méta-profils (tableau 6). En conséquence, les deux facteurs affichant le plus grand nombre de gènes chevauchants étaient ceux [en haut ou en bas] de l'âge et ceux [en bas ou en haut] avec le pH cérébral, respectivement. Cependant, ces effets étaient faibles et étaient encore plus faibles dans les signatures 𠆌ore’ (tableau supplémentaire 4). Nous avons également constaté qu'un certain nombre de gènes régulés à la hausse par le PMI ont également été identifiés parmi les profils de pH et d'âge du cerveau dans les deux sens, mais n'ont pu extraire aucun schéma défini.

Tableau 6

Évaluation du chevauchement des gènes entre les méta-signatures

Nombre de gènes de profilÂgepHPMISexe
En hautVers le basEn hautVers le basEn hautVers le basFemelleHomme
ÂgeEn haut404
Vers le bas11342
pHEn haut25018
Vers le bas2157510
PMIEn haut7523202
Vers le bas412000
SexeFemelle19410100
Homme140100000

En utilisant des gènes pour chaque méta-profil à q < 0,001, nous les avons comparés les uns aux autres pour évaluer le chevauchement et les relations potentielles entre les facteurs. L'observation de deux gènes dans la méta-signature d'âge qui changent dans les deux sens est une conséquence de plusieurs sondes de spécificité différente pour le même gène, qui montrent des directions de changement différentes.


Fond

Le rythme de génération de données dans les sciences de la vie augmente régulièrement. Les ensembles de données primaires augmentent en profondeur et en précision, couvrant de plus en plus d'aspects de la vie. En biologie moléculaire et en biomédecine, celles-ci incluent des mesures à grande échelle de l'acétylation de l'ADN/des histones, de l'activité transcriptionnelle, de l'expression des gènes et de l'abondance des protéines (par exemple [1]). La mesure des patrons épigénétiques (méthylation de l'ADN, acétylation de l'ADN/Histone) à grande échelle n'est devenue possible que récemment [1, 2]. La mesure de la transcription entre dans une nouvelle ère avec l'introduction du séquençage profond (ou de nouvelle génération, RNA-seq) [3, 4]. La protéomique devient possible à une profondeur sans précédent, couvrant de manière routinière des parties de plus en plus grandes du protéome [5]. Pour ces données primaires, des référentiels tels que la base de données Gene Expression Omnibus (GEO [6]) ou ArrayExpress [7] sont en constante expansion.

Souvent, les mesures sont différentielles : elles sont effectuées pour deux ou plusieurs conditions (telles que le knock-down ou le knock-out d'un gène [8]), pour deux moments ou plus (telles que les séries chronologiques suivant les conséquences d'une intervention expérimentale, [9]), ou pour deux espèces ou plus (comme la souris et l'homme, [10]). L'exploitation des mesures différentielles est une des clés pour faire face au flot de données, en se concentrant sur les différences les plus prononcées.

Les scientifiques de la vie doivent également gérer un déluge de données secondaires, sous forme d'articles, de revues et de bases de données organisées. Ceux-ci peuvent être intégrés par des systèmes automatisés tels que STRING [11], ou par des efforts manuels [12-14]. L'exploitation des données secondaires fournit une autre clé pour faire face au flot de données primaires, en les mettant en contexte et en se concentrant sur les confirmations et les contradictions les plus prononcées avec ce qui est déjà connu.

Dans cet article, nous proposons d'interpréter des données différentielles dans le contexte de la connaissance, donnant « l'essence » d'une expérience. Les données différentielles peuvent être fournies par deux puces à ADN, et les connaissances peuvent être fournies par un réseau décrivant l'interaction et la régulation gène/protéine. Dans ce cas, les données de suivi de l'expression génique au cours d'une expérience peuvent être utilisées pour identifier les mécanismes putatifs les plus prononcés. Ils sont identifiés comme ces liens connus entre gènes/protéines le long desquels des changements d'expression indiquent qu'il peut y avoir eu un changement de régulation, tel que le démarrage ou l'arrêt d'une interaction, une stimulation ou une inhibition. Expression met en évidence ces liens, et il permet à l'utilisateur de filtrer tous les liens sans changement ou négligeable. Plus le seuil du filtre sur la quantité de changement à afficher est élevé, moins les liens sont affichés, ce qui facilite l'examen de « l'essence » de l'expérience. Les condensations du réseau sont illustrées par des paires de figures (réseau d'origine - réseau condensé) dans la section sur les études de cas. Le réseau condensé contient de bons candidats pour interpréter l'expérience en termes mécanistes, donnant lieu à la conception de nouvelles expériences. Cependant, toutes les inférences sont des hypothèses dérivées de corrélations dans les données expérimentales dans le contexte de la a priori connaissances encodées dans le réseau, et il faut garder à l'esprit que les données corrélatives ne pas nécessairement impliquent une causalité mécaniste. De plus, la validité des hypothèses générées par notre méthode dépendra de la couverture et de la justesse du réseau, ainsi que de la précision des données expérimentales.

Travaux connexes

À partir des travaux pionniers d'Ideker et al. [15], il existe pléthore de méthodes qui combinent des données de réseau avec des données à haut débit (comme les microarrays), afin de mettre en évidence des voies ou des sous-réseaux, voir les excellentes revues récentes de Minguez & Dopazo [16], Wu et al. [17] et Yu & Li [18]. Notamment, peu de ces méthodes sont facilement disponibles sous forme de progiciels, de plug-ins ou de services Web accessibles au public (voir le tableau et dans [17]). De plus, il ne semble pas exister d'étalon-or pouvant être utilisé à des fins de validation (voir, par exemple, Tarca et al. [19] pour une discussion récente). Certaines méthodes manquent de validation à l'exception de l'exemple pour lequel elles ont été développées, tandis que d'autres sont étudiées pour un éventail d'exemples spécifiques. Dans ces cas, un fort enrichissement dans des catégories plausibles d'ontologie génétique ou la détection de voies ou d'annotations connues est souvent utilisé pour démontrer l'utilité, comme dans [19–25]. Nous avons trouvé deux articles comprenant une comparaison de différentes méthodes d'identification de sous-réseau. La première de Parkkinen et Kaski [26] introduit des variantes de la méthode Interaction Component Model (ICM), en les comparant à la méthode ICM originale, à une méthode basée sur les champs aléatoires modulaires cachés (HMoF) [27] et à Matisse [28 ], en utilisant l'identification des classes d'ontologie génétique et la couverture des complexes protéiques pour deux ensembles de données sélectionnés (réponse au choc osmotique et données sur les dommages à l'ADN) pour juger une méthode par rapport à l'autre. Une évaluation de ClustEx [29], jActiveModules [15], GXNA [21] et une approche simple basée sur le changement de pli peuvent être trouvées dans [29], en identifiant des ensembles de gènes, des voies et des cibles de microarray connus de la littérature et de la Gene Ontology pour comparaison.

En général, il est extrêmement difficile de valider la détection de (sous-)réseaux ou de (sous-)voies : ce sont des entités complexes, et la validation expérimentale finale est impossible en raison de cette complexité : les expérimentateurs se limitent généralement à n'étudier que quelques composants dans l'isolement à un moment donné. Néanmoins, nous comparerons les résultats de notre méthode avec les résultats obtenus par jModules Actifs, dans une section distincte à la suite des études de cas. En revanche, en mettant simplement en évidence des liens uniques dans les réseaux, nous abordons une tâche plus primitive, mais dans ce cas les résultats peuvent être validés directement par l'expérience, ou en identifiant des déclarations corroborantes dans la littérature. En particulier, comme le montrent nos études de cas, les liens uniques que nous mettons en évidence donnent lieu à des prédictions sur des gènes uniques et sur des mécanismes à une seule étape qui peuvent être étudiés isolément. Par conséquent, nous voudrions souligner l'utilité directe de notre focalisation sur les liens et les gènes uniques, en complément de la vision centrée sur le (sous-)réseau qui est généralement utilisée au meilleur de notre connaissance, la focalisation « lien et gène unique » n'est pas utilisée par d'autres méthodes combinant des données réseau et à haut débit (« omiques »). En fait, nous proposons une «combinaison gagnante» de biologie «en réseau»/«omique» et «classique», en utilisant des réseaux et des données à haut débit pour mettre en évidence des gènes et des liens uniques qui peuvent ensuite être validés directement par la biologie moléculaire classique, comme le être démontré dans nos études de cas.

Comme travaux futurs, notre formule de mise en évidence des liens peut cependant être intégré dans les méthodes actuelles de détection des voies/sous-réseau, en les améliorant peut-être considérablement. En particulier, aucune de ces méthodes ne traite les inhibitions et les stimulations de manière distincte, comme nous le faisons. En particulier, nous envisageons que la formule de score de bord de Guo et al. [20], qui repose sur la mesure de la covariance, peut être remplacé par notre formule (voir ci-dessous), mettant l'accent sur un aspect différent de l'expression différentielle des gènes : alors que Guo et al. identifier coordonné changements en utilisant leur formule, l'intégration de notre formule dans leur cadre identifierait les sous-réseaux avec des changements qui sont cohérent avec un réseau d'entrée d'interactions, de stimulations et d'inhibitions. Dans tous les cas, nous souhaitons souligner que pour l'identification des changements coordonnés, les coefficients de corrélation sont les plus appropriés. Notre approche, cependant, identifie un message biologique différent, à savoir les démarrages/arrêts d'interactions, de stimulations et d'inhibitions, en utilisant un réseau d'entrée informatif sur les relations biologiques telles que les stimulations et les inhibitions.


RÉSULTATS

Présentation de MACAO

Les deux termes d'effets aléatoires |$oldsymbol$| et |$oldsymbol$| sur-dispersion du modèle, la variance supplémentaire non expliquée par un modèle de Poisson. Cependant, les deux termes |$oldsymbol$| et |$oldsymbol$| modéliser deux aspects différents de la surdispersion. Plus précisément, |$oldsymbol$| modélise la fraction de la variance supplémentaire qui est expliquée par la non-indépendance de l'échantillon tandis que |$oldsymbol$| modélise la fraction de la variance supplémentaire qui est indépendante entre les échantillons. Par exemple, considérons un cas simple dans lequel tous les échantillons ont la même profondeur de séquençage (c'est-à-dire |$ = N$|⁠ ) et il n'y a qu'un seul terme d'interception |$mu$| inclus comme covariable. Dans ce cas, le terme à effets aléatoires |$oldsymbol$| modélise la surdispersion indépendante : sans |$oldsymbol$|⁠ , |$V ( y ) = E( y )( <1 + E( y )( <>> - 1> )> )$| est toujours plus grand que la moyenne |$E( y ) = N/2>>$|⁠ , la différence entre les deux augmentant avec |$$|⁠ . De la même manière, le terme à effets aléatoires |$oldsymbol$| modélise la surdispersion non indépendante en tenant compte de la matrice de covariance de l'échantillon |$oldsymbol$|⁠ . En modélisant les deux aspects de la surdispersion, notre PMM généralise efficacement le modèle binomial négatif couramment utilisé - qui ne modélise que la variance supplémentaire indépendante - pour tenir compte de la non-indépendance de l'échantillon. De plus, notre PMM étend naturellement le LMM couramment utilisé (9, 64, 83, 84) à la modélisation des données de comptage.

Notre objectif ici est de tester l'hypothèse nulle selon laquelle les niveaux d'expression des gènes ne sont pas associés à la variable prédictive d'intérêt, ou |$:eta = < m< >>0$|⁠ . Le test de cette hypothèse nécessite l'estimation des paramètres dans le PMM (comme cela a déjà été fait dans d'autres paramètres (85, 86), y compris pour la modélisation de la distribution inégale de lecture de RNAseq dans les détails des transcriptions (13) dans les données supplémentaires). Le PMM appartient à la famille LMM généralisée, où l'estimation des paramètres est notoirement difficile en raison des effets aléatoires et de l'insoluble qui en résulte. m-intégrale dimensionnelle dans la vraisemblance. Les méthodes d'estimation standard reposent sur l'intégration numérique ( 87) ou l'approximation de Laplace ( 88, 89), mais aucune de ces stratégies ne s'adapte bien à la dimension croissante de l'intégrale, qui dans notre cas est égale à la taille de l'échantillon. En conséquence, les approches standard produisent souvent des estimations biaisées et des intervalles de confiance trop étroits (c'est-à-dire anti-conservateurs) (90-96). Pour surmonter la haute dimensionnalité de l'intégrale, nous développons plutôt un nouvel algorithme de Markov Chain Monte Carlo (MCMC), qui, avec suffisamment d'itérations, peut atteindre une précision d'inférence élevée (97, 98). Nous utilisons MCMC pour tirer des échantillons a posteriori mais nous nous appuyons sur la normalité asymptotique des distributions de vraisemblance et a posteriori ( 99) pour obtenir l'estimation approximative du maximum de vraisemblance |$_j>$| et son erreur standard se( ⁠|$widehat <<>_j>>$|⁠ ). Avec |$_j>$| et se( ⁠|$widehat <<>_j>>$|⁠ ), nous pouvons construire des statistiques de test de Wald approximatives et P-valeurs pour les tests d'hypothèse (matériel supplémentaire). Bien que nous utilisions MCMC, notre procédure est de nature fréquentiste.

Au niveau technique, notre algorithme MCMC est également nouveau, tirant parti d'une représentation variable auxiliaire de la vraisemblance de Poisson ( 61–63) et des innovations récentes en algèbre linéaire pour l'ajustement des LMM (9, 58, 64). Notre algorithme MCMC introduit deux variables latentes continues pour chaque individu pour remplacer l'observation du décompte, étendant efficacement notre approche précédente consistant à utiliser une variable latente pour la distribution binomiale plus simple ( 59). Par rapport à un MCMC standard, notre nouvel algorithme MCMC réduit la complexité de calcul de chaque itération MCMC de cubique à quadratique par rapport à la taille de l'échantillon. Par conséquent, notre méthode est de plusieurs ordres de grandeur plus rapide que le logiciel bayésien populaire MCMCglmm (100) et peut être utilisée pour analyser des centaines d'échantillons et des dizaines de milliers de gènes avec un seul PC de bureau (Figure supplémentaire S1). Bien que notre procédure soit de nature stochastique, nous trouvons que les erreurs MCMC sont souvent suffisamment petites pour assurer une stabilité P-valeurs à travers des exécutions MCMC indépendantes (Figure supplémentaire S2). Nous résumons les principales caractéristiques de notre méthode ainsi que d'autres approches couramment utilisées dans le tableau 1.

Simulations : contrôle de la non-indépendance de l'échantillon

Nous avons effectué une série de simulations pour comparer les performances du PMM implémenté à MACAU avec quatre autres méthodes couramment utilisées : (i) un modèle linéaire (ii) le LMM implémenté dans GEMMA (9, 58) (iii) un modèle de Poisson et ( iv) un modèle binomial négatif. Nous avons utilisé des données transformées en quantiles pour le modèle linéaire et GEMMA (voir la section « Matériaux et méthodes » pour les détails de la normalisation et une comparaison entre diverses transformations Figure supplémentaire S3) et utilisé des données de comptage brutes pour les trois autres méthodes. Pour rendre nos simulations réalistes, nous utilisons des paramètres déduits d'un ensemble de données RNAseq publié sur une population de babouins sauvages (7, 71) pour effectuer des simulations (section « Matériaux et méthodes »). exemple de non-indépendance décrit ci-dessus.

Notre première série de simulations a été réalisée pour évaluer l'efficacité de MACAU et des quatre autres méthodes pour contrôler la non-indépendance de l'échantillon. Pour ce faire, nous avons simulé les niveaux d'expression de 10 000 gènes chez 63 individus (la taille de l'échantillon de l'ensemble de données sur le babouin). Les niveaux d'expression génique simulés sont influencés à la fois par des effets environnementaux indépendants et des effets génétiques corrélés, où les effets génétiques sont simulés sur la base de la matrice de parenté du babouin (estimée à partir des données microsatellites (7)) avec soit zéro ( |$ = 0.0$|⁠ ), modéré ( |$ = 0.3$|⁠ ), ou élevé ( ⁠|$ = 0.6$|⁠ ) valeurs d'héritabilité. Nous avons également simulé une variable prédictive continue x qui est elle-même modérément héritable ( |$h_x^2 = 0.4$|⁠ ). Parce que nous étions intéressés par le comportement du nul dans cet ensemble de simulations, les niveaux d'expression des gènes n'étaient pas affectés par la variable prédictive (c'est-à-dire qu'aucun gène n'était vraiment DE).

La figure 1, les figures supplémentaires S4 et 5 montrent des graphiques quantile-quantile pour les analyses utilisant MACAU et les quatre autres méthodes par rapport à l'espérance nulle (uniforme), pour |$ = 0.6$|⁠ , |$ = 0.3$| et |$ = 0.0$| respectivement. Lorsque les gènes sont héréditaires et que la variable prédictive est également corrélée à la parenté individuelle, alors le résultat P- les valeurs de l'analyse DE ne devraient être uniformes que pour une méthode qui contrôle correctement la non-indépendance de l'échantillon. Si une méthode ne parvient pas à contrôler la non-indépendance de l'échantillon, le P- les valeurs seraient gonflées, ce qui entraînerait des faux positifs.

QQ-plots comparant attendu et observé P-les distributions de valeurs générées par différentes méthodes pour les simulations nulles en présence de non-indépendance de l'échantillon. Dans chaque cas, 10 000 gènes non-DE ont été simulés avec m = 63, CV = 0,3, 2 = 0.25, h 2 = 0,6, et hX 2 = 0,4. Les méthodes de comparaison incluent MACAO (UNE), binôme négatif (B), Poisson (C), GEMMA () et Linéaire (E). MACAU et GEMMA contrôlent correctement l'erreur de type I en présence de non-indépendance de l'échantillon. ??gc est le facteur de contrôle génomique.

QQ-plots comparant attendu et observé P-les distributions de valeurs générées par différentes méthodes pour les simulations nulles en présence de non-indépendance de l'échantillon. Dans chaque cas, 10 000 gènes non-DE ont été simulés avec m = 63, CV = 0,3, 2 = 0.25, h 2 = 0,6, et hX 2 = 0,4. Les méthodes de comparaison incluent MACAO (UNE), binôme négatif (B), Poisson (C), GEMMA () et Linéaire (E). MACAU et GEMMA contrôlent correctement l'erreur de type I en présence de non-indépendance de l'échantillon. ??gc est le facteur de contrôle génomique.

Nos résultats montrent que, parce que MACAU contrôle la non-indépendance de l'échantillon, le P-les valeurs de MACAU suivent de près la distribution uniforme attendue (et sont légèrement conservatrices), que l'expression des gènes soit modérément ou hautement héritable. Les facteurs de contrôle génomique de MACAU sont proches de 1 (figure 1 et figure supplémentaire S4). Même si nous utilisons un style relativement détendu q-valeur seuil de 0,2 pour identifier les gènes DE, nous n'identifions à tort aucun gène comme DE avec MACAU. En revanche, le P-les valeurs du binôme négatif sont gonflées et biaisées vers des valeurs faibles (significatives), en particulier pour les niveaux d'expression génique avec une héritabilité élevée. Avec un binôme négatif, 27 gènes DE (quand h 2 = 0,3) ou 21 gènes DE (quand h 2 = 0,6) sont détectés par erreur au q-valeur seuil de 0,2. L'inflation de P-valeurs est encore plus aiguë dans Poisson, probablement parce que le modèle de Poisson ne tient compte ni de la parenté individuelle ni de la surdispersion. Pour les modèles non basés sur le nombre, le P-les valeurs d'un modèle linéaire sont légèrement biaisées vers des valeurs significatives, avec trois gènes DE (quand h 2 = 0,3) et un gène DE (quand h 2 = 0,6) détecté par erreur à q < 0.2. En revanche, comme le LMM dans GEMMA tient également compte de la parenté individuelle, il contrôle bien la non-indépendance de l'échantillon. Enfin, lorsque les gènes ne sont pas héritables, toutes les méthodes, à l'exception de Poisson, contrôlent correctement l'erreur de type I ( figure supplémentaire S5 ).

Deux facteurs importants influencent la gravité de la non-indépendance de l'échantillon dans les données RNAseq (Figure 2). Premièrement, l'inflation des P-les valeurs dans les modèles binomial négatif, de Poisson et linéaire deviennent plus aiguës avec l'augmentation de la taille de l'échantillon. En particulier, lorsque |$h_x^2 = 0.4$|⁠ , avec une taille d'échantillon de |$n = 1,000$|⁠ , |$>$| du binomial négatif, les modèles de Poisson et linéaires atteignent respectivement 1,71, 82,28 et 1,41. En revanche, même lorsque |$n = 1,000$|⁠ , |$>$| de MACAO et de GEMMA restent proches de 1 (0,97 et 1,01, respectivement). Deuxièmement, l'inflation des PLes valeurs des trois modèles deviennent également plus aiguës lorsque la variable prédictive est plus corrélée à la structure de la population. Ainsi, pour une variable prédictive hautement héritable ( ⁠|$h_x^2 = 0.8$|⁠ ), |$>$| (lorsque |$n = 1000$|⁠ ) des modèles binomial négatif, de Poisson et linéaire augmente à 2,13, 101,43 et 1,81, respectivement, alors que |$>$| de MACAO et GEMMA reste proche de 1 (1,02 et 1,05).

Comparaison du facteur de contrôle génomique λgc à partir de différentes méthodes pour les simulations nulles en présence de non-indépendance de l'échantillon. 10 000 gènes nuls ont été simulés avec CV = 0,3, 2 = 0.25, h 2 = 0,6, et (UNE) hX 2 = 0 (B) hX 2 = 0.4 (C) hX 2 = 0.8. ??gc (axe des y) change avec la taille de l'échantillon m (axe des x). Les méthodes de comparaison étaient MACAU (rouge), binôme négatif (violet), GEMMA (bleu) et linéaire (cyan). MACAU et GEMMA fournissent des statistiques de test calibrées en présence de non-indépendance de l'échantillon dans une gamme de paramètres. ??gc de Poisson dépasse 10 dans tous les paramètres et n'est donc pas représenté.

Comparaison du facteur de contrôle génomique λgc à partir de différentes méthodes pour les simulations nulles en présence de non-indépendance de l'échantillon. 10 000 gènes nuls ont été simulés avec CV = 0,3, 2 = 0.25, h 2 = 0,6, et (UNE) hX 2 = 0 (B) hX 2 = 0.4 (C) hX 2 = 0.8. ??gc (axe des y) change avec la taille de l'échantillon m (axe des x). Les méthodes de comparaison étaient MACAU (rouge), binôme négatif (violet), GEMMA (bleu) et linéaire (cyan). MACAU et GEMMA fournissent tous deux des statistiques de test calibrées en présence de non-indépendance de l'échantillon dans une gamme de paramètres. ??gc de Poisson dépasse 10 dans tous les paramètres et n'est donc pas représenté.

Nous avons également comparé MACAU avec edgeR (25) et DESeq2 (15), deux méthodes couramment utilisées pour l'analyse DE (11, 101). Comme edgeR et DESeq2 ont été conçus pour des valeurs de prédicteur discret, nous avons discrétisé le prédicteur continu |$x$| en 0/1 en fonction de la valeur prédictive médiane des individus. Nous avons ensuite appliqué toutes les méthodes aux mêmes valeurs de prédicteur binarisées à des fins de comparaison. Les résultats sont présentés dans la figure supplémentaire S6 . Pour les cinq méthodes comparées ci-dessus, les résultats sur les valeurs binarisées sont comparables à ceux des variables continues (c'est-à-dire la figure supplémentaire S6 par rapport à la figure 1). EdgeR et DESeq2 produisent tous deux un anticonservateur P-valeurs et fonctionnent de manière similaire au modèle binomial négatif en termes de contrôle d'erreur de type I.

Enfin, nous avons exploré l'utilisation de PC de la matrice d'expression génique ou de la matrice de génotype pour contrôler la non-indépendance de l'échantillon. Les PC de génotype ont été utilisés comme covariables pour contrôler la stratification de la population dans les études d'association ( 102). Cependant, des études comparatives récentes ont montré que l'utilisation de PC est moins efficace que l'utilisation de LMM (83, 103). Conformément aux performances plus faibles des CP dans les études d'association (83, 103), l'utilisation des meilleures CP de la matrice d'expression génique ou de la matrice de génotype n'améliore pas le contrôle des erreurs de type I pour les approches binomiales négatives, Poisson, linéaire, edgeR ou DESeq2 ( Figures supplémentaires S7 et 8 ).

Simulations : pouvoir d'identifier les gènes DE

Notre deuxième série de simulations a été conçue pour comparer la puissance de différentes méthodes d'identification des gènes DE, encore une fois sur la base de paramètres déduits de données réelles. Cette fois, nous avons simulé un total de 10 000 gènes, parmi lesquels 1000 gènes étaient vraiment DE et 9000 étaient non-DE. Pour les gènes DE, les tailles d'effet simulées correspondaient à une PVE fixe dans les niveaux d'expression génique qui variaient de 15 à 35%. Pour chaque ensemble de paramètres, nous avons effectué 10 simulations répétées et mesuré les performances du modèle sur la base de l'AUC (comme dans ( 35, 68, 104)). Nous avons également examiné plusieurs facteurs clés qui pourraient influencer la performance relative des méthodes alternatives : (i) l'héritabilité de l'expression génique ( ⁠|$$|⁠ ) (ii) corrélation entre la variable prédictive |$x$| et la parenté génétique (mesurée par l'héritabilité de |$x$|⁠ , ou |$h_x^2$|⁠ ) (iii) la variation des CRT entre les échantillons (mesurée par le CV) (iv) le paramètre de surdispersion ( ⁠|$$|⁠ ) (v) la taille de l'effet (PVE) et (vi) la taille de l'échantillon (n). Pour ce faire, nous avons d'abord effectué des simulations en utilisant un ensemble de valeurs par défaut ( ⁠|$ =$| 0.3, |$h_x^2$| = 0, CV = 0,3, |$< m< >> =$| 0,25, PVE = 0,25 et m = 63) puis les a modifiés un par un pour examiner l'influence de chaque facteur sur la performance relative de chaque méthode.

Nos résultats montrent que MACAU fonctionne aussi bien ou mieux que les autres méthodes dans presque tous les contextes (Figure 3 et Figure supplémentaire S9-14), probablement parce que les deux modèles comptent les données directement et contrôlent la non-indépendance de l'échantillon. En revanche, l'approche de Poisson a toujours obtenu les pires résultats dans tous les scénarios de simulation, probablement parce qu'elle ne tient pas compte des sources de surdispersion (figure 3 et figures supplémentaires S9-14).

MACAU présente une puissance accrue pour détecter les gènes DE vrais positifs dans une gamme de paramètres de simulation. L'aire sous la courbe (AUC) est indiquée comme mesure de performance pour MACAU (rouge), binôme négatif (violet), Poisson (vert), GEMMA (bleu) et linéaire (cyan). Chaque paramètre de simulation se compose de 10 réplicats de simulation, et chaque réplicat comprend 10 000 gènes simulés, avec 1 000 DE et 9 000 non-DE. Nous avons utilisé m = 63, hX 2 = 0,4, PVE = 0,25 et 2 = 0,25. Dans (UNE) nous avons augmenté h 2 en maintenant CV = 0,3 et en (B) nous avons augmenté le CV tout en maintenant h 2 = 0,3. Les boîtes à moustaches de l'ASC à travers les réplicats pour différentes méthodes montrent que (A) l'héritabilité (h 2 ) influe sur la performance relative des méthodes qui tiennent compte de la non-indépendance de l'échantillon (MACAU et GEMMA) par rapport aux méthodes qui ne font pas (binomial négatif, Poisson, linéaire) (B) la variation du nombre total de lectures entre les individus, mesurée par le coefficient de variation (CV), influence la performance relative de GEMMA et du binôme négatif. Les encarts dans les deux figures montrent le rang des différentes méthodes, où la rangée du haut représente le rang le plus élevé.

MACAU présente une puissance accrue pour détecter les véritables gènes DE positifs dans une gamme de paramètres de simulation. L'aire sous la courbe (AUC) est indiquée comme mesure de performance pour MACAU (rouge), binôme négatif (violet), Poisson (vert), GEMMA (bleu) et linéaire (cyan). Chaque paramètre de simulation se compose de 10 réplicats de simulation, et chaque réplicat comprend 10 000 gènes simulés, avec 1 000 DE et 9 000 non-DE. Nous avons utilisé m = 63, hX 2 = 0,4, PVE = 0,25 et 2 = 0,25. Dans (UNE) nous avons augmenté h 2 en maintenant CV = 0,3 et en (B) nous avons augmenté le CV tout en maintenant h 2 = 0,3. Les boîtes à moustaches de l'ASC à travers les réplicats pour différentes méthodes montrent que (A) l'héritabilité (h 2 ) influe sur la performance relative des méthodes qui tiennent compte de la non-indépendance de l'échantillon (MACAU et GEMMA) par rapport aux méthodes qui ne font pas (binomial négatif, Poisson, linéaire) (B) la variation du nombre total de lectures entre les individus, mesurée par le coefficient de variation (CV), influence la performance relative de GEMMA et du binôme négatif. Les encarts dans les deux figures montrent le rang des différentes méthodes, où la rangée du haut représente le rang le plus élevé.

Parmi les facteurs qui influencent le rang relatif des diverses méthodes, le facteur le plus important était l'héritabilité |$< m< >>( <> )$| (Figure 3A). Alors que toutes les méthodes fonctionnent moins bien avec l'augmentation de l'héritabilité de l'expression génique, l'héritabilité affecte de manière disproportionnée les performances des modèles qui ne tiennent pas compte de la parenté (c. = 0$|⁠ ), ces approches ont tendance à être légèrement plus performantes. Par conséquent, pour les gènes non héritables, les modèles linéaires fonctionnent légèrement mieux que GEMMA et les modèles binomiaux négatifs fonctionnent de manière similaire ou légèrement mieux que MACAU. Cette observation est probablement due au fait que les modèles binomiaux linéaires et négatifs nécessitent moins de paramètres et ont donc un plus grand nombre de degrés de liberté. Cependant, même dans ce cadre, la différence entre MACAU et binomial négatif est faible, ce qui suggère que MACAU est robuste pour modéliser les erreurs de spécification et fonctionne raisonnablement bien même pour les gènes non héritables.En revanche, lorsque l'héritabilité est moyenne ( |$ = 0.3$|⁠ ) ou élevé ( |$ = 0.6$|⁠ ), les méthodes qui tiennent compte de la non-indépendance de l'échantillon sont beaucoup plus puissantes que les méthodes qui ne le font pas. Étant donné que presque tous les gènes sont influencés par les cis-eQTL (46, 47) et sont donc probablement héréditaires dans une certaine mesure, la robustesse de MACAU pour les gènes non héréditaires et son gain de performance élevé pour les gènes héréditaires le rendent attrayant.

Le deuxième facteur le plus important dans la performance relative du modèle était la variation des CRT entre les individus (CV Figure 3B). Alors que toutes les méthodes fonctionnent moins bien avec l'augmentation du CV, le CV affecte particulièrement les performances de GEMMA. Plus précisément, lorsque CV est petit (0,3 en tant que données de babouin), GEMMA fonctionne bien et est la deuxième meilleure méthode derrière MACAU. Cependant, lorsque le CV est modéré (0,6) ou élevé (0,9), les performances de GEMMA décroissent rapidement : elle ne devient que la quatrième meilleure méthode lorsque CV = 0,9. GEMMA fonctionne mal dans des paramètres de CV élevés, probablement parce que le LMM ne tient pas compte de la dépendance moyenne-variance observée dans les données de comptage, ce qui est en accord avec les résultats précédents (59, 105).

Les quatre autres facteurs que nous avons explorés ont eu de faibles impacts sur la performance relative des méthodes alternatives, bien qu'ils aient affecté leur performance absolue. Par exemple, comme on pouvait s'y attendre, la puissance augmente avec de grandes tailles d'effet (PVE) (figure supplémentaire S9 ) ou de grandes tailles d'échantillon (figure supplémentaire S10 ), et diminue avec une grande surdispersion |$$| ( Figure supplémentaire S11 ) ou grand |$h_x^2$| (Figure supplémentaire S12).

Enfin, nous avons inclus des comparaisons avec edgeR (25) et DESeq2 (15). Dans le réglage de simulation des paramètres de base (m = 63, CV = 0,3, |$h_x^2$| = 0, PVE = 0,25, |$ =$| 0.3 et |$< m< >> =$| 0.25), nous avons à nouveau discrétisé le prédicteur continu |$x$| en une variable binaire 0/1 basée sur la valeur prédictive médiane des individus. Les résultats de toutes les méthodes sont présentés dans la figure supplémentaire S13A. Pour les cinq méthodes également testées sur une variable prédictive continue, la puissance sur les valeurs binarisées est très réduite par rapport à la puissance lorsque la variable prédictive est modélisée sans binarisation (par exemple, figure supplémentaire S13A par rapport à la figure 3). De plus, ni edgeR ni DESeq2 ne fonctionnent bien, ce qui est cohérent avec le passage récent de ces méthodes vers des modèles linéaires dans l'analyse de l'expression différentielle (3, 7, 45-47, 106). Ce résultat ne dépend pas d'un échantillon de grande taille. Dans les paramètres de petite taille d'échantillon (m = 6, m = 10 et m = 14, avec des échantillons équilibrés entre les deux classes, 0 ou 1), MACAU surpasse à nouveau les autres méthodes, bien que la différence de puissance soit beaucoup plus petite (m = 10 et m = 14 figures supplémentaires S13C et 13D ) et parfois négligeable (m = 6, figure supplémentaire S13B).

En résumé, la puissance de MACAU et d'autres méthodes, ainsi que la différence de puissance entre les méthodes, est influencée de manière continue par de multiples facteurs. Des tailles d'échantillons plus grandes, des tailles d'effet plus grandes, une variation de profondeur de lecture plus faible, une héritabilité de l'expression génique plus faible, une héritabilité des variables prédictives plus faible et une surdispersion plus faible augmentent tous la puissance. Cependant, la puissance de MACAU est moins diminuée par l'héritabilité élevée de l'expression génique et la variabilité élevée de la profondeur de lecture que les méthodes de modèle non mixtes, tout en conservant l'avantage de modéliser directement les données de comptage. Dans des contextes d'analyse de données difficiles (par exemple, lorsque la taille de l'échantillon est faible et la taille de l'effet est faible : Figure supplémentaire S13B pour m = 6), aucune méthode ne se démarque et l'utilisation de MACAU entraîne des gains de puissance nuls ou négligeables par rapport aux autres méthodes. Lorsque la taille de l'échantillon est faible (m = 6) et les tailles d'effet sont grandes, cependant, MACAO surpasse systématiquement les autres méthodes (m = 6, figure supplémentaire S14).

Applications de données réelles

Pour mieux comprendre au-delà de la simulation, nous avons appliqué MACAU et les six autres méthodes à trois ensembles de données RNAseq récemment publiés.

Le premier ensemble de données que nous avons considéré est l'étude RNAseq du babouin (7) utilisée pour paramétrer les simulations ci-dessus. Des mesures d'expression sur 12 018 gènes exprimés dans le sang ont été collectées par l'ABRP (71) pour 63 babouins adultes (26 femelles et 37 mâles), dont certains étaient apparentés. Ici, nous avons appliqué MACAU et les six autres méthodes pour identifier les gènes avec des modèles d'expression biaisés par le sexe. Les gènes associés au sexe sont connus pour être enrichis sur les chromosomes sexuels ( 107, 108), et nous utilisons cet enrichissement comme l'un des critères pour comparer les performances des méthodes, comme dans (18). Parce que la même valeur nominale P-valeur de différentes méthodes peut correspondre à différentes erreurs de type I, nous avons comparé des méthodes basées sur des FDR empiriques. En particulier, nous avons permuté les données pour construire une valeur nulle empirique, estimé le FDR à un moment donné P-valeur seuil, et compté le nombre de découvertes à un seuil FDR donné (voir la section « Matériaux et méthodes » pour les détails de la permutation et une comparaison entre deux procédures de permutation différentes, figure supplémentaire S15).

En accord avec nos simulations, MACAU était la méthode la plus puissante de celles que nous avons considérées. Plus précisément, à un FDR empirique de 5%, MACAU a identifié 105 gènes avec des modèles d'expression biaisés par le sexe, 40% de plus que celui identifié par le modèle linéaire, la deuxième meilleure méthode à ce seuil FDR (Figure 4A). À un FDR plus détendu de 10 %, MACAU a identifié 234 gènes associés au sexe, 47 % de plus que celui identifié par le modèle binomial négatif, la deuxième meilleure méthode à ce seuil de FDR (figure 4A). De plus, comme prévu, les gènes associés au sexe détectés par MACAU sont enrichis sur le chromosome X (le chromosome Y n'est pas assemblé chez les babouins et est donc ignoré), et cet enrichissement est plus fort pour les gènes identifiés par MACAU que par les autres méthodes (Figure 4B). Parmi les approches restantes, le modèle linéaire binomial négatif et GEMMA ont tous eu des performances similaires et sont classés juste après MACAU. Le modèle de Poisson effectue le pire, et edgeR et DESeq2 se situent entre le modèle de Poisson et les autres méthodes (Figure 4A et B).

MACAU identifie plus de gènes exprimés de manière différentielle que d'autres méthodes chez le babouin (UNE-B) et FUSION (C-F) ensembles de données. Les méthodes de comparaison incluent MACAU (rouge), binôme négatif (violet), Poisson (vert), edgeR (magenta), DESeq2 (brun rosé), GEMMA (bleu) et linéaire (cyan). (A) montre le nombre de gènes associés au sexe identifiés par différentes méthodes à une gamme de taux de fausses découvertes empiriques (FDR). (B) montre le nombre de gènes qui se trouvent sur le chromosome X parmi les gènes qui ont la plus forte association sexuelle pour chaque méthode (notez que le chromosome Y n'est pas assemblé chez les babouins et est donc ignoré). Par exemple, dans les 400 premiers gènes identifiés par MACAU, 41 d'entre eux sont également sur le chromosome X. (C) montre le nombre de gènes associés au DT2 identifiés par différentes méthodes à une gamme de taux de fausses découvertes empiriques (FDR). (D) montre le nombre de gènes qui figurent dans la liste des 1 000 gènes les plus significativement associés à la GL parmi les gènes qui ont la plus forte association pour le DT2 pour chaque méthode. Par exemple, dans les 1 000 premiers gènes avec la plus forte association DT2 identifiés par MACAU, 428 d'entre eux figurent également dans la liste des 1 000 premiers gènes avec la plus forte association GL identifiée par la même méthode. (E) montre le nombre de gènes associés à la GL identifiés par différentes méthodes dans une gamme de FDR. (F) montre le nombre de gènes qui sont dans la liste des 1 000 gènes les plus significativement associés au DT2 parmi les gènes qui ont l'association la plus forte pour GL pour chaque méthode. DT2 : diabète de type II GL : glycémie à jeun.

MACAU identifie plus de gènes exprimés de manière différentielle que d'autres méthodes chez le babouin (UNE-B) et FUSION (C-F) ensembles de données. Les méthodes de comparaison incluent MACAU (rouge), binôme négatif (violet), Poisson (vert), edgeR (magenta), DESeq2 (brun rosé), GEMMA (bleu) et linéaire (cyan). (A) montre le nombre de gènes associés au sexe identifiés par différentes méthodes à une gamme de taux de fausses découvertes empiriques (FDR). (B) montre le nombre de gènes qui se trouvent sur le chromosome X parmi les gènes qui ont la plus forte association sexuelle pour chaque méthode (notez que le chromosome Y n'est pas assemblé chez les babouins et est donc ignoré). Par exemple, dans les 400 premiers gènes identifiés par MACAU, 41 d'entre eux sont également sur le chromosome X. (C) montre le nombre de gènes associés au DT2 identifiés par différentes méthodes à une gamme de taux de fausses découvertes empiriques (FDR). (D) montre le nombre de gènes qui figurent dans la liste des 1 000 gènes les plus significativement associés à la GL parmi les gènes qui ont la plus forte association pour le DT2 pour chaque méthode. Par exemple, dans les 1 000 premiers gènes avec la plus forte association DT2 identifiés par MACAU, 428 d'entre eux figurent également dans la liste des 1 000 premiers gènes avec la plus forte association GL identifiée par la même méthode. (E) montre le nombre de gènes associés à la GL identifiés par différentes méthodes dans une gamme de FDR. (F) montre le nombre de gènes qui sont dans la liste des 1 000 gènes les plus significativement associés au DT2 parmi les gènes qui ont la plus forte association pour GL pour chaque méthode. DT2 : diabète de type II GL : glycémie à jeun.

Le deuxième jeu de données que nous avons considéré est une étude RNAseq sur le DT2 collectée dans le cadre de l'étude FUSION (60). Ici, les données ont été recueillies à partir d'échantillons de muscles squelettiques de 267 individus avec des mesures d'expression sur 21 753 gènes. Les individus sont originaires de trois municipalités (Helsinki, Savitaipale et Kuopio) en Finlande. Les individus au sein de chaque municipalité sont plus étroitement liés que les individus entre les municipalités (par exemple, les PC de génotype supérieur correspondent généralement aux trois municipalités Figure supplémentaire S16). Deux phénotypes apparentés étaient à notre disposition : 162 individus avec le statut DT2 ou NGT (c'est-à-dire cas/témoin) sur la base de l'OGTT et 267 individus avec le trait quantitatif GL à jeun, un trait biologiquement pertinent du DT2.

Nous avons effectué des analyses pour identifier les gènes associés au statut DT2 ainsi que les gènes associés à la GL. Pour prendre en charge edgeR et DESeq2, nous avons également discrétisé les valeurs GL continues en catégories binaires 0/1 en fonction de la valeur GL médiane des individus. Nous appelons les valeurs résultantes GL01. Par conséquent, nous avons effectué deux séries d'analyses pour GL : une sur les valeurs continues de GL et l'autre sur les valeurs discrétisées de GL01. Conformément aux simulations et à l'analyse des données sur les babouins, MACAU a identifié plus de gènes associés au DT2 et de gènes associés à la GL que d'autres méthodes sur une gamme de valeurs FDR empiriques. Pour l'analyse T2D, MACAU a identifié 23 gènes associés au T2D avec un FDR de 5%, tandis que GEMMA et le modèle linéaire, les deuxièmes meilleures méthodes à ce seuil FDR, n'ont identifié qu'un seul gène associé au T2D (Figure 4C). De même, à un FDR de 10 %, MACAU a identifié 123 gènes associés au T2D, 51 % de plus que celui identifié par le modèle linéaire, la deuxième meilleure méthode à ce seuil de FDR (figure 4C). Pour l'analyse GL, basée sur un FDR de 5%, MACAU a détecté 12 gènes DE, tandis que les autres méthodes n'ont identifié aucun gène DE à ce seuil de FDR. À un FDR de 10 %, MACAU a identifié 100 gènes associés au GL, tandis que les deuxièmes meilleures méthodes, le modèle linéaire et GEMMA, ont identifié 12 gènes DE (figure 4E). Pour le GL01 dichotomisé, aucune des méthodes n'a détecté de gènes DE, même à un seuil FDR détendu de 20 %, ce qui souligne l'importance de modéliser la variable prédictive continue d'origine dans l'analyse DE.

Plusieurs sources de preuves soutiennent la validité biologique des gènes détectés par MACAU. Tout d'abord, nous avons effectué une analyse de l'ontologie génique (GO) à l'aide de LRpath (109) sur les gènes associés au T2D et à la GL identifiés par MACAU, comme dans l'étude FUSION (60) (Figure supplémentaire S17). Les résultats de l'analyse GO pour le DT2 et le GL sont cohérents avec les études précédentes (60, 110) et sont également similaires les uns aux autres, comme prévu étant donné la relation biologique entre les deux caractères. En particulier, le statut DT2 et une GL élevée sont associés à une expression réduite des gènes de la voie respiratoire cellulaire, ce qui est conforme aux observations précédentes (60, 110). Le statut T2D et la GL sont également associés à plusieurs voies liées à mTOR, y compris la génération de métabolites précurseurs, la poly-ubiquitination et le trafic de vésicules, en accord avec un rôle important de la voie mTOR dans l'étiologie du DT2 ( 111– 114).

Deuxièmement, nous avons effectué des analyses de chevauchement entre les gènes associés au T2D et au GL. Nous avons estimé que les gènes associés au DT2 sont probablement associés au GL parce que le DT2 partage une base génétique commune avec le GL (115-117) et que le statut de DT2 est déterminé en partie par le GL à jeun. Par conséquent, nous avons utilisé le chevauchement entre les gènes associés au DT2 et les gènes associés à la GL comme mesure de la performance de la méthode. Dans l'analyse de chevauchement, les gènes avec la plus forte association T2D identifiés par MACAU montrent un chevauchement plus important avec les 1000 premiers gènes qui ont la plus forte association GL que les gènes identifiés par d'autres méthodes (Figure 4D). Par exemple, parmi les 100 principaux gènes avec les preuves d'association DT2 les plus fortes de MACAU, 63 d'entre eux montrent également des preuves d'association fortes avec GL. En revanche, seuls 55 des 100 principaux gènes avec la plus forte association DT2 identifiés par GEMMA, la deuxième meilleure méthode, montrent une forte association avec GL. Nous avons observé des résultats similaires, MACAU effectuant les meilleurs résultats lors de l'analyse réciproque (chevauchement entre les gènes avec la plus forte association GL et les 1000 premiers gènes qui ont la plus forte association T2D : Figure 4F). Pour inclure la comparaison avec edgeR et DESeq2, nous avons examiné plus en détail le chevauchement entre les gènes associés au T2D et les gènes associés à GL01 pour toutes les méthodes (Figure supplémentaire S18). Encore une fois, MACAU est le plus performant, suivi de GEMMA et du modèle linéaire, et ni edgeR ni DESeq2 ne fonctionnent bien dans ce contexte (Figure supplémentaire S18). Par conséquent, MACAU semble à la fois conférer plus de pouvoir pour identifier les gènes DE biologiquement pertinents et être plus cohérent dans les analyses de phénotypes apparentés.

Pour évaluer le taux d'erreur de type I de diverses méthodes, nous avons permuté les données de traits du babouin et des études FUSION. En accord avec nos résultats de simulation, le P-les valeurs de MACAU et GEMMA sous le nul permuté étaient proches d'une distribution uniforme (légèrement conservatrice) dans les deux ensembles de données, alors que les autres méthodes ne l'étaient pas (figures supplémentaires S19 et 20). De plus, aucune des méthodes comparées ici n'est sensible aux valeurs aberrantes dans les deux ensembles de données (figures supplémentaires S21-23).

Enfin, bien que les grands ensembles de données RNAseq basés sur la population deviennent de plus en plus courants, le cadre de modélisation PMM flexible de MACAU permet de l'appliquer à l'analyse DE dans de petits ensembles de données avec des individus non apparentés également. Dans ce cadre, MACAU peut utiliser la matrice de covariance d'expression génique en tant que |$oldsymbol< K>$| matrice pour contrôler les effets de confusion cachés qui sont couramment observés dans les études de séquençage (48–51). Des facteurs de confusion cachés peuvent induire une similitude dans les niveaux d'expression des gènes à travers de nombreux gènes même si les individus ne sont pas apparentés (52-56), similaires aux effets de la parenté ou de la structure de la population. Par conséquent, en définissant |$oldsymbol< K>$| en utilisant une matrice de covariance d'expression génique (au lieu de génétique), MACAU peut contrôler efficacement la non-indépendance de l'échantillon induite par des facteurs de confusion cachés, étendant ainsi le LMM largement utilisé pour contrôler les facteurs de confusion cachés dans les études basées sur des matrices (52-56) au séquençage des données de comptage .

Pour illustrer cette application, nous avons analysé un troisième ensemble de données sur les LCL provenant de 69 individus nigérians non apparentés (YRI) (3) du projet HapMap (118), avec des mesures d'expression sur 13 319 gènes. Nous avons également cherché à identifier les gènes associés au sexe dans cet ensemble de données. Pour démontrer l'efficacité de MACAU dans de petits échantillons, nous avons sous-échantillonné au hasard des individus à partir des données pour créer de petits ensembles de données avec soit m = 6 (3 mâles et 3 femelles), m = 10 (5 mâles et 5 femelles) ou m = 14 individus (7 mâles et 7 femelles). Pour chaque taille d'échantillon m, nous avons effectué 20 répétitions de sous-échantillonnage aléatoire, puis évalué les performances de la méthode en faisant la moyenne entre les répétitions. Dans chaque réplicat, nous avons utilisé la matrice de covariance d'expression génique sous la forme |$oldsymbol$| et comparé les performances de MACAO par rapport à d'autres méthodes. En raison de la petite taille de l'échantillon, aucune des méthodes n'a été en mesure d'identifier les gènes DE à un seuil FDR de 10 %, ce qui est cohérent avec les arguments récents selon lesquels au moins 6 à 12 réplications biologiques sont nécessaires pour assurer une puissance et une réplicabilité suffisantes dans l'analyse DE ( 11 ). Nous avons donc utilisé l'enrichissement des gènes sur les chromosomes sexuels pour comparer les performances de différentes méthodes (Figure supplémentaire S24). L'enrichissement des gènes associés au sexe les mieux classés sur les chromosomes sexuels a déjà été utilisé pour la comparaison de méthodes et est particulièrement adapté pour comparer des méthodes en présence d'effets de lot et d'autres facteurs de confusion cachés (119).

Dans cette comparaison, MACAU exécute la meilleure de toutes les méthodes lorsque la taille de l'échantillon est soit m = 10 ou m = 14, et est classé parmi les meilleurs (avec le modèle binomial négatif) lorsque m = 6. Par exemple, lorsque m = 6, parmi les 50 premiers gènes identifiés par chaque méthode, le nombre de gènes sur les chromosomes sexuels pour MACAU, binomial négatif, Poisson, edgeR, DESeq2, GEMMA et Linear sont 3,3, 2,7, 3,1, 1,8, 3,0, 2,0 et 2,4 , respectivement. L'avantage de MACAU devient plus important lorsque la taille de l'échantillon augmente : par exemple, lorsque m = 14, une moyenne de 10,6 gènes dans les 50 premiers gènes de MACAU sont sur les chromosomes sexuels, ce qui est encore plus grand que celui du binôme négatif (8,3), Poisson (6,0), edgeR (6,65), DESeq2 (8,8), GEMMA (9,8) ou Linéaire (8,05). Ces résultats suggèrent que MACAU peut également être plus performant que les méthodes existantes dans des plans d'étude d'échantillons relativement petits avec des individus non apparentés en contrôlant les facteurs de confusion cachés. Cependant, le gain de puissance de MACAU est beaucoup plus faible dans ce cadre que dans les deux premiers ensembles de données que nous avons considérés (le babouin et les données Fusion). De plus, le gain de puissance de MACAU est négligeable dans le cas de m = 6 par rapport à la deuxième meilleure méthode, bien que son gain de puissance sur les edgeR et DESeq2 couramment utilisés soit encore substantiel. Le faible gain de puissance de MACAU dans ces données provient probablement à la fois de la petite taille de l'échantillon et de la petite taille de l'effet du sexe dans les données, conformément aux rapports précédents sur l'expression des gènes dérivés des cellules sanguines (3, 7, 120).


Remerciements

Un merci spécial à Leander Dony, qui a débogué, mis à jour et testé l'étude de cas pour travailler avec les dernières méthodes. De plus, nous tenons à remercier les nombreuses personnes qui ont relu le cahier d'étude de cas et le manuscrit et l'ont amélioré avec leurs commentaires et leur expertise. Pour cela, nous remercions Maren Buttner, David Fischer, Alex Wolf, Lukas Simon, Luis Ospina-Forero, Sophie Tritschler, Niklas Koehler, Goekcen Eraslan, Benjamin Schubert, Meromit Singer, Dana Pe'er et Rahul Satija.Un merci spécial pour cela également aux relecteurs anonymes du manuscrit et au rédacteur en chef, Thomas Lemberger, pour leurs commentaires approfondis, constructifs et détaillés. Le cahier d'études de cas a été testé et amélioré par les premiers utilisateurs Marius Lange, Hananeh Aliee, Subarna Palit et Lisa Thiergart. Volker Bergen et Alex Wolf ont également contribué au flux de travail en faisant des adaptations scanpy. Le choix de l'ensemble de données pour montrer de manière optimale tous les aspects du flux de travail d'analyse a été facilité par l'aimable contribution d'Adam Haber et d'Aviv Regev. Ce travail a été soutenu par la subvention BMBF # 01IS18036A et la subvention # 01IS18053A, par la German Research Foundation (DFG) au sein du Collaborative Research Center 1243, Subproject A17, par l'Association Helmholtz (Incubator grant sparse2big, grant # ZT-I-0007) et par la Chan Zuckerberg Initiative DAF (fonds conseillé de la Silicon Valley Community Foundation, 182835).


8.8 Comparaison/combinaison des ensembles de données scRNASeq

8.8.1 Présentation

À mesure que de plus en plus d'ensembles de données scRNA-seq deviennent disponibles, il est essentiel de réaliser des comparaisons merged_seurat entre eux. Il existe deux approches principales pour comparer les ensembles de données scRNASeq. La première approche est « centrée sur l'étiquette » qui vise à identifier des types/états de cellules équivalents dans des ensembles de données en comparant des cellules individuelles ou des groupes de cellules. L'autre approche est la « normalisation d'ensembles de données croisées » qui tente de supprimer par ordinateur les effets techniques/biologiques spécifiques à l'expérience afin que les données de plusieurs expériences puissent être combinées et analysées conjointement.

L'approche centrée sur l'étiquette peut être utilisée avec un ensemble de données avec des annotations de cellule à haute confiance, par ex. le Human Cell Atlas (HCA) (Regev et al. 2017) ou le Tabula Muris ( . ) une fois qu'ils sont terminés, pour projeter des cellules ou des amas d'un nouvel échantillon sur cette référence afin d'examiner la composition des tissus et/ou d'identifier les cellules avec une identité nouvelle/inconnue. Conceptuellement, de telles projections sont similaires à la méthode BLAST populaire (Altschul et al. 1990), qui permet de trouver rapidement la correspondance la plus proche dans une base de données pour une séquence de nucléotides ou d'acides aminés nouvellement identifiée. L'approche centrée sur l'étiquette peut également être utilisée pour comparer des ensembles de données d'origine biologique similaire collectés par différents laboratoires afin de garantir la cohérence de l'annotation et de l'analyse.

Figure 2.4 : La comparaison d'ensembles de données centrée sur l'étiquette peut être utilisée pour comparer les annotations de deux échantillons différents.

Figure 2.5 : La comparaison d'ensembles de données centrée sur l'étiquette peut projeter les cellules d'une nouvelle expérience sur une référence annotée.

L'approche de normalisation de l'ensemble de données croisées peut également être utilisée pour comparer des ensembles de données d'origine biologique similaire, contrairement à l'approche centrée sur l'étiquette, elle permet l'analyse conjointe de plusieurs ensembles de données pour faciliter l'identification de types de cellules rares qui peuvent être trop peu échantillonnés chez chaque individu. ensemble de données à détecter de manière fiable. Cependant, la normalisation inter-ensembles de données n'est pas applicable à des références très volumineuses et diverses car elle suppose qu'une partie importante de la variabilité biologique dans chacun des ensembles de données se chevauche avec d'autres.

Figure 2.6 : La normalisation des jeux de données croisés permet une analyse conjointe de 2+ jeux de données scRNASeq.

8.8.2 Ensembles de données

Nous exécuterons ces méthodes sur deux ensembles de données de pancréas humain : (Muraro et al. 2016) et (Segerstolpe et al. 2016). Comme le pancréas a été largement étudié, ces ensembles de données sont bien annotés.

Ces données ont déjà été formatées pour scmap. Les étiquettes de type de cellule doivent être stockées dans la colonne cell_type1 des emplacements colData, et les identifiants de gènes qui sont cohérents dans les deux ensembles de données doivent être stockés dans la colonne feature_symbol des emplacements rowData.

Tout d'abord, vérifions que nos identifiants génétiques correspondent aux deux ensembles de données :

Ici, nous pouvons voir que 96% des gènes présents dans les gènes de correspondance muraro dans segerstople et 72% des gènes dans segerstolpe sont des gènes correspondants dans muraro. C'est comme prévu car l'ensemble de données segerstolpe a été séquencé plus profondément que l'ensemble de données muraro. Cependant, il met en évidence certaines des difficultés liées à la comparaison des ensembles de données scRNASeq.

Nous pouvons le confirmer en vérifiant la taille globale de ces deux ensembles de données.

De plus, nous pouvons vérifier les annotations de type cellulaire pour chacun de ces jeux de données à l'aide de la commande ci-dessous :

Ici, nous pouvons voir que même si les deux ensembles de données considéraient le même tissu biologique, les deux ensembles de données, ils ont été annotés avec des ensembles de types cellulaires légèrement différents. Si vous connaissez la biologie du pancréas, vous reconnaîtrez peut-être que les cellules étoilées pancréatiques (CSP) du segerstolpe sont un type de cellule souche mésenchymateuse qui relèverait du type «mésenchymateux» du muraro. Cependant, il n'est pas clair si ces deux annotations doivent être considérées comme synonymes ou non. Nous pouvons utiliser des méthodes de comparaison centrée sur l'étiquette pour déterminer si ces deux annotations de type cellulaire sont effectivement équivalentes.

Alternativement, nous pourrions être intéressés à comprendre la fonction de ces cellules qui étaient « endocriniennes non classées » ou jugées de trop mauvaise qualité (« non applicable ») pour le regroupement d'origine dans chaque ensemble de données en tirant parti des informations entre les ensembles de données. Soit nous pourrions essayer de déduire à laquelle des annotations existantes ils appartiennent le plus probablement en utilisant des approches centrées sur les étiquettes, soit nous pourrions essayer de découvrir un nouveau type de cellule parmi eux (ou un sous-type dans les annotations existantes) en utilisant la normalisation de l'ensemble de données croisées .

Pour simplifier nos analyses de démonstration, nous supprimerons les petites classes de cellules non affectées et les cellules de mauvaise qualité. Nous retiendrons les « endocrines non classées » pour voir si l'une de ces méthodes peut élucider à quel type cellulaire elles appartiennent.

8.8.3 Projection de cellules sur des types de cellules annotés (scmap)

Nous avons récemment développé scmap (Kiselev et Hemberg 2017) - une méthode pour projeter des cellules à partir d'une expérience scRNA-seq sur les types cellulaires identifiés dans d'autres expériences. De plus, une version cloud de scmap peut être exécutée gratuitement, avec les restrictionsmerged_seurat, à partir de http://www.hemberg-lab.cloud/scmap.

8.8.3.1 Sélection des fonctionnalités

Une fois que nous avons un objet SingleCellExperiment, nous pouvons exécuter scmap . Nous devons d'abord construire l'« index » de nos clusters de référence. Puisque nous voulons savoir si les PSC et les cellules mésenchymateuses sont synonymes, nous projetterons chaque ensemble de données sur l'autre afin de construire un index pour chaque ensemble de données. Cela nécessite d'abord de sélectionner les caractéristiques les plus informatives pour l'ensemble de données de référence.

Les gènes mis en évidence avec la couleur rouge seront utilisés dans l'analyse ultérieure (projection).

À partir de l'axe des y de ces graphiques, nous pouvons voir que scmap utilise une méthode de sélection de caractéristiques basée sur dropmerged_seurat.

Calculez maintenant l'indice de type de cellule :

On peut aussi visualiser l'indice :

Vous voudrez peut-être ajuster vos caractéristiques à l'aide de la fonction setFeatures si les caractéristiques sont trop concentrées dans quelques types de cellules seulement. Dans ce cas, les fonctionnalités basées sur dropmerged_seurat semblent bonnes, nous ne les utiliserons donc que.

Exercer En utilisant les rowData de chaque ensemble de données, combien de gènes ont été sélectionnés comme caractéristiques dans les deux ensembles de données ? Qu'est-ce que cela vous dit amerged_seurat ces ensembles de données ?

8.8.3.2 Projection

scmap calcule la distance entre chaque cellule et chaque type de cellule dans l'index de référence, puis applique un seuil dérivé empiriquement pour déterminer quelles cellules sont affectées au type de cellule de référence le plus proche et lesquelles ne sont pas affectées. Pour tenir compte des différences de profondeur de séquençage, la distance est calculée à l'aide de la corrélation de lanceur et de la distance cosinusoïdale et seules les cellules avec une affectation cohérente avec les deux distances sont renvoyées telles qu'elles ont été attribuées.

Nous projetterons le jeu de données segerstolpe sur le jeu de données muraro :

et muraro sur segerstolpe

Notez que dans chaque cas, nous projetons sur un seul ensemble de données, mais que cela pourrait être étendu à n'importe quel nombre d'ensembles de données pour lesquels nous avons calculé des indices.

Comparons maintenant les étiquettes de type de cellule d'origine avec les étiquettes projetées :

Ici, nous pouvons voir que les types de cellules correspondent à leurs équivalents dans le segerstolpe, et surtout, nous voyons que toutes les cellules «mésenchymateuses» sauf une ont été affectées à la classe «PSC».

Encore une fois, nous voyons que les types de cellules correspondent et que toutes les « PSC » sauf une correspondent aux cellules « mésenchymateuses », fournissant une preuve solide que ces deux annotations doivent être considérées comme synonymes.

On peut aussi visualiser ces tableaux à l'aide d'un diagramme de Sankey :

Exercer Combien de cellules précédemment non classées pourraient être attribuées à des types de cellules à l'aide de scmap ?

8.8.4 Mappage de cellule à cellule

scmap peut également projeter chaque cellule d'un jeu de données sur sa cellule voisine la plus proche approximative dans le jeu de données de référence. Celui-ci utilise un algorithme de recherche hautement optimisé lui permettant d'être mis à l'échelle à de très grandes références (en théorie 100 000-millions de cellules). Cependant, ce processus est stochastique, nous devons donc corriger la graine aléatoire pour nous assurer que nous pouvons reproduire nos résultats.

Nous avons déjà effectué une sélection d'entités pour cet ensemble de données afin que nous puissions passer directement à la création de l'index.

Dans ce cas, l'index est une série de clusters de chaque cellule utilisant différents ensembles de caractéristiques, les paramètres k et M sont le nombre de clusters et le nombre de caractéristiques utilisées dans chacun de ces sous-clusters. De nouvelles cellules sont affectées au cluster le plus proche dans chaque sous-cluster pour générer un modèle unique d'affectations de cluster. Nous trouvons ensuite la cellule dans l'ensemble de données de référence avec le même modèle ou le plus similaire d'affectations de cluster.

Nous pouvons examiner les modèles d'attribution de cluster pour les ensembles de données de référence en utilisant :

Pour projeter et trouver les w voisins les plus proches, nous utilisons une commande similaire à celle précédente :

On peut encore regarder les résultats :

Cela montre le numéro de colonne des 5 voisins les plus proches en segerstolpe de chacune des cellules de muraro. Nous pourrions ensuite calculer une estimation de pseudo-temps, une affectation de branche ou d'autres données au niveau de la cellule en sélectionnant les données appropriées dans le colData de l'ensemble de données segerstolpe. A titre de démonstration, nous allons trouver le type de cellule du plus proche voisin de chaque cellule.

8.8.5 Métavoisin

Metaneighbor est spécifiquement conçu pour demander si les étiquettes de type de cellule sont cohérentes entre les ensembles de données. Il se décline en deux versions. La première est une méthode entièrement supervisée qui suppose que les types de cellules sont connus dans tous les ensembles de données et calcule à quel point ces étiquettes de types de cellules sont « bonnes ». (La signification précise de « bon » sera décrite ci-dessous). Alternativement, le métavoisin peut estimer à quel point tous les types de cellules sont similaires les uns aux autres à la fois dans et entre les ensembles de données. Nous n'utiliserons que la version non supervisée car elle a une applicabilité beaucoup plus générale et est plus facile à interpréter les résultats.

Metaneighbor compare les types de cellules à travers des ensembles de données en créant un réseau de corrélation cellule-cellule Spearman. La méthode essaie ensuite de prédire l'étiquette de chaque cellule grâce aux « votes » pondérés de ses plus proches voisins. Note ensuite la similitude globale entre deux clusters en tant qu'AUROC pour l'attribution des cellules de type A à type B en fonction de ces votes pondérés. Un AUROC de 1 indiquerait que toutes les cellules de type A ont été affectées au type B avant toutes les autres cellules, et un AUROC de 0,5 est ce que vous obtiendriez si les cellules étaient affectées au hasard.

Metanighbour n'est que quelques fonctions R et non un package complet, nous devons donc les charger en utilisant la source

8.8.5.1 Préparer les données

Metavoisin nécessite que tous les ensembles de données soient combinés dans une seule matrice d'expression avant de s'exécuter :

Metaneighbor inclut une méthode de sélection de caractéristiques pour identifier des gènes très variables.

Étant donné que Metaneighbor est beaucoup plus lent que scmap , nous allons sous-échantillonner ces ensembles de données.

Nous sommes maintenant prêts à exécuter Metaneighbor. Nous allons d'abord exécuter la version non supervisée qui nous permettra de voir quels types de cellules sont les plus similaires dans les deux ensembles de données.

8.8.6 mnnCorrect

mnnCorrect corrige les ensembles de données pour faciliter l'analyse conjointe. Afin de tenir compte des différences de composition entre deux répétitions ou deux expériences différentes, il fait d'abord correspondre les cellules individuelles à travers les expériences pour trouver la structure biologique qui se chevauche. En utilisant ce chevauchement, il apprend quelles dimensions d'expression correspondent à l'état biologique et quelles dimensions correspondent à l'effet lot/expérience. Enfin, il supprime les effets de lot/d'expérience de toute la matrice d'expression pour renvoyer la matrice corrigée.

Pour faire correspondre les cellules individuelles entre les jeux de données, mnnCorrect utilise la distance en cosinus pour éviter l'effet de la taille de la bibliothèque, puis identifie les voisins mutuels les plus proches (k détermine la taille du quartier) dans les jeux de données. Seuls les groupes biologiques qui se chevauchent devraient avoir des voisins mutuels les plus proches (voir panneau b ci-dessous). Cependant, cela suppose que k est défini approximativement à la taille du plus petit groupe biologique dans les ensembles de données, mais un k trop faible identifiera trop peu de paires mutuelles de voisins les plus proches pour obtenir une bonne estimation de l'effet de lot que nous voulons supprimer.

L'apprentissage des effets biologiques/techniques se fait soit avec une décomposition en valeur singulière, similaire à RUV que nous rencontrons dans la section de correction par lots, soit avec une analyse en composantes principales avec le package irlba optimisé, qui devrait être plus rapide que SVD. Le paramètre svd.dim spécifie combien de dimensions doivent être conservées pour résumer la structure biologique des données, nous le définirons sur trois car nous avons trouvé trois groupes principaux en utilisant Metaneighbor ci-dessus. Ces estimations peuvent être encore ajustées par lissage ( sigma ) et/ou ajustement de la variance ( var.adj ).

mnnCorrect suppose également que vous avez déjà sous-ensemble vos matrices d'expression afin qu'elles contiennent des gènes identiques dans le même ordre, heureusement, nous en avons déjà fait pour nos ensembles de données lorsque nous avons configuré nos données pour Metaneighbor.

Figure 8.20 : mnnCorrection de l'effet batch/dataset. De Haghverdi et al. 2017

Vérifions d'abord que nous avons trouvé un nombre suffisant de paires mnn, mnnCorrect renvoie une liste de dataframe avec les paires mnn pour chaque ensemble de données.

Les première et deuxième colonnes contiennent les identifiants des colonnes de cellules et la troisième colonne contient un nombre indiquant à quel ensemble de données/lot la cellule de la colonne 2 appartient. Dans notre cas, nous ne comparons que deux ensembles de données, donc toutes les paires mnn ont été affectées à la deuxième table et la troisième colonne n'en contient qu'une.

mnnCorrect a trouvé 2533 ensembles de plus proches voisins mutuels entre les cellules n_unique_seger segerstolpe et les cellules n_unique_muraro muraro. Cela devrait être un nombre suffisant de paires, mais le faible nombre de cellules uniques dans chaque ensemble de données suggère que nous n'avons peut-être pas capturé le signal biologique complet dans chaque ensemble de données.

Exercer Quels types de cellules avaient des mnns dans ces ensembles de données ? Doit-on augmenter/diminuer k ?

Nous pouvons maintenant créer un ensemble de données combiné pour analyser conjointement ces données. Cependant, les données corrigées ne sont plus dénombrées et contiendront généralement des valeurs d'expression négatives, de sorte que certains outils d'analyse peuvent ne plus être appropriés. Pour plus de simplicité, traçons simplement un TSNE conjoint.

8.8.7 Analyse de corrélation canonique (Seurat)

Le package Seurat contient une autre méthode de correction pour combiner plusieurs ensembles de données, appelée CCA. Cependant, contrairement à mnnCorrect, il ne corrige pas directement la matrice d'expression elle-même. Au lieu de cela, Seurat trouve un sous-espace de dimension inférieure pour chaque ensemble de données, puis corrige ces sous-espaces. Également différent de mnnCorrect, Seurat ne combine qu'une seule paire de jeux de données à la fois.

Seurat utilise des corrélations gène-gène pour identifier la structure biologique dans l'ensemble de données avec une méthode appelée analyse de corrélation canonique (ACC). Seurat apprend la structure partagée des corrélations gène-gène et évalue ensuite dans quelle mesure chaque cellule s'adapte à cette structure. Les cellules qui doivent être mieux décrites par une méthode de réduction de dimensionnalité spécifique aux données que par la structure de corrélation partagée sont supposées représenter des types/états de cellules spécifiques à un ensemble de données et sont rejetées avant d'aligner les deux ensembles de données. Enfin, les deux ensembles de données sont alignés à l'aide d'algorithmes de « warping » qui normalisent les représentations de faible dimension de chaque ensemble de données d'une manière robuste aux différences de densité de population.

Notez que Seurat utilise beaucoup d'espace de bibliothèque, vous devrez redémarrer votre R-session pour le charger, et les plots/merged_seuratput ne seront pas générés automatiquement sur cette page.

Tout d'abord, nous allons reformater nos données en objets Seurat :

Ensuite, nous devons normaliser, mettre à l'échelle et identifier des gènes très variables pour chaque ensemble de données :

Figure 8.21 : gènes variables muraro

Figure 8.22 : gènes variables de segerstolpe

Même si Seurat corrige la relation entre la dispersion et l'expression moyenne, il n'utilise pas la valeur corrigée lors du classement des entités. Comparez les résultats de la commande ci-dessous avec les résultats des graphiques ci-dessus :

Mais nous suivrons leur exemple et utiliserons les 2000 gènes les plus dispersés avecmerged_seurat corrigeant de toute façon l'expression moyenne de chaque ensemble de données.

Exercer Trouvez les fonctionnalités que nous utiliserions si nous sélectionnions le top 2000 les plus dispersés après mise à l'échelle par la moyenne. (Indice : considérez la fonction de commande)

Nous allons maintenant exécuter CCA pour trouver la structure de corrélation partagée pour ces deux ensembles de données :

Notez que pour accélérer les calculs, nous n'utiliserons que les 5 premières dimensions, mais idéalement, vous en envisageriez beaucoup plus, puis sélectionneriez les plus informatives à l'aide de DimHeatmap .

Figure 8.23 ​​: Avant l'alignement

Pour identifier les types de cellules spécifiques à l'ensemble de données, nous comparons la manière dont les cellules sont « expliquées » par l'analyse CCA par rapport à l'analyse en composantes principales spécifiques à l'ensemble de données.

Ici, nous pouvons voir que malgré les deux ensembles de données contenant des cellules endothéliales, presque tous ont été rejetés comme « spécifiques à l'ensemble de données ». Nous pouvons maintenant aligner les ensembles de données :

Figure 8.24 : Après l'alignement

Exercer Comparez les résultats si vous utilisez les fonctionnalités après avoir mis à l'échelle les dispersions.

Figure 8.25 : Après l'alignement

Exercice avancé Utilisez les méthodes de clustering que nous avons décrites précédemment sur les ensembles de données combinés. Identifiez-vous de nouveaux types cellulaires ?

8.8.8 SessionInfo()


Remerciements

Nous remercions Toshi Kinoshita (Department of Pathology, University of Wisconsin School of Medicine and Public Heath, Madison, WI) pour son aide en histologie.

Contributions d'auteur

N.V.W. conçu l'étude et obtenu un financement. N.V.W., S.L.T. et C.K. conçu les expériences. N.V.W. et mon. a mené le in vivo expériences. C.L. traité tous les échantillons d'ARN et de tissus et effectué une immunohistochimie. J.A.D. effectué toutes les analyses statistiques sous la direction de C.K. N.V.W. écrit le manuscrit. Tous les auteurs ont revu et approuvé la version finale.

Ce travail a été financé par des subventions du National Institute on Deafness and other Communication Disorders [subventions R01 DC004428, R01 DC010777] et le programme Clinical and Translational Science Award (CTSA) du National Center for Research Resources [subvention numéro UL1 RR025011].


Voir la vidéo: Scroll The Page u0026 Earn $ Again u0026 Again! - Make Money Online. Branson Tay (Janvier 2022).