Une étude unicellulaire de l’épithélium de l’intestin grêle

Souris

Tous les travaux sur les souris ont été réalisés conformément aux comités institutionnels de soins et d’utilisation des animaux (IACUC) et aux directives pertinentes du Broad Institute et du Massachusetts Institute of Technology, avec les protocoles 0055-05-15 et 0612-058-18, respectivement. Pour toutes les expériences, les souris ont été réparties au hasard dans les groupes de traitement après avoir été appariées en fonction du sexe et de l’âge à des souris de type sauvage C57BL/6J ou Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP) âgées de 7 à 10 semaines, femelles ou mâles, obtenues auprès du Jackson Laboratory (Bar Harbour) ou à des souris Gfi1beGFP/+ (Gfi1b-GFP)43. Les souris ont été hébergées dans des conditions exemptes de pathogènes spécifiques dans les animaleries du Broad Institute, du Massachusetts Institute of Technology ou de la Harvard T. H. Chan School of Public Health.

Infection à Salmonella enterica et H. polygyrus. Des souris C57BL/6J (Jackson Laboratory) ont été infectées avec 200 larves de troisième stade de H. polygyrus ou 108 Salmonella enterica, maintenues dans des conditions exemptes de pathogènes spécifiques au Massachusetts General Hospital (Charlestown), avec le protocole 2003N000158. H. polygyrus a été propagé comme décrit précédemment44. Les souris ont été euthanasiées 3 et 10 jours après l’infection par H. polygyrus. Pour Salmonella enterica, les souris ont été infectées avec une souche SL1344 de S. Typhimurium naturellement résistante à la streptomycine (108 cellules) comme décrit précédemment44, et ont été euthanasiées 48 h après l’infection.

Dissociation cellulaire et isolement des cryptes

Isolation des cryptes. L’intestin grêle de souris C57BL/6J de type sauvage, Lgr5-GFP ou Gfi1b-GFP a été isolé et rincé dans du PBS froid. Le tissu a été ouvert longitudinalement et tranché en petits fragments d’environ 2 mm de long. Le tissu a été incubé dans du PBS EDTA 20 mM sur glace pendant 90 minutes, en le secouant toutes les 30 minutes. Le tissu a ensuite été secoué vigoureusement et le surnageant a été recueilli comme fraction 1 dans un nouveau tube conique. Le tissu a été incubé dans du EDTA-PBS frais et une nouvelle fraction a été recueillie toutes les 30 minutes. Les fractions ont été collectées jusqu’à ce que le surnageant soit presque entièrement constitué de cryptes. La fraction finale (enrichie en cryptes) a été lavée deux fois dans du PBS, centrifugée à 300g pendant 3 min, et dissociée avec TrypLE express (Invitrogen) pendant 1 min à 37 °C. La suspension unicellulaire a ensuite été passée à travers un filtre de 40-μm et colorée pour FACS pour scRNA-seq (ci-dessous) ou utilisée pour la culture d’organoïdes. Nous avons confirmé la robustesse de cette méthode en testant d’autres méthodes d’isolement de cellules uniques – soit  » entières  » (grattage de la muqueuse épithéliale), soit  » enrichies en villosités  » (fraction 1 ; voir ci-dessus) – et nous avons constaté qu’en raison du taux de mortalité élevé (via l’anoikis) des cellules différenciées post-mitotiques (dont le composant principal est les entérocytes matures), la suspension de cellules uniques enrichie en crypte représente fidèlement la composition des types de cellules de l’intestin grêle (données non présentées).

Isolation des épithéliums associés aux follicules. Les cellules épithéliales de l’épithélium associé aux follicules ont été isolées en extrayant de petites sections (0,2-0,5 cm) contenant les plaques de Peyer de l’intestin grêle de souris C57Bl/6J ou Gfi1beGFP/+.

Tri des cellules

