En enkelsidig undersökning av tunntarmens epitel

Möss

Alt arbete med möss utfördes i enlighet med Institutional Animal Care and Use Committees (IACUC) och relevanta riktlinjer vid Broad Institute och Massachusetts Institute of Technology, med protokoll 0055-05-15 respektive 0612-058-18. För alla experiment tilldelades möss slumpmässigt till behandlingsgrupper efter matchning med avseende på kön och ålder av 7-10 veckor gamla hon- eller hanmöss av vildtyp C57BL/6J eller Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP), erhållna från Jackson Laboratory (Bar Harbour) eller Gfi1beGFP/+ (Gfi1b-GFP) möss43. Möss hölls under specifika patogenfria förhållanden vid djuranläggningar vid Broad Institute, Massachusetts Institute of Technology eller Harvard T. H. Chan School of Public Health.

Salmonella enterica och H. polygyrus-infektion. C57BL/6J-möss (Jackson Laboratory) infekterades med 200 larver i tredje stadiet av H. polygyrus eller 108 Salmonella enterica, som hölls under specifika patogenfria förhållanden vid Massachusetts General Hospital (Charlestown), med protokoll 2003N000158. H. polygyrus förökades enligt tidigare beskrivning44. Möss avlivades 3 och 10 dagar efter H. polygyrus-infektion. För Salmonella enterica infekterades möss med en naturligt streptomycinresistent SL1344-stam av S. Typhimurium (108 celler) enligt tidigare beskrivning44 och avlivades 48 timmar efter infektionen.

Celldissociation och kryptisolering

Kryptisolering. Tunntarmen från C57BL/6J-möss av vildtyp, Lgr5-GFP eller Gfi1b-GFP isolerades och sköljdes i kall PBS. Vävnaden öppnades i längsled och skars i små fragment som var ungefär 2 mm långa. Vävnaden inkuberades i 20 mM EDTA-PBS på is i 90 minuter och skakades var 30:e minut. Vävnaden skakades sedan kraftigt och supernatanten samlades upp som fraktion 1 i ett nytt koniskt rör. Vävnaden inkuberades i ny EDTA-PBS och en ny fraktion samlades upp var 30:e minut. Fraktioner samlades upp tills supernatanten nästan helt bestod av kryptor. Den sista fraktionen (berikad på kryptor) tvättades två gånger i PBS, centrifugerades vid 300 g i 3 minuter och dissocierades med TrypLE express (Invitrogen) i 1 minut vid 37 °C. Encellssuspensionen passerade sedan genom ett 40-μm-filter och färgades för FACS för scRNA-seq (nedan) eller användes för organoidkultur. Vi bekräftade denna metods robusthet genom att testa ytterligare metoder för isolering av enstaka celler – antingen ”hela” (skrapning av epitelfoderet) eller ”villusberikade” (fraktion 1; se ovan) – och fann att den kryptberikade encellssuspensionen, på grund av den höga dödligheten (via anoikis) av postmitotiska differentierade celler (vars primära komponent är mogna enterocyter), troget representerar sammansättningen av tunntarmcellernas celltyper (data inte visad).

Follikelassocierad epitelisolering. Epitelceller från de follikelassocierade epitelcellerna isolerades genom att extrahera små sektioner (0,2-0,5 cm) innehållande Peyer’s patches från tunntarmen från C57Bl/6J- eller Gfi1beGFP/+-möss.

Cellsortering

