Um levantamento monocelular do epitélio do intestino delgado

Ratos

Todo o trabalho do mouse foi realizado de acordo com os Comitês Institucionais de Cuidados e Uso de Animais (IACUC) e com diretrizes relevantes no Broad Institute e no Massachusetts Institute of Technology, com os protocolos 0055-05-15 e 0612-058-18, respectivamente. Para todos os experimentos, os ratos foram designados aleatoriamente para grupos de tratamento após a combinação para o sexo e idade de 7-10 semanas de idade, do tipo selvagem feminino ou masculino C57BL/6J ou Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP), obtidos dos ratos Jackson Laboratory (Bar Harbour) ou Gfi1beGFP/+ (Gfi1b-GFP)43. Os ratos foram alojados sob condições específicas sem patógenos nas instalações do Broad Institute, Massachusetts Institute of Technology ou Harvard T. H. Chan School of Public Health.

Salmonella enterica e infecção por H. polygyrus. Ratos C57BL/6J (Laboratório Jackson) foram infectados com 200 larvas do terceiro estágio de H. polygyrus ou 108 Salmonella enterica, mantidas sob condições específicas sem patógenos no Massachusetts General Hospital (Charlestown), com protocolo 2003N000158. O H. polygyrus foi propagado como descrito anteriormente44. Os ratos foram eutanizados 3 e 10 dias após a infecção pelo H. polygyrus. Para Salmonella enterica, ratos foram infectados com uma estirpe SL1344 de S. Typhimurium (108 células) naturalmente resistente à estreptomicina, como descrito anteriormente44, e foram eutanizados 48 h após a infecção.

Separação de células e isolamento de cripta

Separação de cripta. O intestino delgado de ratos do tipo C57BL/6J, Lgr5-GFP ou Gfi1b-GFP foi isolado e enxaguado em PBS fria. O tecido foi aberto longitudinalmente e cortado em pequenos fragmentos de aproximadamente 2 mm de comprimento. O tecido foi incubado em 20 mM EDTA-PBS sobre gelo durante 90 min, agitando a cada 30 min. O tecido foi então sacudido vigorosamente e o sobrenadante foi coletado como fração 1 em um novo tubo cônico. O tecido foi incubado em EDTA-PBS fresco e uma nova fração foi coletada a cada 30 min. As frações foram coletadas até que o sobrenadante consistisse quase inteiramente de criptas. A fração final (enriquecida para criptas) foi lavada duas vezes em PBS, centrifugada a 300g por 3 min, e dissociada com TrypLE express (Invitrogen) por 1 min a 37 °C. A suspensão monocelular foi então passada por um filtro de 40-μm e corada para FACS para scRNA-seq (abaixo) ou utilizada para cultura de organóides. Confirmamos a robustez deste método testando métodos adicionais de isolamento de células únicas – “inteiras” (raspando o revestimento epitelial) ou “enriquecidas com vilosidades” (fração 1; ver acima) – e constatamos que, devido à alta taxa de mortalidade (via anoikis) de células diferenciadas pós-mitóticas (cujo componente primário é o enterócito maduro), a suspensão criptográfica de células únicas representa fielmente a composição dos tipos de células do intestino delgado (dados não mostrados).

Isolamento de epitélio associado ao folículo. As células epiteliais dos epitélios associados ao folículo foram isoladas pela extração de pequenos cortes (0,2-0,5 cm) contendo manchas de Peyer do intestino delgado de ratos C57Bl/6J ou Gfi1beGFP/+.

Selecção de células

Para experiências com scRNA-seq à base de placas, foi utilizada uma máquina FACS (Astrios) para separar uma única célula em cada poço de uma placa PCR de 96 poços contendo 5 μl de tampão TCL com 1% de 2-mercaptoetanol. Para o isolamento EpCAM+, as células foram coradas para 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience); para células epiteliais específicas, também foram coradas para CD24+/- (eBioscience) e c-Kit+/- (eBioscience). Para enriquecer para populações específicas de células epiteliais intestinais, as células foram isoladas de ratos Lgr5-GFP, coradas com os anticorpos acima mencionados e fechadas em GFP-high (células-tronco), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (linhagens secretoras) ou GFP-/CD24-/EpCAM+ (células epiteliais). Para uma melhor recuperação das células Paneth, permitimos parâmetros de dispersão lateral e de dispersão frontal mais elevados em combinação com o CD24+/c-Kit+ para verificar a recuperação das células Paneth nas células EpCAM+. Para isolamento do tufo 2, células epiteliais de três ratos diferentes foram coradas como acima, mas usando EpCAM+/CD45+ para classificar 2.000 células individuais. Usamos uma porta de classificação indulgente para garantir que obtivéssemos um número suficiente dessas raras células de tufo-2, o que levou a uma maior taxa de contaminação de células T, que removemos em nossa análise de células únicas usando agrupamento não supervisionado.