Pour les expériences scRNA-seq pleine longueur sur plaque, une machine FACS (Astrios) a été utilisée pour trier une seule cellule dans chaque puits d’une plaque PCR 96 puits contenant 5 μl de tampon TCL avec 1% de 2-mercaptoéthanol. Pour isoler les EpCAM+, les cellules ont été colorées pour 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience) ; pour les cellules épithéliales spécifiques, nous avons également coloré pour CD24+/- (eBioscience) et c-Kit+/- (eBioscience). Pour enrichir des populations spécifiques de cellules épithéliales intestinales, les cellules ont été isolées à partir de souris Lgr5-GFP, colorées avec les anticorps mentionnés ci-dessus et sélectionnées sur GFP-high (cellules souches), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (lignées sécrétoires) ou GFP-/CD24-/EpCAM+ (cellules épithéliales). Pour une meilleure récupération des cellules de Paneth, nous avons autorisé des paramètres de diffusion latérale et de diffusion avant plus élevés en combinaison avec CD24+/c-Kit+ pour vérifier la récupération des cellules de Paneth dans les cellules EpCAM+. Pour l’isolement des cellules tuft-2, les cellules épithéliales de trois souris différentes ont été colorées comme ci-dessus, mais en utilisant EpCAM+/CD45+ pour trier 2 000 cellules uniques. Nous avons utilisé une porte de tri indulgente pour nous assurer d’obtenir un nombre suffisant de ces rares cellules tuft-2, ce qui a conduit à un taux de contamination plus élevé de cellules T, que nous avons éliminées dans notre analyse unicellulaire à l’aide d’un clustering non supervisé.

Pour le tri scRNA-seq complet, la plaque à 96 puits a été fermée hermétiquement avec un Microseal F et centrifugée à 800g pendant 1 min. La plaque a été immédiatement congelée sur de la glace sèche et conservée à -80 °C jusqu’à ce qu’elle soit prête pour le nettoyage du lysat. Les cellules de la population en vrac ont été triées dans un tube Eppendorf contenant 100 μl de solution de TCL avec 1% de 2-mercaptoéthanol et conservées à -80 °C.

Pour le scRNA-seq basé sur les gouttelettes, les cellules ont été triées avec les mêmes paramètres que pour le scRNA-seq basé sur les plaques, mais elles ont été triées dans un tube Eppendorf contenant 50 μl de BSA-PBS à 0,4 % et stockées sur la glace jusqu’à ce qu’elles passent à la plateforme de cellules uniques GemCode.

ScRNA-seq basé sur les plaques

Cellules uniques. Librairies ont été préparés en utilisant un protocole SMART-Seq2 modifié16. En bref, le nettoyage du lysat d’ARN a été effectué à l’aide de billes RNAClean XP (Agencourt), suivi d’une transcription inverse avec la transcriptase inverse Maxima (Life Technologies) et d’une amplification par transcription complète (WTA) avec KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) pendant 21 cycles. Les produits de l’ATA ont été purifiés avec des billes Ampure XP (Beckman Coulter), quantifiés avec le Qubit dsDNA HS Assay Kit (ThermoFisher) et évalués avec une puce à ADN haute sensibilité (Agilent). Les bibliothèques RNA-seq ont été construites à partir des produits purifiés de l’ATA à l’aide du kit de préparation de bibliothèques d’ADN Nextera XT (Illumina). Sur chaque plaque, la population et les contrôles sans cellule ont été traités en utilisant la même méthode que pour les cellules uniques. Les bibliothèques ont été séquencées sur un NextSeq 500 d’Illumina.

Échantillons en vrac. Les échantillons de population en vrac ont été traités en extrayant l’ARN avec le kit RNeasy Plus Micro (Qiagen) selon les recommandations du fabricant, puis en procédant au protocole SMART-Seq2 modifié après le nettoyage du lysat, comme décrit ci-dessus.

Droplet-based scRNA-seq

Les cellules uniques ont été traitées par la plateforme de cellules uniques GemCode en utilisant les kits GemCode Gel Bead, Chip et Library (10X Genomics, Pleasanton) selon le protocole du fabricant. En bref, les cellules uniques ont été triées dans du BSA-PBS à 0,4 %. 6 000 cellules ont été ajoutées à chaque canal avec un taux de récupération moyen de 1 500 cellules. Les cellules ont ensuite été réparties dans des billes de gel en émulsion dans l’instrument GemCode, où ont eu lieu la lyse cellulaire et la transcription inverse de l’ARN avec code-barres, puis l’amplification, le cisaillement et la fixation de l’adaptateur 5′ et de l’index d’échantillon. Les librairies ont été séquencées sur un Illumina NextSeq 500.

Immunofluorescence et smFISH