För plattbaserade full-längds scRNA-seq-experiment användes en FACS-maskin (Astrios) för att sortera en enskild cell till varje brunn i en 96-brunns PCR-platta som innehöll 5 μl TCL-buffert med 1 % 2-mercaptoetanol. För isolering av EpCAM+ färgades cellerna för 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience); för specifika epitelceller färgades vi också för CD24+/- (eBioscience) och c-Kit+/- (eBioscience). För att anrika specifika tarmepitelcellspopulationer isolerades celler från Lgr5-GFP-möss, färgades med de antikroppar som nämns ovan och gated på GFP-high (stamceller), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (sekretoriska linjer) eller GFP-/CD24-/EpCAM+ (epitelceller). För bättre återhämtning av Paneth-cellerna tillät vi högre parametrar för side scatter och forward scatter i kombination med CD24+/c-Kit+ för att verifiera återhämtning av Paneth-cellerna i EpCAM+-celler. För tuft-2-isolering färgades epitelceller från tre olika möss på samma sätt som ovan, men med EpCAM+/CD45+ för att sortera 2 000 enskilda celler. Vi använde en mild sorteringsgrind för att säkerställa att vi fick tillräckligt många av dessa sällsynta tuft-2-celler, vilket ledde till en högre kontamineringsgrad av T-celler, som vi tog bort i vår singelcellsanalys med hjälp av oövervakad klustring.

För sortering av scRNA-seq i full längd förseglades 96-brunnsplattan tätt med en Microseal F och centrifugerades vid 800 g i 1 minut. Plattan frystes omedelbart på torris och förvarades vid -80 °C tills den var klar för lysatrensning. Bulkpopulationens celler sorterades i ett Eppendorfrör som innehöll 100 μl lösning av TCL med 1 % 2-mercaptoetanol och förvarades vid -80 °C.

För droppbaserad scRNA-seq sorterades cellerna med samma parametrar som för plattbaserad scRNA-seq, men sorterades i ett Eppendorf-rör som innehöll 50 μl 0,4 % BSA-PBS och förvarades på is tills de gick vidare till GemCode-plattformen för enstaka celler.

Plattbaserad scRNA-seq

Enstaka celler. Bibliotek framställdes med hjälp av ett modifierat SMART-Seq2-protokoll16. I korthet utfördes RNA-lysatrening med RNAClean XP-beads (Agencourt) följt av omvänd transkription med Maxima Reverse Transcriptase (Life Technologies) och förstärkning av hela transkriptioner (WTA) med KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) i 21 cykler. WTA-produkterna renades med Ampure XP-beads (Beckman Coulter), kvantifierades med Qubit dsDNA HS Assay Kit (ThermoFisher) och bedömdes med ett högkänsligt DNA-chip (Agilent). RNA-seq-bibliotek konstruerades från renade WTA-produkter med hjälp av Nextera XT DNA Library Preperation Kit (Illumina). På varje platta bearbetades populationen och kontrollerna utan celler med samma metod som för de enskilda cellerna. Biblioteken sekvenserades på en Illumina NextSeq 500.

Plumpprover. Bulkpopulationsprover bearbetades genom att extrahera RNA med RNeasy Plus Micro Kit (Qiagen) enligt tillverkarens rekommendationer och sedan fortsätta med det modifierade SMART-Seq2-protokollet efter lysatrensning, enligt beskrivningen ovan.

Dropletbaserad scRNA-seq

Enskilda celler bearbetades genom GemCode Single Cell Platform med hjälp av GemCode Gel Bead, Chip and Library Kits (10X Genomics, Pleasanton) enligt tillverkarens protokoll. I korthet sorterades enskilda celler i 0,4 % BSA-PBS. 6 000 celler lades till i varje kanal med en genomsnittlig återvinningsgrad på 1 500 celler. Cellerna fördelades sedan i Gel Beads in Emulsion i GemCode-instrumentet, där celllys och streckkodad omvänd transkription av RNA ägde rum, följt av amplifiering, skjuvning och 5′-adapter- och provindexinfästning. Biblioteken sekvenserades på en Illumina NextSeq 500.

Immunofluorescens och smFISH