Para a classificação de scRNA-seq de comprimento total, a placa de 96 poços foi selada firmemente com um Microseal F e centrifugada a 800g durante 1 min. A placa foi imediatamente congelada em gelo seco e mantida a -80 °C até estar pronta para a limpeza do lisado. As células da população a granel foram classificadas em um tubo Eppendorf contendo 100 μl solução de TCL com 1% de 2-mercaptoetanol e armazenadas a -80 °C.

Para scRNA-seq à base de gotas, as células foram classificadas com os mesmos parâmetros que para scRNA-seq à base de placas, mas foram classificadas em um tubo Eppendorf contendo 50 μl de 0,4% de BSA-PBS e armazenadas em gelo até prosseguir para a plataforma monocelular GemCode.

ScRNA-seq à base de placas

Células únicas. As bibliotecas foram preparadas usando um protocolo SMART-Seq2 modificado16. Em resumo, a limpeza do lisado do RNA foi realizada usando contas RNAClean XP (Agencourt) seguidas por transcrição reversa com Maxima Reverse Transcriptase (Life Technologies) e amplificação de transcrição inteira (WTA) com KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) durante 21 ciclos. Os produtos WTA foram purificados com contas Ampure XP (Beckman Coulter), quantificados com o Qubit dsDNA HS Assay Kit (ThermoFisher) e avaliados com um chip de DNA de alta sensibilidade (Agilent). As bibliotecas de RNA-seq foram construídas a partir de produtos WTA purificados usando o Nextera XT DNA Library Preperation Kit (Illumina). Em cada placa, a população e os controles sem células foram processados usando o mesmo método que para as células individuais. As bibliotecas foram sequenciadas num Illumina NextSeq 500.

Amostras de Bulk. As amostras de população a granel foram processadas extraindo RNA com o RNeasy Plus Micro Kit (Qiagen) de acordo com as recomendações do fabricante, e então procedendo com o protocolo SMART-Seq2 modificado após a limpeza do lisado, como descrito acima.

Droplet-based scRNA-seq

Células únicas foram processadas através da plataforma de células únicas GemCode usando os GemCode Gel Bead, Chip and Library Kits (10X Genomics, Pleasanton) de acordo com o protocolo do fabricante. Em resumo, as células únicas foram classificadas em 0,4% BSA-PBS. 6.000 células foram adicionadas a cada canal com uma taxa média de recuperação de 1.500 células. As células foram então divididas em Gel Beads em Emulsão no instrumento GemCode, onde ocorreu lise celular e transcrição reversa de RNA em código de barras, seguida de amplificação, cisalhamento e adaptador 5′ e fixação de índice de amostra. As bibliotecas foram sequenciadas em uma Illumina NextSeq 500.

Immunofluorescência e smFISH

Immunofluorescência. A coloração dos tecidos do intestino delgado foi conduzida conforme descrito anteriormente34. Em resumo, os tecidos foram fixados por 14 h em formalina, embutidos em parafina e cortados em cortes de 5-μm de espessura. As secções foram desparafinizadas usando técnicas padrão, incubadas com anticorpos primários durante a noite a 4 °C e depois com anticorpos secundários à temperatura ambiente durante 30 min. As lâminas foram montadas com Slowfade Mountant + DAPI (Life Technologies, S36964) e seladas.

smFISH. Um Kit Flourescente Multiplex RNAScope (Advanced Cell Diagnostics) foi utilizado de acordo com as recomendações do fabricante com as seguintes alterações. O tempo de ebulição da recuperação do alvo foi ajustado para 12 min e a incubação com Protease IV a 40 °C foi ajustada para 8 min. As lâminas foram montadas com Slowfade Mountant+DAPI (Life Technologies, S36964) e seladas.

