Un estudio unicelular del epitelio del intestino delgado

Ratones

Todo el trabajo con ratones se llevó a cabo de acuerdo con los Comités Institucionales de Cuidado y Uso de Animales (IACUC) y con las directrices pertinentes en el Instituto Broad y el Instituto Tecnológico de Massachusetts, con los protocolos 0055-05-15 y 0612-058-18, respectivamente. Para todos los experimentos, los ratones se asignaron aleatoriamente a los grupos de tratamiento después de emparejar por el sexo y la edad a ratones C57BL/6J hembra o macho de 7-10 semanas de edad de tipo salvaje o Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP), obtenidos del Jackson Laboratory (Bar Harbour) o ratones Gfi1beGFP/+ (Gfi1b-GFP)43. Los ratones se alojaron en condiciones libres de patógenos específicos en las instalaciones de animales del Instituto Broad, el Instituto Tecnológico de Massachusetts o la Escuela de Salud Pública T. H. Chan de Harvard.

Infección por Salmonella enterica y H. polygyrus. Los ratones C57BL/6J (Jackson Laboratory) fueron infectados con 200 larvas de tercer estadio de H. polygyrus o 108 de Salmonella enterica, mantenidos en condiciones libres de patógenos específicos en el Massachusetts General Hospital (Charlestown), con el protocolo 2003N000158. H. polygyrus se propagó como se ha descrito previamente44. Los ratones fueron eutanasiados 3 y 10 días después de la infección por H. polygyrus. Para Salmonella enterica, los ratones se infectaron con una cepa de S. Typhimurium naturalmente resistente a la estreptomicina (108 células) como se ha descrito previamente44, y se les aplicó la eutanasia 48 h después de la infección.

Disociación de células y aislamiento de criptas

Aislamiento de criptas. Se aisló el intestino delgado de ratones C57BL/6J de tipo salvaje, Lgr5-GFP o Gfi1b-GFP y se enjuagó en PBS frío. El tejido se abrió longitudinalmente y se cortó en pequeños fragmentos de aproximadamente 2 mm de longitud. El tejido se incubó en EDTA-PBS 20 mM en hielo durante 90 minutos, agitando cada 30 minutos. A continuación se agitó enérgicamente el tejido y se recogió el sobrenadante como fracción 1 en un nuevo tubo cónico. El tejido se incubó en EDTA-PBS fresco y se recogió una nueva fracción cada 30 minutos. Se recogieron fracciones hasta que el sobrenadante estaba formado casi exclusivamente por criptas. La fracción final (enriquecida en criptas) se lavó dos veces en PBS, se centrifugó a 300g durante 3 minutos y se disoció con TrypLE express (Invitrogen) durante 1 minuto a 37 °C. A continuación, la suspensión unicelular se pasó por un filtro de 40 μm y se tiñó para FACS para scRNA-seq (más abajo) o se utilizó para el cultivo de organoides. Confirmamos la solidez de este método probando otros métodos de aislamiento de células individuales -ya sea «entero» (raspando el revestimiento epitelial) o «enriquecido con vellosidades» (fracción 1; véase más arriba)- y descubrimos que, debido a la alta tasa de mortalidad (a través de la anoikis) de las células diferenciadas posmitóticas (cuyo componente principal son los enterocitos maduros), la suspensión de células individuales enriquecida con criptas representa fielmente la composición de los tipos de células del intestino delgado (datos no mostrados).

Aislamiento de epitelios asociados al folículo. Las células epiteliales de los epitelios asociados a los folículos se aislaron extrayendo pequeñas secciones (0,2-0,5 cm) que contenían manchas de Peyer del intestino delgado de ratones C57Bl/6J o Gfi1beGFP/+.

Selección de células