Immunofluorescens. Färgning av tunntarmsvävnader utfördes enligt tidigare beskrivning34. I korthet fixerades vävnaderna i 14 timmar i formalin, bäddades in i paraffin och skars i 5 μm tjocka sektioner. Sektionerna avparaffiniserades med standardteknik, inkuberades med primära antikroppar över natten vid 4 °C och sedan med sekundära antikroppar vid rumstemperatur i 30 minuter. Objektglasen monterades med Slowfade Mountant + DAPI (Life Technologies, S36964) och förseglades.

smFISH. Ett RNAScope Multiplex Flourescent Kit (Advanced Cell Diagnostics) användes enligt tillverkarens rekommendationer med följande ändringar. Kokationstiden för target retrieval justerades till 12 minuter och inkubationen med proteas IV vid 40 °C justerades till 8 minuter. Objektglasen monterades med Slowfade Mountant+DAPI (Life Technologies, S36964) och förseglades.

Kombinerad immunofluorescens och smFISH. Detta genomfördes genom att först utföra smFISH enligt beskrivningen ovan, med följande ändringar. Efter Amp 4 tvättades vävnadssektionerna i tvättbuffert, inkuberades med primära antikroppar över natten vid 4 °C, tvättades i 1× TBST tre gånger och inkuberades sedan med sekundära antikroppar i 30 minuter i rumstemperatur. Objektglasen monterades med Slowfade Mountant + DAPI (Life Technologies, S36964) och förseglades.

Bildanalys

Bilder av vävnadssektioner togs med ett konfokalmikroskop Fluorview FV1200 med Kalman och sekventiell laseremission för att minska brus och signalöverlappning. Skalstreck lades till varje bild med hjälp av den konfokala programvaran FV10-ASW 3.1 Viewer. Bilderna överlagrades och visualiserades med hjälp av programvaran Image J45.

Antikroppar och prober

Intestinal organoidkulturer

Efter isolering av kryptan resuspenderades singelcellssuspensionen i Matrigel (BD Bioscience) med 1 μM Jagged-1 peptid (Ana-Spec). Ungefär 300 kryptor inbäddade i 25 μl Matrigel såddes på varje brunn i en 24-brunnsplatta. När matrigelen hade stelnat inkuberades den i 600 μl kulturmedium (Advanced DMEM/F12, Invitrogen) med streptomycin/penicillin och glutamatax och kompletterades med EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), Y-276432 dihydrochloridmonohydrat (10 μM, Tochris), N-acetyl-1-cystein (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) och Wnt3A (25 ng ml-1, R&D Systems). Nya medier byttes ut dag 3, och organoiderna passagerades genom dissociation med TrypLE och återuppsattes i ny Matrigel dag 6 med ett delningsförhållande på 1:3. För utvalda experiment behandlades organoiderna dessutom med RANKL (100 ng ml-1, Biolegends). Behandlade organoider dissocierades och utsattes för scRNA-seq med båda metoderna.

Kvantitativ PCR

cDNA från 16 heltranskriptomamplifierade enskilda celler av tuft-1, tuft-2 och slumpmässiga EpCam+ från de fullängdsbaserade scRNA-seq-plattorna användes för relativ qPCR. Genuttrycket analyserades genom kvantitativ realtids-PCR på en LightCycler 480 Instrument II (Roche) med LightCycler 480 SYBR green mix (Roche) med följande primeruppsättningar: HPRT1-F, GTTAAGCAGCAGTACAGCCCCAAA; HPRT1-R, AGGGCATATCATCCAACACAACAAACTT; UBC-F, CAGCCGTATATATCTTCCCAGACT; UBC-R, CTCAGAGGGATGCCAGAGTAATCTA; tslp-F, TACTCTCTCAATCCTATCCCTCTGGCTG; Tlsp-R, CCATTTCCTCTGAGTACCGTCATTTC; Alpi-F, TCCTACACCACCTCCATTCTCTATGG, Alpi-R, CCGCCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCATCTACACCACCATC; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. qPCR-primers utformades för en exon-exon-gräns i alla transkript.