Imunofluorescência combinada e smFISH. Isto foi implementado pela primeira vez executando o smFISH como descrito acima, com as seguintes alterações. Após o Amp 4, as secções de tecido foram lavadas em tampão de lavagem, incubadas com anticorpos primários durante a noite a 4 °C, lavadas em 1× TBST três vezes e depois incubadas com anticorpos secundários durante 30 min à temperatura ambiente. As lâminas foram montadas com Slowfade Mountant + DAPI (Life Technologies, S36964) e seladas.

Image analysis

Images of tissue sections were taken with a confocal microscope Fluorview FV1200 using Kalman and sequential laser emission to reduce noise and signal overlap. Barras de escala foram adicionadas a cada imagem usando o software confocal FV10-ASW 3.1 Viewer. As imagens foram sobrepostas e visualizadas usando o software Image J45.

Antibodies and probes

Culturas de organóides intestinais

Solamento de cripta de seguimento, a suspensão monocelular foi ressuspendida em Matrigel (BD Bioscience) com 1 μM Jagged-1 peptide (Ana-Spec). Cerca de 300 criptas incrustadas em 25 μl de Matrigel foram semeadas em cada poço de uma placa de 24 poços. Uma vez solidificado, o Matrigel foi incubado em 600 μl em meio de cultura (Advanced DMEM/F12, Invitrogen) com estreptomicina/penicilina e glutamatax e suplementado com EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), Y-276432 di-hidrocloreto monohidratado (10 μM, Tochris), N-acetil-1-cisteína (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) e Wnt3A (25 ng ml-1, R&D Systems). Os meios frescos foram substituídos no dia 3, e os organóides foram passados por dissociação com TrypLE e ressuspendidos no novo Matrigel no dia 6 com uma razão de divisão 1:3. Para experiências seleccionadas, os organóides foram adicionalmente tratados com RANKL (100 ng ml-1, Biolegends). Os organóides tratados foram dissociados e submetidos a scRNA-seq usando ambos os métodos.

PCR quantitativa

cDNA de 16 células unitárias de tufo-1, tufo-2 e EpCam+ aleatórias, a partir das placas scRNA-seq baseadas no comprimento total, foram usadas para o qPCR relativo. A expressão gênica foi analisada por PCR quantitativa em tempo real em um LightCycler 480 Instrument II (Roche) usando a mistura verde do LightCycler 480 SYBR (Roche) com os seguintes conjuntos de primer: HPRT1-F, GTTAAGCAGTACAGCCCCAAA; HPRT1-R, AGGGGCATATCCAACAACAAACTT; UBC-F, CAGCCGTATATCTTCCCAGACT; UBC-R, CTCAGAGGGATGCCAGTAATCTA; tslp-F, TACTCTCAATCCTATCCCTGGCTGCTG; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC; Alpi-F, TCCTACACCTCCATTCTCTCTATGGG, Alpi-R, CCGCCTGCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCTACACCATC; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. Os primers qPCR foram projetados para um limite exon-exon em todas as transcrições.

Análise computacional

P>Pré-processamento de dados scRNA-seq baseados em gotas. Desmultiplexação, alinhamento com a transcriptoma mm10 e identificador molecular único (UMI)-collapsing foram realizados usando o kit de ferramentas Cellranger (versão 1.0.1) fornecido pela 10X Genomics. Para cada célula, quantificamos o número de genes para os quais pelo menos uma leitura foi mapeada, e depois excluímos todas as células com menos de 800 genes detectados. Os valores de expressão Ei,j para o gene i na célula j foram calculados dividindo as contagens UMI para o gene i pela soma das contagens UMI na célula j, para normalizar para diferenças de cobertura, e depois multiplicando por 10.000 para criar valores do tipo TPM, e finalmente calculando log2(TPM + 1). A correção em lote foi realizada usando ComBat46 como implementado no pacote R sva47, usando o modo de ajuste paramétrico padrão. A saída foi uma matriz de expressão corrigida, que foi usada como entrada para análise posterior.

Seleção de genes variáveis foi realizada ajustando um modelo linear generalizado à relação entre o coeficiente de variação quadrático e o nível médio de expressão no espaço logarítmico, e selecionando genes que se desviaram significativamente (P < 0.05) da curva ajustada48.