Para los experimentos de scRNA-seq de longitud completa basados en placas, se utilizó una máquina FACS (Astrios) para clasificar una sola célula en cada pocillo de una placa PCR de 96 pocillos que contenía 5 μl de tampón TCL con un 1% de 2-mercaptoetanol. Para el aislamiento de EpCAM+, las células se tiñeron para 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience); para las células epiteliales específicas, también se tiñeron para CD24+/- (eBioscience) y c-Kit+/- (eBioscience). Para enriquecer las poblaciones específicas de células epiteliales intestinales, se aislaron células de ratones Lgr5-GFP, se tiñeron con los anticuerpos mencionados anteriormente y se marcaron con GFP-high (células madre), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (linajes secretores) o GFP-/CD24-/EpCAM+ (células epiteliales). Para una mejor recuperación de las células Paneth, permitimos parámetros de dispersión lateral y frontal más altos en combinación con CD24+/c-Kit+ para verificar la recuperación de las células Paneth en las células EpCAM+. Para el aislamiento de tuft-2, las células epiteliales de tres ratones diferentes se tiñeron como en el caso anterior, pero utilizando EpCAM+/CD45+ para clasificar 2.000 células individuales. Utilizamos una puerta de clasificación indulgente para asegurarnos de obtener un número suficiente de estas raras células tuft-2, lo que condujo a una mayor tasa de contaminación de células T, que eliminamos en nuestro análisis de células individuales utilizando la agrupación no supervisada.

Para la clasificación de scRNA-seq de longitud completa, la placa de 96 pocillos se selló herméticamente con un Microseal F y se centrifugó a 800g durante 1 min. La placa se congeló inmediatamente en hielo seco y se mantuvo a -80 °C hasta que estuvo lista para la limpieza del lisado. Las células de la población a granel se clasificaron en un tubo Eppendorf que contenía 100 μl de solución de TCL con 2-mercaptoetanol al 1% y se almacenó a -80 °C.

Para el scRNA-seq basado en gotas, las células se clasificaron con los mismos parámetros que para el scRNA-seq basado en placas, pero se clasificaron en un tubo Eppendorf que contenía 50 μl de BSA-PBS al 0,4% y se almacenaron en hielo hasta que se procedió a la plataforma de células individuales GemCode.

ScRNA-seq basado en placas

Células individuales. Las bibliotecas se prepararon utilizando un protocolo modificado de SMART-Seq216. En resumen, la limpieza del lisado de ARN se llevó a cabo utilizando perlas RNAClean XP (Agencourt), seguido de la transcripción inversa con Maxima Reverse Transcriptase (Life Technologies) y la amplificación de la transcripción completa (WTA) con KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) durante 21 ciclos. Los productos de la WTA se purificaron con perlas Ampure XP (Beckman Coulter), se cuantificaron con Qubit dsDNA HS Assay Kit (ThermoFisher) y se evaluaron con un chip de ADN de alta sensibilidad (Agilent). Las bibliotecas de ARN-seq se construyeron a partir de los productos de WTA purificados utilizando el Nextera XT DNA Library Preperation Kit (Illumina). En cada placa, la población y los controles sin células se procesaron utilizando el mismo método que para las células individuales. Las bibliotecas se secuenciaron en un Illumina NextSeq 500.

Muestras a granel. Las muestras de población a granel se procesaron extrayendo el ARN con el RNeasy Plus Micro Kit (Qiagen) según las recomendaciones del fabricante y, a continuación, se procedió a realizar el protocolo SMART-Seq2 modificado tras la limpieza del lisado, tal y como se ha descrito anteriormente.

ScRNA-seq basado en gotas

Las células individuales se procesaron a través de la plataforma GemCode Single Cell utilizando los kits GemCode Gel Bead, Chip y Library (10X Genomics, Pleasanton) según el protocolo del fabricante. En resumen, las células individuales se clasificaron en BSA-PBS al 0,4%. Se añadieron 6.000 células a cada canal con una tasa de recuperación media de 1.500 células. A continuación, las células se dividieron en perlas de gel en emulsión en el instrumento GemCode, donde se produjo la lisis celular y la transcripción inversa del ARN con código de barras, seguida de la amplificación, el cizallamiento y la fijación del adaptador de 5′ y el índice de muestra. Las bibliotecas se secuenciaron en un Illumina NextSeq 500.