Immunofluorescence. La coloration des tissus de l’intestin grêle a été réalisée comme décrit précédemment34. En bref, les tissus ont été fixés pendant 14 h dans du formol, inclus dans de la paraffine et coupés en sections de 5-μm d’épaisseur. Les sections ont été déparaffinées à l’aide de techniques standard, incubées avec des anticorps primaires pendant une nuit à 4 °C, puis avec des anticorps secondaires à température ambiante pendant 30 min. Les lames ont été montées avec le Slowfade Mountant + DAPI (Life Technologies, S36964) et scellées.

smFISH. Un kit de florescence RNAScope Multiplex (Advanced Cell Diagnostics) a été utilisé selon les recommandations du fabricant avec les modifications suivantes. Le temps d’ébullition pour la récupération de la cible a été ajusté à 12 min et l’incubation avec la protéase IV à 40 °C a été ajustée à 8 min. Les lames ont été montées avec le Slowfade Mountant+DAPI (Life Technologies, S36964) et scellées.

Immunofluorescence combinée et smFISH. Ceci a été mis en œuvre en effectuant d’abord le smFISH comme décrit ci-dessus, avec les changements suivants. Après l’Amp 4, les sections de tissu ont été lavées dans un tampon de lavage, incubées avec des anticorps primaires pendant une nuit à 4 °C, lavées dans 1× TBST trois fois, puis incubées avec des anticorps secondaires pendant 30 min à température ambiante. Les lames ont été montées avec le Slowfade Mountant + DAPI (Life Technologies, S36964) et scellées.

Analyse des images

Les images des sections de tissus ont été prises avec un microscope confocal Fluorview FV1200 en utilisant Kalman et l’émission laser séquentielle pour réduire le bruit et le chevauchement des signaux. Des barres d’échelle ont été ajoutées à chaque image à l’aide du logiciel confocal FV10-ASW 3.1 Viewer. Les images ont été superposées et visualisées à l’aide du logiciel Image J45.

Anticorps et sondes

Cultures d’organoïdes intestinaux

Après isolement des cryptes, la suspension unicellulaire a été remise en suspension dans du Matrigel (BD Bioscience) avec 1 μM de peptide Jagged-1 (Ana-Spec). Environ 300 cryptes incluses dans 25 μl de Matrigel ont été ensemencées dans chaque puits d’une plaque à 24 puits. Une fois solidifié, le Matrigel a été incubé dans 600 μl de milieu de culture (Advanced DMEM/F12, Invitrogen) avec streptomycine/pénicilline et glutamatax et supplémenté avec EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), Y-276432 dihydrochlorure monohydraté (10 μM, Tochris), N-acétyl-1-cystéine (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) et Wnt3A (25 ng ml-1, R&D Systems). Le milieu frais a été remplacé au jour 3, et les organoïdes ont été passés par dissociation avec TrypLE et remis en suspension dans du nouveau Matrigel au jour 6 avec un rapport de division 1:3. Pour certaines expériences, les organoïdes ont été traités en plus avec du RANKL (100 ng ml-1, Biolegends). Les organoïdes traités ont été dissociés et soumis à scRNA-seq en utilisant les deux méthodes.

PCR quantitative

L’ADNc de 16 cellules uniques de tuft-1, tuft-2 et EpCam+ aléatoires, amplifiées par transcriptome entier, provenant des plaques scRNA-seq basées sur la longueur totale, a été utilisé pour la qPCR relative. L’expression des gènes a été analysée par PCR quantitative en temps réel sur un instrument II LightCycler 480 (Roche) en utilisant le mélange vert SYBR LightCycler 480 (Roche) avec les ensembles d’amorces suivants : HPRT1-F, GTTAAGCAGTACAGCCCCAAA ; HPRT1-R, AGGGCATATCCAACAACAAACTT ; UBC-F, CAGCCGTATCTTCCCAGACT ; UBC-R, CTCAGAGGGATGCCAGTAATCTA ; tslp-F, TACTCTCAATCCTATCCCTGGCTG ; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC ; Alpi-F, TCCTACACCTCCATTCTCTATGG, Alpi-R, CCGCCTGCTGCTTGTAG ; Dclk1-F, GGGTGAGAACCATCTACCATC ; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. Les amorces qPCR ont été conçues pour une limite exon-exon dans tous les transcrits.