Pré-processamento de dados SMART-Seq2 scRNA-seq. Os arquivos BAM foram convertidos para FASTQs fundidos e desmultiplexados usando o pacote de software Bcl2Fastq v2.17.1.14 fornecido pelo Illumina-provided. As leituras do final do Paire foram mapeadas para o transcriptome do mouse UCSC mm10 usando Bowtie49 com os parâmetros ‘-q –phred33-quals -n 1 -e 9999999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6’, o que permite o alinhamento de seqüências com um mismatch. Os níveis de expressão dos genes foram quantificados usando valores TPM calculados pelo RSEM50 v1.2.3 em modo paired-end. Para cada célula, quantificamos o número de genes para os quais pelo menos uma leitura foi mapeada, e então excluímos todas as células com menos de 3.000 genes detectados ou um mapeamento transcriptométrico de menos de 40%. Identificamos então genes altamente variáveis como descrito acima.

Dimensionalidade reduzida usando PCA e t-SNE. Restringimos a matriz de expressão aos subconjuntos de genes variáveis e células de alta qualidade acima mencionados, e então centralizamos e escalamos valores antes de inseri-los na análise de componentes principais (PCA), que foi implementada usando a função R prcomp do pacote stats para o conjunto de dados SMART-seq2. Para o conjunto de dados baseado em gotas usamos uma aproximação randomizada para PCA, implementada usando a função rpca do pacote rsvd R, com o parâmetro k ajustado para 100. Esta aproximação de baixo nível foi usada porque é várias ordens de grandeza mais rápidas de calcular para matrizes muito largas. Dado que muitos componentes principais explicam muito pouco da variância, a relação sinal/ruído pode ser substancialmente melhorada selecionando um subconjunto de n componentes principais ‘significativos’. Após o PCA, componentes principais significativos foram identificados usando o teste de permutação51 , implementado usando a função permutaçãoPA do pacote jackstraw R. Este teste identificou 13 e 15 componentes principais significativos nos conjuntos de dados 10X e SMART-Seq2 da Fig. 1b e Dados Estendidos da Fig. 2a, respectivamente. As pontuações destes componentes principais significativos foram usadas como entrada para análise posterior.

Para visualização, a dimensionalidade dos conjuntos de dados foi ainda mais reduzida usando a versão aproximada ‘Barnes-hut’ do t-SNE52,53. Isto foi implementado usando a função Rtsne do pacote Rtsne R usando 20.000 iterações e uma configuração de perplexidade que variou de 10 a 30 dependendo do tamanho do conjunto de dados.

Identificando trajetórias de diferenciação celular usando mapas de difusão

Antes de executar a redução da dimensionalidade do mapa de difusão, selecionamos genes altamente variáveis nos dados da seguinte forma Primeiramente, ajustamos um modelo nulo de variabilidade de expressão de genes de células de base nos dados, usando uma relação de poder-direito entre o coeficiente de variação e a média da contagem do IMC de todos os genes expressos, semelhante ao trabalho anterior54. Em seguida, calculamos para cada gene a diferença entre o valor de seu coeficiente de variação observado e o esperado pelo modelo nulo (CVdiff). O histograma de CVdiff exibia uma cauda ‘gorda’. Calculamos a média μ e o desvio padrão σ desta distribuição, e selecionamos todos os genes para os quais CVdiff > μ + 1.67σ, produzindo 761 genes para análise posterior.

Realizamos a redução da dimensionalidade utilizando a abordagem de difusão-mapa22. Em resumo, uma matriz de transição célula-célula foi computada usando um kernel gaussiano, com a largura do kernel ajustada para a vizinhança local de cada célula55. Esta matriz foi convertida para uma matriz Markoviana após a normalização. Os autovectores vi (i = 0, 1, 2, …) direitos desta matriz foram computados e ordenados em ordem decrescente de valor próprio λi (i = 0, 1, 2, …), depois de excluir o autovector v0 ‘top’, correspondente a λ0 = 1 (que reflecte a restrição de normalização da matriz Markoviana). Os restantes autovectores vi (i = 1, 2, …) definem a incorporação do mapa de difusão e são referidos como componentes de difusão (DCk, k = 1, 2, …). Notamos uma lacuna espectral entre λ4 e λ5, e portanto reteve-se DC1-DC4 tanto para o conjunto de dados inicial (Extended Data Fig. 4) quanto para os dados extraídos de regiões intestinais distintas (Fig. 2c).

Remover células imunes contaminantes e doublets