Beräkningsanalys

Förbehandling av droppbaserade scRNA-seq-data. De-multiplexering, anpassning till mm10-transkriptomet och UMI-kollapsning (Unique Molecular Identifier) utfördes med hjälp av Cellranger-verktygslådan (version 1.0.1) som tillhandahölls av 10X Genomics. För varje cell kvantifierade vi antalet gener för vilka minst en avläsning kartlades och exkluderade sedan alla celler med färre än 800 detekterade gener. Uttrycksvärden Ei,j för gen i i cell j beräknades genom att dividera UMI-räkningar för gen i med summan av UMI-räkningarna i cell j, för att normalisera för skillnader i täckning, och sedan multiplicera med 10 000 för att skapa TPM-liknande värden, och slutligen beräkna log2(TPM + 1). Batchkorrigering utfördes med ComBat46 som implementerats i R-paketet sva47, med standardläget för parametrisk justering. Utgången var en korrigerad uttrycksmatris, som användes som indata för vidare analys.

Selektion av variabla gener utfördes genom att anpassa en generaliserad linjär modell till förhållandet mellan den kvadrerade variationskoefficienten och den genomsnittliga uttrycksnivån i logaritmiskt utrymme, och genom att välja ut gener som avvek signifikant (P < 0,05) från den anpassade kurvan48.

Förbehandling av SMART-Seq2 scRNA-seq-data. BAM-filer konverterades till sammanslagna, avmultiplexade FASTQs med hjälp av Illuminas programvarupaket Bcl2Fastq v2.17.1.14. Paired-end reads mappades till UCSC mm10 musens transkriptom med hjälp av Bowtie49 med parametrarna ”-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6”, vilket möjliggör anpassning av sekvenser med en felmatchning. Genernas uttrycksnivåer kvantifierades med hjälp av TPM-värden som beräknades av RSEM50 v1.2.3 i parvis slutläge. För varje cell kvantifierade vi antalet gener för vilka minst en läsning kartlades och exkluderade sedan alla celler med antingen färre än 3 000 detekterade gener eller en transkriptomkartläggning på mindre än 40 %. Vi identifierade sedan mycket variabla gener enligt beskrivningen ovan.

Dimensionalitetsreduktion med hjälp av PCA och t-SNE. Vi begränsade uttrycksmatrisen till de undergrupper av variabla gener och celler av hög kvalitet som nämns ovan och centrerade och skalade sedan värdena innan de matades in i huvudkomponentanalysen (PCA), som genomfördes med hjälp av R-funktionen prcomp från stats-paketet för SMART-seq2-dataset. För det droppbaserade datasetet använde vi en randomiserad approximation av PCA, som genomfördes med hjälp av funktionen rpca från R-paketet rsvd, med parametern k satt till 100. Denna approximation med låg rang användes eftersom den är flera storleksordningar snabbare att beräkna för mycket breda matriser. Med tanke på att många huvudkomponenter förklarar en mycket liten del av variansen kan signal-brusförhållandet förbättras avsevärt genom att välja en delmängd av n ”signifikanta” huvudkomponenter. Efter PCA identifierades signifikanta huvudkomponenter med hjälp av permutationstestet51 , som genomfördes med hjälp av funktionen permutationPA från R-paketet jackstraw. Detta test identifierade 13 och 15 signifikanta huvudkomponenter i 10X- och SMART-Seq2-datasetterna i fig. 1b respektive Extended Data fig. 2a. Poäng från endast dessa signifikanta huvudkomponenter användes som indata för vidare analys.

För visualisering reducerades datasetens dimensionalitet ytterligare med hjälp av den ungefärliga versionen ”Barnes-hut” av t-SNE52,53. Detta implementerades med hjälp av funktionen Rtsne från Rtsne R-paketet med 20 000 iterationer och en perplexitetsinställning som varierade från 10 till 30 beroende på datasetets storlek.