Analyse computationnelle

Pré-traitement des données scRNA-seq basées sur les gouttelettes. Le démultiplexage, l’alignement sur le transcriptome mm10 et le collapsage par identifiant moléculaire unique (UMI) ont été réalisés à l’aide de la boîte à outils Cellranger (version 1.0.1) fournie par 10X Genomics. Pour chaque cellule, nous avons quantifié le nombre de gènes pour lesquels au moins une lecture était cartographiée, puis nous avons exclu toutes les cellules ayant moins de 800 gènes détectés. Les valeurs d’expression Ei,j pour le gène i dans la cellule j ont été calculées en divisant les comptes UMI pour le gène i par la somme des comptes UMI dans la cellule j, afin de normaliser les différences de couverture, puis en multipliant par 10 000 pour créer des valeurs de type TPM, et enfin en calculant log2(TPM + 1). La correction par lots a été effectuée à l’aide de ComBat46 tel qu’implémenté dans le paquet R sva47, en utilisant le mode d’ajustement paramétrique par défaut. La sortie était une matrice d’expression corrigée, qui a été utilisée comme entrée pour les analyses ultérieures.

La sélection des gènes variables a été effectuée en ajustant un modèle linéaire généralisé à la relation entre le coefficient de variation au carré et le niveau d’expression moyen dans l’espace logarithmique, et en sélectionnant les gènes qui s’écartaient significativement (P < 0,05) de la courbe ajustée48.

Pré-traitement des données SMART-Seq2 scRNA-seq. Les fichiers BAM ont été convertis en FASTQ fusionnés et démultiplexés à l’aide du progiciel Bcl2Fastq v2.17.1.14 fourni par Illumina. Les lectures en paires ont été mappées sur le transcriptome de la souris UCSC mm10 à l’aide de Bowtie49 avec les paramètres ‘-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6’, qui permet l’alignement des séquences avec un mismatch. Les niveaux d’expression des gènes ont été quantifiés en utilisant les valeurs TPM calculées par RSEM50 v1.2.3 en mode paired-end. Pour chaque cellule, nous avons quantifié le nombre de gènes pour lesquels au moins une lecture était cartographiée, puis nous avons exclu toutes les cellules pour lesquelles moins de 3 000 gènes avaient été détectés ou dont la cartographie du transcriptome était inférieure à 40 %. Nous avons ensuite identifié les gènes hautement variables comme décrit ci-dessus.

Réduction de la dimensionnalité en utilisant PCA et t-SNE. Nous avons restreint la matrice d’expression aux sous-ensembles de gènes variables et aux cellules de haute qualité notés ci-dessus, puis nous avons centré et mis à l’échelle les valeurs avant de les introduire dans l’analyse en composantes principales (ACP), qui a été mise en œuvre à l’aide de la fonction R prcomp du paquet stats pour le jeu de données SMART-seq2. Pour l’ensemble de données basé sur les gouttelettes, nous avons utilisé une approximation aléatoire de l’ACP, mise en œuvre à l’aide de la fonction rpca du paquet R rsvd, avec le paramètre k fixé à 100. Cette approximation à faible rang a été utilisée car elle est plusieurs ordres de grandeur plus rapide à calculer pour des matrices très larges. Étant donné que de nombreuses composantes principales n’expliquent que très peu de la variance, le rapport signal/bruit peut être amélioré de manière substantielle en sélectionnant un sous-ensemble de n composantes principales « significatives ». Après l’ACP, les composantes principales significatives ont été identifiées à l’aide du test de permutation51, mis en œuvre à l’aide de la fonction permutationPA du paquet R jackstraw. Ce test a identifié 13 et 15 composantes principales significatives dans les ensembles de données 10X et SMART-Seq2 de la figure 1b et des données étendues de la figure 2a, respectivement. Les scores de ces seules composantes principales significatives ont été utilisés comme entrée pour une analyse plus poussée.

Pour la visualisation, la dimensionnalité des ensembles de données a été réduite davantage en utilisant la version approximative  » Barnes-hut  » de t-SNE52,53. Cela a été mis en œuvre à l’aide de la fonction Rtsne du paquet Rtsne en utilisant 20 000 itérations et un paramètre de perplexité qui variait de 10 à 30 selon la taille de l’ensemble de données.