Inmunofluorescencia y smFISH

Inmunofluorescencia. La tinción de los tejidos del intestino delgado se realizó como se ha descrito previamente34. En resumen, los tejidos se fijaron durante 14 h en formol, se incrustaron en parafina y se cortaron en secciones de 5 μm de grosor. Las secciones se desparafinaron utilizando técnicas estándar, se incubaron con anticuerpos primarios durante la noche a 4 °C y luego con anticuerpos secundarios a temperatura ambiente durante 30 minutos. Los portaobjetos se montaron con Slowfade Mountant + DAPI (Life Technologies, S36964) y se sellaron.

smFISH. Se utilizó un RNAScope Multiplex Flourescent Kit (Advanced Cell Diagnostics) según las recomendaciones del fabricante con las siguientes alteraciones. El tiempo de ebullición de la recuperación de la diana se ajustó a 12 minutos y la incubación con la proteasa IV a 40 °C se ajustó a 8 minutos. Los portaobjetos se montaron con Slowfade Mountant+DAPI (Life Technologies, S36964) y se sellaron.

Inmunofluorescencia combinada y smFISH. Esto se implementó realizando primero smFISH como se ha descrito anteriormente, con los siguientes cambios. Después de la Amp 4, las secciones de tejido se lavaron en tampón de lavado, se incubaron con anticuerpos primarios durante la noche a 4 °C, se lavaron en 1× TBST tres veces y luego se incubaron con anticuerpos secundarios durante 30 min a temperatura ambiente. Los portaobjetos se montaron con Slowfade Mountant + DAPI (Life Technologies, S36964) y se sellaron.

Análisis de imágenes

Las imágenes de las secciones de tejido se tomaron con un microscopio confocal Fluorview FV1200 utilizando Kalman y emisión láser secuencial para reducir el ruido y la superposición de señales. Se añadieron barras de escala a cada imagen utilizando el software confocal FV10-ASW 3.1 Viewer. Las imágenes se superpusieron y visualizaron utilizando el software Image J45.

Anticuerpos y sondas

Cultivos de organoides intestinales

Después del aislamiento de las criptas, la suspensión unicelular se resuspendió en Matrigel (BD Bioscience) con 1 μM de péptido Jagged-1 (Ana-Spec). Se sembraron aproximadamente 300 criptas incrustadas en 25 μl de Matrigel en cada pocillo de una placa de 24 pocillos. Una vez solidificado, el Matrigel se incubó en 600 μl de medio de cultivo (DMEM/F12 avanzado, Invitrogen) con estreptomicina/penicilina y glutamatax y se complementó con EGF (100 ng ml-1, Peprotech), R-espondina-1 (600 ng ml-1, R&D) noggin (100 ng ml-1, Prepotech), dihidrocloruro monohidrato de Y-276432 (10 μM, Tochris), N-acetil-1-cisteína (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) y Wnt3A (25 ng ml-1, R&D Systems). Se sustituyó el medio fresco el día 3, y los organoides se pasaron por disociación con TrypLE y se volvieron a suspender en nuevo Matrigel el día 6 con una proporción de división de 1:3. En determinados experimentos, los organoides se trataron adicionalmente con RANKL (100 ng ml-1, Biolegends). Los organoides tratados se disociaron y se sometieron a scRNA-seq utilizando ambos métodos.

PCR cuantitativa

El ADNc de 16 células individuales de tuft-1, tuft-2 y EpCam+ al azar de las placas de scRNA-seq basadas en la longitud completa se utilizó para la qPCR relativa. La expresión génica se analizó mediante PCR cuantitativa en tiempo real en un LightCycler 480 Instrument II (Roche) utilizando la mezcla LightCycler 480 SYBR green (Roche) con los siguientes conjuntos de cebadores: HPRT1-F, GTTAAGCAGTACAGCCCCAAA; HPRT1-R, AGGGCATATCCAACAAACTT; UBC-F, CAGCCGTATCTTCCCAGACT; UBC-R, CTCAGAGGATGCCAGTAATCTA; tslp-F, TACTCTCAATCCTCCGCTG; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC; Alpi-F, TCCTACACCTCCATTCTATGG, Alpi-R, CCGCCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCTACACCATC; Dclk1-R, CCAGCTTCTTAAAGGCTCGAT. Los cebadores de qPCR se diseñaron para un límite de exón-exón en todos los transcritos.