Identifiera celldifferentieringsbanor med hjälp av diffusionskartor

För att köra dimensionalitetsreducering av diffusionskartor valde vi ut mycket variabla gener i data på följande sätt. Vi anpassade först en nollmodell för baslinjens cell-cell genuttrycksvariabilitet i data, med hjälp av ett power-law-förhållande mellan variationskoefficienten och medelvärdet av UMI-räkningarna för alla uttryckta gener, i likhet med tidigare arbete54. Därefter beräknade vi för varje gen skillnaden mellan värdet av dess observerade variationskoefficient och det värde som förväntas av nollmodellen (CVdiff). Histogrammet för CVdiff uppvisade en ”fet” svans. Vi beräknade medelvärdet μ och standardavvikelsen σ för denna fördelning och valde ut alla gener för vilka CVdiff > μ + 1,67σ, vilket gav 761 gener för vidare analys.

Vi reducerade dimensionaliteten med hjälp av diffusionsmappmetoden22. I korthet beräknades en cell-cellövergångsmatris med hjälp av en gaussisk kärna, där kärnans bredd anpassades till varje cells lokala grannskap55. Denna matris omvandlades till en markovsk matris efter normalisering. De högra egenvektorerna vi (i = 0, 1, 2, …) i denna matris beräknades och sorterades i ordning efter minskande egenvärde λi (i = 0, 1, 2, …), efter att ha uteslutit den ”översta” egenvektorn v0, som motsvarar λ0 = 1 (vilket avspeglar normaliseringsbegränsningen i den markovska matrisen). De återstående egenvektorerna vi (i = 1, 2, …) definierar inbäddningen av diffusionskartan och kallas diffusionskomponenter (DCk, k = 1, 2, …). Vi noterade ett spektralt gap mellan λ4 och λ5 och behöll därför DC1-DC4 för både det ursprungliga datasetet (Extended Data Fig. 4) och de data som extraherats från distinkta tarmregioner (Fig. 2c).

Fråntagande av kontaminerande immunceller och dubbletter

Trots att cellerna sorterades före sekvensering med hjälp av EpCAM, observerades ett litet antal kontaminerande immunceller i 10X-datasetet. Dessa 264 celler avlägsnades genom en första omgång av oövervakad klustring (densitetsbaserad klustring av t-SNE-kartan med hjälp av dbscan56 från R-paketet fpc) eftersom de bildade ett extremt distinkt kluster. För SMART-Seq2-datasetet var flera celler utfallande när det gäller bibliotekskomplexitet, vilket möjligen kan motsvara mer än en enskild cell per sekvenseringsbibliotek (”dubbletter”). Dessa celler togs sedan bort genom att beräkna den övre 1 %-kvantilen av fördelningen av gener som upptäcktes per cell och ta bort alla celler i denna kvantil.

Klusteranalys

För att klustra enskilda celler efter deras uttryck använde vi en oövervakad klusteransats, baserad på Infomap-grafklusteralgoritmen9 , som följer ansatser för CyTOF-data från enstaka celler57 och scRNA-seq10. I korthet konstruerade vi en k-närmaste granne-graf på data genom att för varje cellpar använda det euklidiska avståndet mellan poängen för de signifikanta huvudkomponenterna för att identifiera k närmaste grannar. Parametern k valdes för att vara förenlig med datasetets storlek. Specifikt sattes k till 200 och 80 för det droppbaserade datasetet med 7 216 celler (fig. 1b) respektive för SMART-Seq2-dataset med 1 522 celler (Extended Data fig. 2a). RANKL-behandlade organoider innehöll 5 434 celler och k sattes till 200. Datasetetet för Salmonella och H. polygyrus innehöll 9 842 celler och k sattes till 500. För klusteranalyser inom celltyper, särskilt de enteroendokrina och tuftcellsundergrupperna, använde vi Pearsonkorrelationsavståndet i stället för det euklidiska avståndet och ställde in k = 15, k = 30 och k = 40 för de enteroendokrina undertyperna (533 celler) och för de 166 och 102 tuftcellerna i 10X- och SMART-Seq2-dataseten, respektive. Grafen för närmaste grannar beräknades med hjälp av funktionen nng från R-paketet cccd. Grafen med k närmaste grannar användes sedan som indata till Infomap9, som implementerades med hjälp av funktionen infomap.community från R-paketet igraph.