Identification des trajectoires de différenciation cellulaire à l’aide de cartes de diffusion

Avant d’exécuter la réduction de la dimensionnalité des cartes de diffusion, nous avons sélectionné les gènes hautement variables dans les données comme suit. Nous avons d’abord ajusté un modèle nul pour la variabilité de base de l’expression des gènes entre les cellules dans les données, en utilisant une relation de puissance-loi entre le coefficient de variation et la moyenne des comptes UMI de tous les gènes exprimés, similaire à des travaux antérieurs54. Ensuite, nous avons calculé pour chaque gène la différence entre la valeur de son coefficient de variation observé et celle attendue par le modèle nul (CVdiff). L’histogramme du CVdiff présentait une queue « grasse ». Nous avons calculé la moyenne μ et l’écart-type σ de cette distribution, et avons sélectionné tous les gènes pour lesquels le CVdiff > μ + 1,67σ, ce qui a donné 761 gènes pour une analyse plus approfondie.

Nous avons effectué une réduction de la dimensionnalité en utilisant l’approche de la carte de diffusion22. En bref, une matrice de transition cellule-cellule a été calculée en utilisant un noyau gaussien, avec la largeur du noyau ajustée au voisinage local de chaque cellule55. Cette matrice a été convertie en une matrice markovienne après normalisation. Les vecteurs propres de droite vi (i = 0, 1, 2, …) de cette matrice ont été calculés et triés par ordre décroissant de valeur propre λi (i = 0, 1, 2, …), après exclusion du vecteur propre  » supérieur  » v0, correspondant à λ0 = 1 (ce qui reflète la contrainte de normalisation de la matrice markovienne). Les autres vecteurs propres vi (i = 1, 2, …) définissent l’intégration de la carte de diffusion et sont appelés composantes de diffusion (DCk, k = 1, 2, …). Nous avons remarqué un écart spectral entre λ4 et λ5, et avons donc retenu DC1-DC4 à la fois pour l’ensemble de données initial (Données étendues Fig. 4) et pour les données extraites de régions intestinales distinctes (Fig. 2c).

Enlever les cellules immunitaires contaminantes et les doublets

Bien que les cellules aient été triées avant le séquençage à l’aide d’EpCAM, un petit nombre de cellules immunitaires contaminantes a été observé dans l’ensemble de données 10X. Ces 264 cellules ont été éliminées par un premier tour de clustering non supervisé (clustering basé sur la densité de la carte t-SNE à l’aide de dbscan56 du package R fpc) car elles formaient un cluster extrêmement distinct. Pour l’ensemble de données SMART-Seq2, plusieurs cellules étaient aberrantes en termes de complexité de bibliothèque, ce qui pouvait éventuellement correspondre à plus d’une cellule individuelle par bibliothèque de séquençage (« doublets »). Ces cellules ont ensuite été supprimées en calculant le quantile supérieur de 1 % de la distribution des gènes détectés par cellule et en supprimant toutes les cellules dans ce quantile.

Analyse de cluster

Pour regrouper les cellules individuelles par leur expression, nous avons utilisé une approche de clustering non supervisée, basée sur l’algorithme de clustering de graphes Infomap9, suivant les approches pour les données CyTOF de cellules individuelles57 et scRNA-seq10. En bref, nous avons construit un graphe des k plus proches voisins sur les données en utilisant, pour chaque paire de cellules, la distance euclidienne entre les scores des composantes principales significatives pour identifier les k plus proches voisins. Le paramètre k a été choisi pour être cohérent avec la taille de l’ensemble de données. Plus précisément, k a été fixé à 200 et 80 pour l’ensemble de données basé sur les gouttelettes de 7 216 cellules (Fig. 1b) et pour l’ensemble de données SMART-Seq2 de 1 522 cellules (Données étendues Fig. 2a), respectivement. Les organoïdes traités au RANKL contenaient 5 434 cellules et k a été fixé à 200 ; l’ensemble de données sur Salmonella et H. polygyrus contenait 9 842 cellules et k a été fixé à 500. Pour les analyses de regroupement au sein des types de cellules, en particulier les sous-ensembles de cellules entéroendocrines et de cellules de touffes, nous avons utilisé la distance de corrélation de Pearson au lieu de la distance euclidienne, et avons fixé k = 15, k = 30 et k = 40 pour les sous-types entéroendocrines (533 cellules), et pour les 166 et 102 cellules de touffes dans les ensembles de données 10X et SMART-Seq2, respectivement. Le graphique des plus proches voisins a été calculé à l’aide de la fonction nng du paquet R cccd. Le graphe des k plus proches voisins a ensuite été utilisé comme entrée pour Infomap9, mis en œuvre à l’aide de la fonction infomap.community du package R igraph.