Embora as células tenham sido classificadas antes do sequenciamento usando EpCAM, um pequeno número de células imunes contaminantes foi observado no conjunto de dados 10X. Estas 264 células foram removidas por uma ronda inicial de agrupamento não supervisionado (agrupamento baseado na densidade do mapa t-SNE usando dbscan56 do pacote R fpc) porque formaram um agrupamento extremamente distinto. Para o conjunto de dados SMART-Seq2, várias células foram aberrantes em termos de complexidade de biblioteca, o que poderia corresponder a mais de uma célula individual por biblioteca de sequenciamento (“doublets”). Estas células foram então removidas calculando o quantil topo de 1% da distribuição dos genes detectados por célula e removendo quaisquer células neste quantil.

Análise de agrupamento

Para agrupar células únicas pela sua expressão, utilizamos uma abordagem de agrupamento não supervisionada, baseada no algoritmo Infomap de agrupamento de gráficos9, seguindo abordagens para dados CyTOF de célula única57 e scRNA-seq10. Em resumo, construímos um gráfico de k-nearest-neighbour nos dados usando, para cada par de células, a distância euclidiana entre os escores de componentes principais significativos para identificar k vizinhos mais próximos. O parâmetro k foi escolhido para ser consistente com o tamanho do conjunto de dados. Especificamente, k foi definido como 200 e 80 para o conjunto de dados baseado em gotas de 7.216 células (Fig. 1b) e para o conjunto de dados SMART-Seq2 de 1.522 células (Extended Data Fig. 2a), respectivamente. Os organóides tratados com RANKL continham 5.434 células e k foi definido para 200; o conjunto de dados de Salmonella e H. polygyrus continha 9.842 células e k foi definido para 500. Para análises de agrupamento dentro dos tipos de células, especificamente os subconjuntos de células enteroendócrinas e tufo, utilizamos a distância de correlação de Pearson ao invés da distância Euclidiana, e estabelecemos k = 15, k = 30 e k = 40 para os subtipos enteroendócrinos (533 células), e para as 166 e 102 células de tufo nos conjuntos de dados 10X e SMART-Seq2, respectivamente. O gráfico de vizinhança mais próxima foi calculado usando a função nng do pacote R cccd. O gráfico k-nearest-neighbour foi então usado como entrada para o Infomap9, implementado usando a função infomap.community do pacote igraph R.

Clusters detectados foram mapeados para tipos de células ou estados intermediários usando marcadores conhecidos para subtipos de células epiteliais intestinais. (Dados Estendidos Fig. 1g, Dados Estendidos Fig. 2a). Para a sub-análise celular enteroendocrina (EEC) (Fig. 3), qualquer grupo de grupos de progenitores EEC com correlações parciais médias entre pontuações significativas dos componentes principais de r > 0,85 foi fundido, resultando em quatro grupos. Nós rotulamos esses quatro grupos de progenitores ‘A’ com base em altos níveis de Ghrl, ou progenitor (precoce), (médio) ou (tardio) (nessa ordem) com base em níveis decrescentes de genes do caule (Slc12a2, Ascl2, Axin2) e do ciclo celular e níveis crescentes de fatores regulatórios conhecidos da EEC (Neurod1, Neurod2 e Neurog3) (Dados Estendidos Fig. 5c, Tabela Suplementar 6). Para o conjunto de dados SMART-Seq2, dois clusters que expressam altos níveis de genes marcadores de células-tronco (Extended Data Fig. 2a) foram fundidos para formar um cluster ‘stem’ e dois outros clusters foram fundidos para formar um cluster ‘TA’.

Para a análise de agrupamento do conjunto de dados do epitélio associado ao folículo de 4.700 células, as células microfoldadas eram muito raras (0,38%) e assim o método ClusterDP58 foi usado para identificá-las porque teve um desempenho empírico melhor que o algoritmo do gráfico k-nearest-neighbour neste conjunto de dados. Assim como nos métodos k-nearest-neighbour, ClusterDP foi executado usando scores significativos (P < 0,05) dos componentes principais (19 neste caso) como input, e foi implementado usando as funções findClusters e densityClust do pacote densityClust R usando os parâmetros rho = 1.1 e delta = 0,25.

Extraindo tipos de células raras para análise posterior