Detekterade kluster kartlades till celltyper eller intermediära tillstånd med hjälp av kända markörer för tarmepitelcellsubtyper. (Extended Data Fig. 1g, Extended Data Fig. 2a). För den enteroendokrina (EEC) cellunderanalysen (fig. 3) slogs alla grupper av EEC-progenitorkluster med genomsnittliga parvisa korrelationer mellan signifikanta huvudkomponentpoäng på r > 0,85 samman, vilket resulterade i fyra kluster. Vi betecknade dessa fyra kluster som progenitor ”A” på grundval av höga nivåer av Ghrl, eller progenitor (tidig), (medel) eller (sen) (i den ordningen) på grundval av minskande nivåer av stam- (Slc12a2, Ascl2, Axin2) och cellcykelgener och ökande nivåer av kända EEG-regulatoriska faktorer (Neurod1, Neurod2 och Neurog3) (Extended Data Fig. 5c, Supplementary Table 6). För SMART-Seq2-dataset slogs två kluster som uttrycker höga nivåer av stamcellsmarkörgener (Extended Data Fig. 2a) samman för att bilda ett ”stam”-kluster och två andra kluster slogs samman för att bilda ett ”TA”-kluster.

För klusteranalysen av datasetet för follikelassocierat epitel med 4 700 celler var mikroveckcellerna mycket sällsynta (0,38 %) och därför användes ClusterDP-metoden58 för att identifiera dem eftersom den empiriskt sett presterade bättre än algoritmen k-nearest-neighbour-graf på detta dataset. Precis som med k-nearest-neighbour-metoderna kördes ClusterDP med signifikanta (P < 0,05) huvudkomponentpoäng (19 i det här fallet) som indata, och implementerades med hjälp av funktionerna findClusters och densityClust från R-paketet densityClust med parametrarna rho = 1.1 och delta = 0,25.

Extrahera sällsynta celltyper för vidare analys

Den inledande klusterbildningen av datasetet för hela tarmen (7 216 celler; fig. 1b) visade ett kluster med 310 EEC-celler och 166 tuftceller. Tufftcellerna togs ”som de var” för underanalysen (fig. 4a, b), medan EEC-cellerna kombinerades med ett andra kluster av 239 EEC-celler som identifierades i det regionala datasetet (fig. 2a, höger) för att få totalt 549 EEC-celler. En grupp på 16 celler gav uttryck för EEC-markörerna Chga och Chgb tillsammans med markörer för Paneth-celler, inklusive Lyz1, Defa5 och Defa22, och tolkades därför som dubbleringar och togs bort från analysen, varvid 533 EEC-celler återstod, som låg till grund för analysen i fig. 3. För att jämföra uttrycksprofiler för enterocyter från den proximala och distala tunntarmen (Fig. 2b) användes de 1 041 enterocyter som identifierats från 11 665 celler i det regionala datasetet (Fig. 2a).

Definiering av signaturer för celltyper