Les clusters détectés ont été cartographiés sur des types de cellules ou des états intermédiaires en utilisant des marqueurs connus pour les sous-types de cellules épithéliales intestinales. (Données étendues Fig. 1g, Données étendues Fig. 2a). Pour la sous-analyse des cellules entéroendocrines (EEC) (Fig. 3), tout groupe de clusters de progéniteurs EEC présentant des corrélations moyennes par paire entre les scores significatifs des composantes principales de r > 0,85 a été fusionné, ce qui a donné lieu à quatre clusters. Nous avons étiqueté ces quatre grappes progéniteur « A » sur la base de niveaux élevés de Ghrl, ou progéniteur (précoce), (moyen) ou (tardif) (dans cet ordre) sur la base de niveaux décroissants de gènes de tige (Slc12a2, Ascl2, Axin2) et de cycle cellulaire et de niveaux croissants de facteurs régulateurs connus de la CEE (Neurod1, Neurod2 et Neurog3) (Données étendues Fig. 5c, Tableau supplémentaire 6). Pour l’ensemble de données SMART-Seq2, deux clusters exprimant des niveaux élevés de gènes marqueurs de cellules souches (Extended Data Fig. 2a) ont été fusionnés pour former un cluster  » stem  » et deux autres clusters ont été fusionnés pour former un cluster  » TA « .

Pour l’analyse de cluster de l’ensemble de données de l’épithélium associé aux follicules de 4 700 cellules, les cellules microfold étaient très rares (0,38 %) et la méthode ClusterDP58 a donc été utilisée pour les identifier car elle a donné de meilleurs résultats empiriques que l’algorithme de graphe des plus proches voisins k sur cet ensemble de données. Comme pour les méthodes de k-plus proche voisin, ClusterDP a été exécuté en utilisant des scores de composantes principales significatives (P < 0,05) (19 dans ce cas) comme entrée, et a été mis en œuvre en utilisant les fonctions findClusters et densityClust du paquet R densityClust en utilisant les paramètres rho = 1.1 et delta = 0,25.

Extraction des types de cellules rares pour une analyse plus approfondie

Le regroupement initial de l’ensemble de données de l’intestin entier (7 216 cellules ; figure 1b) a montré un regroupement de 310 cellules EEC et 166 cellules de touffe. Les cellules de touffe ont été prises  » telles quelles  » pour la sous-analyse (Fig. 4a, b), tandis que les cellules CEE ont été combinées avec un deuxième groupe de 239 cellules CEE identifiées dans l’ensemble de données régionales (Fig. 2a, à droite) pour un total de 549 cellules CEE. Un groupe de 16 cellules co-exprimait les marqueurs EEC Chga et Chgb avec les marqueurs des cellules de Paneth, notamment Lyz1, Defa5 et Defa22, et ont donc été interprétées comme des doublets et retirées de l’analyse, laissant 533 cellules EEC, qui ont servi de base à l’analyse de la figure 3. Pour comparer les profils d’expression des entérocytes de l’intestin grêle proximal et distal (Fig. 2b), les 1 041 entérocytes identifiés à partir de 11 665 cellules dans l’ensemble de données régionales (Fig. 2a) ont été utilisés.

Définir les signatures de type cellulaire

Pour identifier les gènes maximalement spécifiques des types cellulaires, nous avons effectué des tests d’expression différentielle entre chaque paire de clusters pour toutes les comparaisons par paires possibles. Ensuite, pour un cluster donné, les gènes de signature putatifs ont été filtrés en utilisant la valeur Q FDR maximale et classés par le log2(fold change) minimal. Le changement de pli minimum et la valeur Q maximum représentent la taille d’effet la plus faible dans toutes les comparaisons par paires ; il s’agit donc d’un critère rigoureux. Les gènes de signature de type cellulaire présentés dans la figure 1c, les données étendues de la figure 2b, les données étendues de la figure 8e et les tableaux supplémentaires 2-4 et 8 ont été obtenus en utilisant un FDR maximum de 0,05 et un log2(fold change) minimum de 0,5. Dans le cas des signatures de types cellulaires post-mitotiques, tous les gènes ont passé ce seuil dans les ensembles de données 3′ (Fig. 1c) et complets (Extended Data Fig. 2b).