A agregação inicial do conjunto de dados de todo o intestino (7.216 células; Fig. 1b) mostrou um cluster de 310 células EEC e 166 células de tufo. As células de tufo foram tomadas “como estão” para a subanálise (Fig. 4a, b), enquanto as células de EEC foram combinadas com um segundo agrupamento de 239 células de EEC que foram identificadas no conjunto de dados regionais (Fig. 2a, à direita) para um total de 549 células de EEC. Um grupo de 16 células co-expressas de marcadores de EEC Chga e Chgb com marcadores de células Paneth, incluindo Lyz1, Defa5 e Defa22, foram interpretadas como doublets e retiradas da análise, deixando 533 células de EEC, que foram a base para a análise na Fig. 3. Para comparar perfis de expressão de enterócitos do intestino delgado proximal e distal (Fig. 2b), foram utilizados os 1.041 enterócitos identificados a partir de 11.665 células do conjunto de dados regional (Fig. 2a).

Definindo assinaturas do tipo celular

Para identificar genes máximos específicos para os tipos celulares, fizemos testes de expressão diferencial entre cada par de agrupamentos para todas as comparações possíveis em pares. Então, para um dado cluster, os genes de assinaturas putativas foram filtrados usando o valor máximo FDR Q e classificados pelo log mínimo2(fold change). O mínimo de mudança de dobra e o máximo valor Q representam o tamanho do efeito mais fraco em todas as comparações de pares; é, portanto, um critério rigoroso. Os genes de assinatura tipo célula mostrados na Fig. 1c, Dados Estendidos Fig. 2b, Dados Estendidos Fig. 8e e Tabelas Suplementares 2-4 e 8 foram obtidos usando um FDR máximo de 0,05 e um log2(fold change) mínimo de 0,5. No caso de assinaturas pós-mitóticas do tipo célula, todos os genes passaram esse limite tanto no conjunto de dados 3′ (Fig. 1c) quanto no conjunto de dados full-length (Extended Data Fig. 2b).

No caso de genes de assinaturas para subtipos dentro de tipos de células (Fig. 3b, Fig. 4b, Extended Data Fig. 7b), um valor P combinado (nos testes de pares) para enriquecimento foi calculado usando o método de Fisher – um critério mais brando do que simplesmente tomar o valor P máximo – e um valor FDR Q máximo de 0,01 foi usado, juntamente com um corte de log2(fold change) mínimo de 0,25 para subtipos de células de tufo (Fig. 4b, Dados Estendidos Fig. 7b, Tabela Complementar 7) e de 0,1 para subtipos CEE (Fig. 3b, Tabela Complementar 6). Todos os genes na assinatura da célula de tufo passaram por este corte tanto no conjunto de dados 3′ (Fig. 4b) como no conjunto de dados completo (Extended Data Fig. 7b), enquanto que as assinaturas dos subtipos EEC foram definidas usando apenas 3′. Devido ao baixo número de células (n = 18), o valor P combinado de Fisher também foi utilizado para a assinatura in vivo de células em microfibras, com um corte FDR de 0,001 (Fig. 5d, Tabela Suplementar 8). Os genes marcadores foram classificados por log mínimo2(fold change). Os testes de expressão diferencial foram realizados usando o teste Mann-Whitney U (também conhecido como teste Wilcoxon rank-sum) implementado usando a função R wilcox.test. Para os experimentos de infecção (Fig. 6) usamos um modelo de ‘obstáculo’ em duas partes para controlar tanto a qualidade técnica quanto a variação de mouse para mouse. Isto foi implementado usando o pacote R MAST59, e os valores de P para expressão diferencial foram computados usando o teste de razão de verossimilhança. A correção de múltiplas hipóteses foi realizada através do controle do FDR60 usando a função R p.adjust.

Células escoreadoras usando conjuntos de genes de assinatura

Para obter um escore para um conjunto específico de n genes em uma determinada célula, foi definido um conjunto de genes ‘background’ para controlar as diferenças na cobertura de seqüenciamento e complexidade da biblioteca entre células de forma similar ao ref. 12. O conjunto de genes de fundo foi selecionado para ser similar aos genes de interesse em termos de nível de expressão. Especificamente, foram seleccionados os 10n vizinhos mais próximos no espaço bidimensional definido pela expressão média e frequência de detecção em todas as células. O escore de assinatura para aquela célula foi então definido como a expressão média dos n genes de assinatura naquela célula, menos a expressão média dos 10n genes de fundo naquela célula.

Etimações das frequências de amostragem do tipo de célula