För att identifiera maximalt specifika gener för celltyper körde vi differentialuttryckstester mellan varje par av kluster för alla möjliga parvisa jämförelser. För ett givet kluster filtrerades sedan de förmodade signaturgenerna med hjälp av det maximala FDR Q-värdet och rangordnades efter minsta log2(fold change). Den minsta vikförändringen och det maximala Q-värdet representerar den svagaste effektstorleken över alla parvisa jämförelser; det är därför ett strängt kriterium. Celltypssignaturgener som visas i fig. 1c, Extended Data fig. 2b, Extended Data fig. 8e och de kompletterande tabellerna 2-4 och 8 erhölls med hjälp av en maximal FDR på 0,05 och en minsta log2(fold change) på 0,5. När det gäller signaturer för postmitotiska celltyper passerade alla gener detta tröskelvärde i både 3′ (fig. 1c) och fullständiga dataset (Extended Data fig. 2b).

I fallet med signaturgener för subtyper inom celltyper (fig. 3b, fig. 4b, Extended Data fig. 7b) beräknades ett kombinerat P-värde (över de parvisa testerna) för anrikning med hjälp av Fishers metod – ett mildare kriterium än att helt enkelt ta det maximala P-värdet – och ett maximalt FDR Q-värde på 0,01 användes, tillsammans med ett gränsvärde för minsta log2(fold change) på 0,25 för subtyper av tuftceller (Fig. 4b, Extended Data Fig. 7b, Supplementary Table 7) och på 0,1 för EEC-subtyper (Fig. 3b, Supplementary Table 6). Alla gener i tuftcellssignaturen passerade detta gränsvärde i både 3′ (fig. 4b) och fullständiga dataset (Extended Data Fig. 7b), medan signaturer för EEC-subtyper definierades endast med hjälp av 3′. På grund av det låga antalet celler (n = 18) användes också Fishers kombinerade P-värde för in vivo mikroveckcellssignaturen, med en FDR-gräns på 0,001 (fig. 5d, kompletterande tabell 8). Markörgener rangordnades efter minsta log2(fold change). Test av differentiellt uttryck utfördes med Mann-Whitney U-testet (även känt som Wilcoxon rank-sum-testet) som implementerades med R-funktionen wilcox.test. För infektionsexperimenten (fig. 6) använde vi en tvådelad ”hindermodell” för att kontrollera både den tekniska kvaliteten och variationen från mus till mus. Denna implementerades med hjälp av R-paketet MAST59, och P-värden för differentiellt uttryck beräknades med hjälp av likelihood-ratio-testet. Korrigering av multipel hypotesprövning utfördes genom att kontrollera FDR60 med hjälp av R-funktionen p.adjust.

Scoring av celler med hjälp av signaturgenuppsättningar

För att erhålla en poäng för en specifik uppsättning n gener i en given cell definierades en ”bakgrundsgenuppsättning” för att kontrollera skillnaderna i sekvenseringstäckning och bibliotekskomplexitet mellan cellerna på ett liknande sätt som i ref. 12. Bakgrundsgenuppsättningen valdes ut för att likna de intressanta generna när det gäller uttrycksnivå. Närmare bestämt valdes de 10n närmaste grannarna i det tvådimensionella rum som definieras av medeluttryck och detektionsfrekvens i alla celler. Signaturpoängen för den cellen definierades sedan som medeluttrycket av de n signaturgenerna i den cellen, minus medeluttrycket av de 10n bakgrundsgenerna i den cellen.

Skattningar av provtagningsfrekvenser för celltyper

För varje celltyp modelleras sannolikheten för att observera minst n celler i ett prov av storleken k med hjälp av den kumulativa fördelningsfunktionen för en negativ binomial NBcdf(k, n, p), där p är den relativa abundansen av denna celltyp. För m celltyper med samma parameter p är den totala sannolikheten att se varje typ minst n gånger NBcdf(k; n, p)m. En sådan analys kan utföras med användarspecificerade parametrar på http://satijalab.org/howmanycells.

EEC dendrogram