Análisis computacional

Preprocesamiento de los datos de scRNA-seq basados en gotas. El desmultiplexado, la alineación con el transcriptoma de mm10 y el colapso del identificador molecular único (UMI) se realizaron utilizando el kit de herramientas Cellranger (versión 1.0.1) proporcionado por 10X Genomics. Para cada célula, cuantificamos el número de genes para los que se mapeó al menos una lectura, y luego excluimos todas las células con menos de 800 genes detectados. Los valores de expresión Ei,j para el gen i en la célula j se calcularon dividiendo los recuentos de UMI para el gen i por la suma de los recuentos de UMI en la célula j, para normalizar las diferencias de cobertura, y luego multiplicando por 10.000 para crear valores similares a TPM, y finalmente calculando log2(TPM + 1). La corrección por lotes se llevó a cabo utilizando ComBat46 tal y como se implementa en el paquete R sva47, utilizando el modo de ajuste paramétrico por defecto. El resultado fue una matriz de expresión corregida, que se utilizó como entrada para el análisis posterior.

La selección de genes variables se realizó ajustando un modelo lineal generalizado a la relación entre el coeficiente de variación al cuadrado y el nivel de expresión medio en el espacio logarítmico, y seleccionando los genes que se desviaban significativamente (P < 0,05) de la curva ajustada48.

Preprocesamiento de los datos de scRNA-seq de SMART-Seq2. Los archivos BAM se convirtieron en FASTQs fusionados y desmultiplexados utilizando el paquete de software Bcl2Fastq v2.17.1.14 proporcionado por Illumina. Las lecturas de extremo emparejado se asignaron al transcriptoma del ratón UCSC mm10 utilizando Bowtie49 con los parámetros ‘-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6’, que permite la alineación de secuencias con un desajuste. Los niveles de expresión de los genes se cuantificaron utilizando los valores de TPM calculados por RSEM50 v1.2.3 en modo paired-end. Para cada célula, cuantificamos el número de genes para los que se mapeó al menos una lectura, y luego excluimos todas las células con menos de 3.000 genes detectados o con un mapeo del transcriptoma inferior al 40%. A continuación, identificamos los genes altamente variables como se ha descrito anteriormente.

Reducción de la dimensionalidad mediante PCA y t-SNE. Restringimos la matriz de expresión a los subconjuntos de genes variables y células de alta calidad señalados anteriormente, y luego centramos y escalamos los valores antes de introducirlos en el análisis de componentes principales (PCA), que se implementó utilizando la función R prcomp del paquete stats para el conjunto de datos SMART-seq2. Para el conjunto de datos basado en gotas, utilizamos una aproximación aleatoria al PCA, implementada mediante la función rpca del paquete R rsvd, con el parámetro k fijado en 100. Esta aproximación de bajo rango se utilizó porque es varios órdenes de magnitud más rápida de calcular para matrices muy amplias. Dado que muchos componentes principales explican muy poco de la varianza, la relación señal-ruido puede mejorarse sustancialmente seleccionando un subconjunto de n componentes principales «significativos». Tras el ACP, se identificaron los componentes principales significativos mediante la prueba de permutación51 , implementada con la función permutationPA del paquete R jackstraw. Esta prueba identificó 13 y 15 componentes principales significativos en los conjuntos de datos 10X y SMART-Seq2 de la Fig. 1b y Extended Data Fig. 2a, respectivamente. Las puntuaciones de sólo estos componentes principales significativos se utilizaron como entrada para el análisis posterior.