Para cada tipo de célula a probabilidade de observar pelo menos n células numa amostra de tamanho k é modelada usando a função de distribuição cumulativa de um binômio negativo NBcdf(k, n, p), onde p é a abundância relativa deste tipo de célula. Para m tipos de células com o mesmo parâmetro p, a probabilidade geral de ver cada tipo pelo menos n vezes é NBcdf(k; n, p)m. Tal análise pode ser realizada com parâmetros especificados pelo usuário em http://satijalab.org/howmanycells.

EEC dendrogram

Vetores de expressão médios foram calculados para todos os 12 subconjuntos de EEC, usando valores log2(TPM + 1), e restritos ao subconjunto de 1.361 genes identificados como variável significativa entre os subconjuntos EEC (P < 0.05), como descrito acima. Os vetores de expressão médios incluindo esses genes foram hierarquicamente agrupados usando o pacote R pvclust (Spearman distance, ward.D2 clustering method), que fornece estimativas de confiança bootstrap em cada nó de dendrograma como um valor empírico de P acima de 100.000 ensaios (Dados Estendidos Fig. 6a).

Factores de transcrição específicos do tipo Célula, GPCRs e proteínas repetidas ricas em leucina

Uma lista de todos os genes identificados como agindo como factores de transcrição em ratos foi obtida de AnimalTFDB61. O conjunto de GPCRs foi obtido da base de dados 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). Anotações funcionais para cada proteína (Extended Data Fig. 2d) foram obtidas da British Pharmacological Society (BPS) e da International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). A lista de proteínas de repetição rica em leucina foi retirada da ref. 62. Para mapear os nomes dos genes humano e do rato, ortologues humanos e do rato foram descarregados do Ensembl (última versão 86; http://www.ensembl.org/biomart/martview), e sinónimos dos genes humano e do rato do NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). Para cada gene repetido humano rico em leucina, todos os sinônimos humanos foram mapeados para o gene ortológico em camundongos usando a lista de ortologues, e os nomes dos genes de camundongos foram mapeados para aqueles nos dados de célula única usando a lista de sinônimos.

Correntes de transcrição enriquecidos com leucina, GPCRs e proteínas repetidas ricas em leucina foram então identificadas através da intersecção da lista de genes enriquecidos em cada tipo de célula com as listas de fatores de transcrição, GPCRs e proteínas repetidas ricas em leucina definidas acima. Os genes enriquecidos por tipo de célula foram definidos usando o conjunto de dados SMART-Seq2 como aqueles com um log2(fold change) mínimo de 0 e um FDR máximo de 0,5, retendo um máximo de 10 genes por tipo de célula nos Dados Estendidos Fig. 2e, f (listas completas são fornecidas na Tabela Complementar 5). Além disso, foi identificado um painel mais extenso de GPCRs específicos do tipo de célula (Extended Data Fig. 2d) através da seleção de um limiar mais indulgente. Isto foi conseguido comparando cada tipo de célula com todas as outras células, ao invés das comparações em pares descritas na seção anterior, e selecionando todos os genes GPCR que foram diferentemente expressos (FDR < 0,001).

Teste para alterações nas proporções de tipo de célula

Modelizamos o número detectado de cada tipo de célula em cada mouse analisado como uma variável de contagem aleatória usando um processo de Poisson. A taxa de detecção é então modelada fornecendo o número total de células perfiladas em um determinado mouse como uma variável de offset, com a condição de cada mouse (tratamento ou controle) fornecida como uma covariada. O modelo foi ajustado usando o glm de comando R do pacote de estatísticas. O valor de P para a significância do efeito produzido pelo tratamento foi avaliado usando um teste de Wald sobre o coeficiente de regressão.

Para a avaliação da significância das distribuições espaciais dos subconjuntos EEC (Fig. 3e), a comparação envolveu mais de dois grupos. Em particular, nossa hipótese nula foi que a proporção de cada subconjunto de EEC detectado nas três regiões intestinais (duodeno, jejuno e íleo) foi igual. Para testar esta hipótese, utilizamos a análise de variância (ANOVA) com um teste χ2 no modelo de Poisson encaixado descrito acima, implementado usando a função anova do pacote stats.

Enriquecimento do conjunto genético e análise ontológica do gene

Análise ontológica do gene foi realizada usando o pacote goseq R63, usando genes significativamente expressos de forma diferente (FDR < 0.05) como genes alvo, e todos os genes expressos com log2(TPM + 1) > 3 em pelo menos dez células como fundo.

Disponibilidade de dados

Disponibilidade de código