Genomsnittliga uttrycksvektorer beräknades för alla 12 EEC-undergruppskluster, med hjälp av log2(TPM + 1)-värden, och begränsades till den delmängd av 1 361 gener som identifierades som signifikant varierande mellan EEC-undergrupperna (P < 0,05), enligt beskrivning ovan. De genomsnittliga uttrycksvektorerna inklusive dessa gener klustrades hierarkiskt med hjälp av R-paketet pvclust (Spearman-distans, ward.D2-klustermetod), som ger bootstrap-konfidensskattningar för varje dendrogramknut som ett empiriskt P-värde över 100 000 försök (Extended Data Fig. 6a).

Celltypspecifika transkriptionsfaktorer, GPCR:er och proteiner med leukinrika upprepningar

En förteckning över alla gener som identifierats som verksamma som transkriptionsfaktorer hos möss hämtades från AnimalTFDB61. Uppsättningen av GPCR:er hämtades från UniProt-databasen (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). Funktionella annotationer för varje protein (Extended Data Fig. 2d) hämtades från British Pharmacological Society (BPS) och International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). Förteckningen över proteiner med leukinrika upprepningar hämtades från ref. 62. För att kartlägga gennamnen från människa till mus hämtades ortologer från människa och mus från Ensembl (senaste versionen 86; http://www.ensembl.org/biomart/martview), och synonymer till gener från människa och mus från NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). För varje mänsklig leucinrik upprepningsgen kartlades alla mänskliga synonymer till den ortologa genen i musen med hjälp av ortologförteckningen, och musens gennamn kartlades till namnen i encellsdata med hjälp av synonymförteckningen.

Celltypberikade transkriptionsfaktorer, GPCR:er och leucinrika upprepningsproteiner identifierades sedan genom att korsa listan över gener som är berikade i varje celltyp med listorna över transkriptionsfaktorer, GPCR:er och leucinrika upprepningsproteiner som definierats ovan. Celltypberikade gener definierades med hjälp av SMART-Seq2-dataset som gener med en minsta log2(fold change) på 0 och en högsta FDR på 0,5, med högst 10 gener per celltyp i Extended Data Fig. 2e, f (fullständiga listor finns i Supplementary Table 5). Dessutom identifierades en mer omfattande panel av celltypspecifika GPCR:er (Extended Data Fig. 2d) genom att välja ett mildare tröskelvärde. Detta uppnåddes genom att jämföra varje celltyp med alla andra celler, i stället för de parvisa jämförelser som beskrivs i föregående avsnitt, och genom att välja alla GPCR-gener som var differentiellt uttryckta (FDR < 0,001).

Testning av förändringar i celltypsandelar

Vi modellerar det påvisade antalet av varje celltyp i varje analyserad mus som en slumpmässig räknevariabel med hjälp av en Poissonprocess. Upptäcktshastigheten modelleras sedan genom att tillhandahålla det totala antalet celler som profileras i en viss mus som en offsetvariabel, med varje mus tillstånd (behandling eller kontroll) som en kovariabel. Modellen anpassades med hjälp av R-kommandot glm från stats-paketet. P-värdet för betydelsen av den effekt som behandlingen gav upphov till bedömdes med hjälp av ett Wald-test på regressionskoefficienten.

För bedömningen av betydelsen av rumsliga fördelningar av EEC-undergrupper (fig. 3e) omfattade jämförelsen mer än två grupper. I synnerhet var vår nollhypotes att andelen av varje EEC-undergrupp som upptäcktes i de tre tarmregionerna (duodenum, jejunum och ileum) var lika stor. För att testa denna hypotes använde vi variansanalys (ANOVA) med ett χ2-test på Poisson-modellens anpassning som beskrivs ovan, implementerad med hjälp av funktionen anova från stats-paketet.

Anrikning av genuppsättningar och analys av genontologi

Genontologianalysen utfördes med hjälp av R-paketet goseq63, med hjälp av signifikant differentiellt uttryckta gener (FDR < 0.05) som målgener, och alla gener uttryckta med log2(TPM + 1) > 3 i minst tio celler som bakgrund.

Datatillgänglighet

Kodtillgänglighet