Para la visualización, la dimensionalidad de los conjuntos de datos se redujo aún más utilizando la versión aproximada ‘Barnes-hut’ de t-SNE52,53. Esto se implementó utilizando la función Rtsne del paquete Rtsne R utilizando 20.000 iteraciones y un ajuste de perplejidad que varió de 10 a 30 dependiendo del tamaño del conjunto de datos.

Identificación de trayectorias de diferenciación celular utilizando mapas de difusión

Antes de ejecutar la reducción de la dimensionalidad del mapa de difusión seleccionamos genes altamente variables en los datos de la siguiente manera. En primer lugar, ajustamos un modelo nulo para la variabilidad de la expresión génica célula-célula en los datos, utilizando una relación de ley de potencia entre el coeficiente de variación y la media de los recuentos de UMI de todos los genes expresados, de forma similar a trabajos anteriores54. A continuación, calculamos para cada gen la diferencia entre el valor de su coeficiente de variación observado y el esperado por el modelo nulo (CVdiff). El histograma de CVdiff mostraba una cola «gorda». Calculamos la media μ y la desviación estándar σ de esta distribución, y seleccionamos todos los genes para los que CVdiff > μ + 1,67σ, obteniendo 761 genes para su posterior análisis.

Realizamos la reducción de la dimensionalidad utilizando el enfoque del mapa de difusión22. En resumen, se calculó una matriz de transición célula-célula utilizando un núcleo gaussiano, con el ancho del núcleo ajustado a la vecindad local de cada célula55. Esta matriz se convirtió en una matriz markoviana después de la normalización. Los vectores propios de la derecha vi (i = 0, 1, 2, …) de esta matriz se calcularon y clasificaron en orden de valor propio decreciente λi (i = 0, 1, 2, …), tras excluir el vector propio «superior» v0, correspondiente a λ0 = 1 (que refleja la restricción de normalización de la matriz markoviana). Los restantes vectores propios vi (i = 1, 2, …) definen la incrustación del mapa de difusión y se denominan componentes de difusión (DCk, k = 1, 2, …). Observamos una brecha espectral entre λ4 y λ5, y por lo tanto retuvimos DC1-DC4 tanto para el conjunto de datos inicial (Extended Data Fig. 4) como para los datos extraídos de distintas regiones intestinales (Fig. 2c).

Eliminación de células inmunes contaminantes y dobletes

Aunque las células fueron clasificadas antes de la secuenciación utilizando EpCAM, se observó un pequeño número de células inmunes contaminantes en el conjunto de datos 10X. Estas 264 células fueron eliminadas mediante una ronda inicial de clustering no supervisado (clustering basado en la densidad del mapa t-SNE utilizando dbscan56 del paquete R fpc) porque formaban un cluster extremadamente distinto. Para el conjunto de datos SMART-Seq2, varias celdas eran valores atípicos en términos de complejidad de la biblioteca, que podrían corresponder a más de una célula individual por biblioteca de secuenciación («dobletes»). Estas células se eliminaron calculando el cuantil superior del 1% de la distribución de genes detectados por célula y eliminando cualquier célula en este cuantil.

Análisis de clústeres

Para agrupar las células individuales por su expresión, utilizamos un enfoque de agrupación no supervisado, basado en el algoritmo de agrupación de gráficos Infomap9, siguiendo los enfoques para los datos CyTOF de células individuales57 y scRNA-seq10. En resumen, construimos un gráfico de k-próximos en los datos utilizando, para cada par de células, la distancia euclidiana entre las puntuaciones de los componentes principales significativos para identificar a los k vecinos más cercanos. El parámetro k se eligió para que fuera coherente con el tamaño del conjunto de datos. En concreto, k se fijó en 200 y 80 para el conjunto de datos basado en gotas de 7.216 células (Fig. 1b) y para el conjunto de datos SMART-Seq2 de 1.522 células (Extended Data Fig. 2a), respectivamente. Los organoides tratados con RANKL contenían 5.434 células y k se fijó en 200; el conjunto de datos de Salmonella y H. polygyrus contenía 9.842 células y k se fijó en 500. Para los análisis de conglomerados dentro de los tipos de células, específicamente los subconjuntos de células enteroendocrinas y de penacho, utilizamos la distancia de correlación de Pearson en lugar de la distancia euclidiana, y establecimos k = 15, k = 30 y k = 40 para los subtipos enteroendocrinos (533 células), y para las 166 y 102 células de penacho en los conjuntos de datos 10X y SMART-Seq2, respectivamente. El gráfico del vecino más cercano se calculó utilizando la función nng del paquete R cccd. El gráfico k-cercano más cercano se utilizó entonces como entrada para Infomap9, implementado utilizando la función infomap.community del paquete R igraph.