Dans le cas des gènes de signature pour les sous-types au sein des types cellulaires (Fig. 3b, Fig. 4b, Extended Data Fig. 7b), une valeur P combinée (sur l’ensemble des tests par paire) pour l’enrichissement a été calculée à l’aide de la méthode de Fisher – un critère plus indulgent que la simple prise en compte de la valeur P maximale – et une valeur Q FDR maximale de 0,01 a été utilisée, ainsi qu’un seuil de log2(fold change) minimal de 0,25 pour les sous-types de cellules tuft (Fig. 4b, Extended Data Fig. 7b, tableau supplémentaire 7) et de 0,1 pour les sous-types EEC (Fig. 3b, tableau supplémentaire 6). Tous les gènes de la signature des cellules touffues ont passé ce seuil dans les ensembles de données 3′ (Fig. 4b) et complets (Données étendues Fig. 7b), tandis que les signatures des sous-types EEC ont été définies en utilisant 3′ seulement. En raison du faible nombre de cellules (n = 18), la valeur P combinée de Fisher a également été utilisée pour la signature cellulaire microfold in vivo, avec un seuil FDR de 0,001 (Fig. 5d, Tableau supplémentaire 8). Les gènes marqueurs ont été classés par log2(fold change) minimum. Les tests d’expression différentielle ont été effectués à l’aide du test U de Mann-Whitney (également connu sous le nom de test de Wilcoxon rank-sum) mis en œuvre à l’aide de la fonction R wilcox.test. Pour les expériences d’infection (Fig. 6), nous avons utilisé un modèle de « hurdle » en deux parties pour contrôler à la fois la qualité technique et la variation entre les souris. Ce modèle a été mis en œuvre à l’aide du paquet R MAST59, et les valeurs P pour l’expression différentielle ont été calculées à l’aide du test du rapport de vraisemblance. La correction des tests d’hypothèses multiples a été effectuée en contrôlant le FDR60 à l’aide de la fonction R p.adjust.

Cotation des cellules à l’aide d’ensembles de gènes signatures

Pour obtenir un score pour un ensemble spécifique de n gènes dans une cellule donnée, un ensemble de gènes  » de fond  » a été défini pour contrôler les différences de couverture de séquençage et de complexité de la bibliothèque entre les cellules d’une manière similaire à la réf. 12. L’ensemble de gènes de fond a été sélectionné pour être similaire aux gènes d’intérêt en termes de niveau d’expression. Plus précisément, les 10n voisins les plus proches dans l’espace bidimensionnel défini par l’expression moyenne et la fréquence de détection dans toutes les cellules ont été sélectionnés. Le score de signature pour cette cellule a ensuite été défini comme l’expression moyenne des n gènes de signature dans cette cellule, moins l’expression moyenne des 10n gènes de fond dans cette cellule.

Estimations des fréquences d’échantillonnage des types de cellules

Pour chaque type de cellule, la probabilité d’observer au moins n cellules dans un échantillon de taille k est modélisée à l’aide de la fonction de distribution cumulative d’une binomiale négative NBcdf(k, n, p), où p est l’abondance relative de ce type de cellule. Pour m types de cellules avec le même paramètre p, la probabilité globale de voir chaque type au moins n fois est NBcdf(k ; n, p)m. Une telle analyse peut être effectuée avec des paramètres spécifiés par l’utilisateur à http://satijalab.org/howmanycells.

Dendrogramme EEC

Les vecteurs d’expression moyens ont été calculés pour l’ensemble des 12 clusters de sous-ensembles EEC, en utilisant les valeurs log2(TPM + 1), et restreints au sous-ensemble de 1 361 gènes identifiés comme significativement variables entre les susbsets EEC (P < 0,05), comme décrit ci-dessus. Les vecteurs d’expression moyens incluant ces gènes ont été regroupés hiérarchiquement à l’aide du package R pvclust (distance de Spearman, méthode de regroupement ward.D2), qui fournit des estimations de confiance bootstrap sur chaque nœud du dendrogramme sous forme de valeur P empirique sur 100 000 essais (Extended Data Fig. 6a).