Los clústeres detectados se asignaron a tipos de células o estados intermedios utilizando marcadores conocidos para los subtipos de células epiteliales intestinales. (Datos extendidos Fig. 1g, Datos extendidos Fig. 2a). Para el subanálisis de las células enteroendocrinas (EEC) (Fig. 3), cualquier grupo de clusters de progenitores EEC con correlaciones medias por pares entre las puntuaciones de componentes principales significativas de r > 0,85 se fusionó, dando lugar a cuatro clusters. Etiquetamos estos cuatro clusters como progenitor «A» sobre la base de altos niveles de Ghrl, o progenitor (temprano), (medio) o (tardío) (en ese orden) sobre la base de niveles decrecientes de genes del tallo (Slc12a2, Ascl2, Axin2) y del ciclo celular y niveles crecientes de factores reguladores de EEC conocidos (Neurod1, Neurod2 y Neurog3) (Datos Extendidos Fig. 5c, Tabla Suplementaria 6). Para el conjunto de datos SMART-Seq2, dos clusters que expresaban altos niveles de genes marcadores de células madre (Extended Data Fig. 2a) se fusionaron para formar un cluster «stem» y otros dos clusters se fusionaron para formar un cluster «TA».

Para el análisis de clústeres del conjunto de datos del epitelio asociado al folículo de 4.700 células, las células del micropliegue eran muy escasas (0,38%) y por ello se utilizó el método ClusterDP58 para identificarlas, ya que se comportó empíricamente mejor que el algoritmo del gráfico k-cercano en este conjunto de datos. Al igual que con los métodos de k-nearest-neighbour, ClusterDP se ejecutó utilizando puntuaciones de componentes principales significativas (P < 0,05) (19 en este caso) como entrada, y se implementó utilizando las funciones findClusters y densityClust del paquete R densityClust utilizando los parámetros rho = 1.1 y delta = 0,25.

Extracción de tipos celulares raros para su posterior análisis

La agrupación inicial del conjunto de datos del intestino completo (7.216 células; Fig. 1b) mostró un grupo de 310 células EEC y 166 células de penacho. Las células del penacho se tomaron «tal cual» para el subanálisis (Fig. 4a, b), mientras que las células EEC se combinaron con un segundo grupo de 239 células EEC que se identificaron en el conjunto de datos regionales (Fig. 2a, derecha) para un total de 549 células EEC. Un grupo de 16 células coexpresaban los marcadores EEC Chga y Chgb con los marcadores de las células Paneth, incluyendo Lyz1, Defa5 y Defa22, por lo que se interpretaron como dobletes y se eliminaron del análisis, dejando 533 células EEC, que fueron la base para el análisis de la Fig. 3. Para comparar los perfiles de expresión de los enterocitos del intestino delgado proximal y distal (Fig. 2b), se utilizaron los 1.041 enterocitos identificados a partir de 11.665 células en el conjunto de datos regionales (Fig. 2a).

Definición de las firmas del tipo de célula

Para identificar los genes máximamente específicos para los tipos de células, realizamos pruebas de expresión diferencial entre cada par de clústeres para todas las posibles comparaciones por pares. A continuación, para un clúster determinado, se filtraron los genes putativos de firma utilizando el valor Q de FDR máximo y se clasificaron por el log2 mínimo (cambio de pliegue). El cambio de pliegue mínimo y el valor Q máximo representan el tamaño del efecto más débil en todas las comparaciones por pares; por lo tanto, es un criterio estricto. Los genes de la firma del tipo de célula que se muestran en la Fig. 1c, en la Fig. 2b de Datos Ampliados, en la Fig. 8e de Datos Ampliados y en las Tablas Suplementarias 2-4 y 8 se obtuvieron utilizando un FDR máximo de 0,05 y un log2(fold change) mínimo de 0,5. En el caso de las firmas de tipos celulares postmitóticos, todos los genes superaron este umbral tanto en los conjuntos de datos de 3′ (Fig. 1c) como en los de longitud completa (Extended Data Fig. 2b).

En el caso de los genes de firma para subtipos dentro de los tipos celulares (Fig. 3b, Fig. 4b, Extended Data Fig. 7b), se calculó un valor P combinado (a través de las pruebas por pares) para el enriquecimiento utilizando el método de Fisher -un criterio más indulgente que simplemente tomar el valor P máximo- y se utilizó un valor FDR Q máximo de 0,01, junto con un corte de log2(fold change) mínimo de 0,25 para los subtipos de células del penacho (Fig. 4b, Extended Data Fig. 7b, Supplementary Table 7) y de 0,1 para los subtipos EEC (Fig. 3b, Supplementary Table 6). Todos los genes de la firma de las células del penacho superaron este corte tanto en los conjuntos de datos de 3′ (Fig. 4b) como en los de longitud completa (Extended Data Fig. 7b), mientras que las firmas de los subtipos EEC se definieron utilizando sólo 3′. Debido al bajo número de células (n = 18), también se utilizó el valor P combinado de Fisher para la firma de microcélulas in vivo, con un corte FDR de 0,001 (Fig. 5d, Tabla Suplementaria 8). Los genes marcadores se clasificaron por log2(fold change) mínimo. Las pruebas de expresión diferencial se llevaron a cabo utilizando la prueba U de Mann-Whitney (también conocida como la prueba de suma de rangos de Wilcoxon) implementada utilizando la función de R wilcox.test. Para los experimentos de infección (Fig. 6) utilizamos un modelo de «obstáculo» de dos partes para controlar tanto la calidad técnica como la variación entre ratones. Esto se implementó utilizando el paquete R MAST59, y los valores P para la expresión diferencial se calcularon utilizando la prueba de relación de verosimilitud. La corrección de las pruebas de hipótesis múltiples se realizó controlando el FDR60 utilizando la función de R p.adjust.

Puntuación de las células utilizando conjuntos de genes de firma

Para obtener una puntuación para un conjunto específico de n genes en una célula determinada, se definió un conjunto de genes de «fondo» para controlar las diferencias en la cobertura de secuenciación y la complejidad de la biblioteca entre las células de una manera similar a la ref. 12. El conjunto de genes de fondo se seleccionó para que fuera similar a los genes de interés en términos de nivel de expresión. En concreto, se seleccionaron los 10n vecinos más cercanos en el espacio bidimensional definido por la expresión media y la frecuencia de detección en todas las células. La puntuación de la firma para esa célula se definió entonces como la expresión media de los n genes de la firma en esa célula, menos la expresión media de los 10n genes de fondo en esa célula.

Estimaciones de las frecuencias de muestreo del tipo de célula

Para cada tipo de célula, la probabilidad de observar al menos n células en una muestra de tamaño k se modela utilizando la función de distribución acumulativa de una binomial negativa NBcdf(k, n, p), donde p es la abundancia relativa de este tipo de célula. Para m tipos de células con el mismo parámetro p, la probabilidad global de ver cada tipo al menos n veces es NBcdf(k; n, p)m. Este análisis puede realizarse con parámetros especificados por el usuario en http://satijalab.org/howmanycells.

Dendrograma de CEE