Facteurs de transcription spécifiques du type cellulaire, RCPG et protéines répétées riches en leucine

Une liste de tous les gènes identifiés comme agissant en tant que facteurs de transcription chez la souris a été obtenue auprès d’AnimalTFDB61. L’ensemble des RCPG a été obtenu à partir de la base de données UniProt (http://www.uniprot.org/uniprot/?query=family%3A%22g+protein+coupled+receptor%22+AND+organism%3A%22Mouse+%5B10090%5D%22+AND+reviewed%3Ayes&sort=score). Les annotations fonctionnelles de chaque protéine (Extended Data Fig. 2d) ont été obtenues auprès de la British Pharmacological Society (BPS) et de l’International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). La liste des protéines à répétition riche en leucine a été tirée de la réf. 62. Pour établir la correspondance entre les noms de gènes humains et de gènes de souris, les orthologues humains et de souris ont été téléchargés à partir d’Ensembl (dernière version 86 ; http://www.ensembl.org/biomart/martview), et les synonymes des gènes humains et de souris à partir du NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). Pour chaque gène humain à répétition riche en leucine, tous les synonymes humains ont été mis en correspondance avec le gène orthologue chez la souris à l’aide de la liste des orthologues, et les noms des gènes de la souris ont été mis en correspondance avec ceux des données unicellulaires à l’aide de la liste des synonymes.

Les facteurs de transcription, les RCPG et les protéines à répétition riche en leucine enrichis par type cellulaire ont ensuite été identifiés en croisant la liste des gènes enrichis dans chaque type cellulaire avec les listes des facteurs de transcription, des RCPG et des protéines à répétition riche en leucine définies ci-dessus. Les gènes enrichis par type de cellule ont été définis à l’aide de l’ensemble de données SMART-Seq2 comme étant ceux dont le log2(fold change) minimum était de 0 et le FDR maximum de 0,5, en retenant un maximum de 10 gènes par type de cellule dans les figures 2e, f des données étendues (les listes complètes sont fournies dans le tableau supplémentaire 5). En outre, un panel plus étendu de RCPG spécifiques du type de cellule a été identifié (Données étendues Fig. 2d) en sélectionnant un seuil plus indulgent. Cela a été réalisé en comparant chaque type de cellule à toutes les autres cellules, au lieu des comparaisons par paires décrites dans la section précédente, et en sélectionnant tous les gènes de RCPG qui étaient exprimés de manière différentielle (FDR < 0,001).

Tester les changements dans les proportions de types de cellules

Nous modélisons le nombre détecté de chaque type de cellule dans chaque souris analysée comme une variable de comptage aléatoire en utilisant un processus de Poisson. Le taux de détection est ensuite modélisé en fournissant le nombre total de cellules profilées dans une souris donnée comme variable de décalage, avec la condition de chaque souris (traitement ou contrôle) fournie comme covariable. Le modèle a été ajusté à l’aide de la commande R glm du paquet stats. La valeur P de l’importance de l’effet produit par le traitement a été évaluée à l’aide d’un test de Wald sur le coefficient de régression.

Pour l’évaluation de l’importance des distributions spatiales des sous-ensembles d’EEC (figure 3e), la comparaison a impliqué plus de deux groupes. En particulier, notre hypothèse nulle était que la proportion de chaque sous-ensemble EEC détecté dans les trois régions intestinales (duodénum, jéjunum et iléon) était égale. Pour tester cette hypothèse, nous avons utilisé l’analyse de la variance (ANOVA) avec un test du χ2 sur l’ajustement du modèle de Poisson décrit ci-dessus, mis en œuvre à l’aide de la fonction anova du paquet stats.

Enrichissement des ensembles de gènes et analyse de l’ontologie des gènes

L’analyse de l’ontologie des gènes a été réalisée à l’aide du paquet R goseq63, en utilisant les gènes significativement exprimés de manière différentielle (FDR < 0.05) comme gènes cibles, et tous les gènes exprimés avec log2(TPM + 1) > 3 dans au moins dix cellules comme arrière-plan.

Data availability

Disponibilité du code