Se calcularon los vectores de expresión media para los 12 clústeres de subconjuntos de CEE, utilizando los valores log2(TPM + 1), y se restringieron al subconjunto de 1.361 genes identificados como significativamente variables entre los susbsets de CEE (P < 0,05), como se describió anteriormente. Los vectores de expresión promedio que incluyen estos genes se agruparon jerárquicamente utilizando el paquete R pvclust (distancia de Spearman, método de agrupación ward.D2), que proporciona estimaciones de confianza bootstrap en cada nodo del dendrograma como un valor P empírico sobre 100.000 ensayos (Extended Data Fig. 6a).

Factores de transcripción específicos del tipo de célula, GPCRs y proteínas de repetición ricas en leucina

De AnimalTFDB61 se obtuvo una lista de todos los genes identificados que actúan como factores de transcripción en ratones. El conjunto de GPCRs se obtuvo de la base de datos 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). Las anotaciones funcionales para cada proteína (Extended Data Fig. 2d) se obtuvieron de la Sociedad Farmacológica Británica (BPS) y de la Unión Internacional de Farmacología Básica y Clínica (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). La lista de proteínas de repetición ricas en leucina se tomó de la ref. 62. Para mapear los nombres de los genes humanos a los de los ratones, los ortólogos humanos y de los ratones se descargaron de Ensembl (última versión 86; http://www.ensembl.org/biomart/martview), y los sinónimos de los genes humanos y de los ratones del NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). Para cada gen humano de repetición rica en leucina, todos los sinónimos humanos se mapearon con el gen ortólogo en el ratón utilizando la lista de ortólogos, y los nombres de los genes de ratón se mapearon con los de los datos de una sola célula utilizando la lista de sinónimos.

Los factores de transcripción enriquecidos en el tipo de célula, los GPCRs y las proteínas de repetición ricas en leucina se identificaron entonces intersectando la lista de genes enriquecidos en cada tipo de célula con las listas de factores de transcripción, GPCRs y proteínas de repetición ricas en leucina definidas anteriormente. Los genes enriquecidos por tipo de célula se definieron utilizando el conjunto de datos SMART-Seq2 como aquellos con un log2(fold change) mínimo de 0 y un FDR máximo de 0,5, reteniendo un máximo de 10 genes por tipo de célula en Extended Data Fig. 2e, f (las listas completas se proporcionan en la Tabla Suplementaria 5). Además, se identificó un panel más extenso de GPCRs específicos del tipo celular (Extended Data Fig. 2d) seleccionando un umbral más indulgente. Esto se consiguió comparando cada tipo celular con todas las demás células, en lugar de las comparaciones por pares descritas en la sección anterior, y seleccionando todos los genes GPCR que se expresaban de forma diferencial (FDR < 0,001).

Prueba de los cambios en las proporciones de los tipos celulares

Modelamos el número detectado de cada tipo celular en cada ratón analizado como una variable de recuento aleatoria utilizando un proceso de Poisson. La tasa de detección se modela entonces proporcionando el número total de células perfiladas en un ratón dado como una variable de compensación, con la condición de cada ratón (tratamiento o control) proporcionada como una covariable. El modelo se ajustó utilizando el comando de R glm del paquete stats. El valor P para la significación del efecto producido por el tratamiento se evaluó mediante una prueba de Wald sobre el coeficiente de regresión.

Para la evaluación de la significación de las distribuciones espaciales de los subconjuntos de CEE (Fig. 3e), la comparación implicó más de dos grupos. En particular, nuestra hipótesis nula era que la proporción de cada subconjunto de CEE detectado en las tres regiones intestinales (duodeno, yeyuno e íleon) era igual. Para probar esta hipótesis, utilizamos el análisis de la varianza (ANOVA) con una prueba χ2 sobre el ajuste del modelo de Poisson descrito anteriormente, implementado mediante la función anova del paquete stats.

Enriquecimiento del conjunto de genes y análisis de ontología de genes

El análisis de ontología de genes se llevó a cabo mediante el paquete R goseq63, utilizando genes significativamente expresados de forma diferencial (FDR < 0.05) como genes objetivo, y todos los genes expresados con log2(TPM + 1) > 3 en al menos diez células como fondo.

Disponibilidad de datos

Disponibilidad de código