Jednobuněčný průzkum epitelu tenkého střeva

Myši

Všechny práce na myších byly provedeny v souladu s institucionálními výbory pro péči o zvířata a jejich použití (IACUC) a s příslušnými pokyny Broad Institute a Massachusetts Institute of Technology, s protokoly 0055-05-15, resp. 0612-058-18. Pro všechny experimenty byly myši náhodně přiřazeny do léčebných skupin po porovnání pohlaví a věku 7-10týdenních samic nebo samců divokého typu C57BL/6J nebo myší Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP), získaných od Jackson Laboratory (Bar Harbour), nebo myší Gfi1beGFP/+ (Gfi1b-GFP)43 . Myši byly umístěny v podmínkách bez specifických patogenů ve zvířecích zařízeních Broad Institute, Massachusetts Institute of Technology nebo Harvard T. H. Chan School of Public Health.

Infekce salmonelou enterica a H. polygyrus. Myši C57BL/6J (Jackson Laboratory) byly infikovány 200 larvami třetího stadia H. polygyrus nebo 108 larvami Salmonella enterica, chovanými v podmínkách prostých specifických patogenů v Massachusetts General Hospital (Charlestown), podle protokolu 2003N000158. H. polygyrus byl rozmnožen podle předchozího popisu44. Myši byly utraceny 3 a 10 dní po infekci H. polygyrus. V případě Salmonella enterica byly myši infikovány přirozeně streptomycin-rezistentním kmenem S. Typhimurium SL1344 (108 buněk), jak bylo popsáno dříve44, a byly eutanazovány 48 hodin po infekci.

Disociace buněk a izolace krypt

Izolace krypt. Tenké střevo myší C57BL/6J divokého typu, Lgr5-GFP nebo Gfi1b-GFP bylo izolováno a opláchnuto ve studeném PBS. Tkáň byla podélně otevřena a nakrájena na malé fragmenty o délce zhruba 2 mm. Tkáň byla inkubována ve 20 mM EDTA-PBS na ledu po dobu 90 minut, přičemž se každých 30 minut protřepala. Tkáň byla poté silně protřepána a supernatant byl odebrán jako frakce 1 do nové kónické zkumavky. Tkáň byla inkubována v čerstvém EDTA-PBS a každých 30 min byla odebrána nová frakce. Frakce se sbíraly, dokud supernatant neobsahoval téměř výhradně krypty. Poslední frakce (obohacená o krypty) byla dvakrát promyta v PBS, centrifugována při 300 g po dobu 3 min a disociována pomocí TrypLE express (Invitrogen) po dobu 1 min při 37 °C. Suspenze jednotlivých buněk byla poté propasírována přes 40 μm filtr a obarvena pro FACS pro scRNA-seq (níže) nebo použita pro organoidní kulturu. Robustnost této metody jsme potvrdili testováním dalších metod izolace jednotlivých buněk – buď „celých“ (seškrábnutím epitelové výstelky), nebo „obohacených vilek“ (frakce 1; viz výše) – a zjistili jsme, že vzhledem k vysoké úmrtnosti (prostřednictvím anoikis) postmitotických diferencovaných buněk (jejichž primární složkou jsou zralé enterocyty) suspenze jednotlivých buněk obohacená kryptami věrně reprezentuje složení typů buněk tenkého střeva (údaje nejsou uvedeny).

Izolace epitelů spojených s folikuly. Epitelové buňky z epitelů asociovaných s folikuly byly izolovány extrakcí malých řezů (0,2-0,5 cm) obsahujících Peyerovy skvrny z tenkého střeva myší C57Bl/6J nebo Gfi1beGFP/+.

Třídění buněk

Pro experimenty s plnohodnotnou scRNA-seq na destičkách byl použit přístroj FACS (Astrios) k třídění jedné buňky do každé jamky 96jamkové destičky PCR obsahující 5 μl TCL pufru s 1% 2-merkaptoetanolem. Pro izolaci EpCAM+ byly buňky barveny na 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience); pro specifické epiteliální buňky jsme také barvili na CD24+/- (eBioscience) a c-Kit+/- (eBioscience). Pro obohacení specifických populací střevních epiteliálních buněk byly izolovány buňky z myší Lgr5-GFP, obarveny výše uvedenými protilátkami a označeny GFP-high (kmenové buňky), GFP-low (TA), GFP-/CD24+/c-Kit+/- (sekreční linie) nebo GFP-/CD24-/EpCAM+ (epiteliální buňky). Pro lepší obnovu Panethových buněk jsme povolili vyšší parametry bočního rozptylu a předního rozptylu v kombinaci s CD24+/c-Kit+ k ověření obnovy Panethových buněk u EpCAM+ buněk. Pro izolaci tuft-2 byly epiteliální buňky ze tří různých myší obarveny stejně jako výše, ale s použitím EpCAM+/CD45+ pro třídění 2 000 jednotlivých buněk. Použili jsme mírnou třídicí bránu, abychom zajistili získání dostatečného počtu těchto vzácných buněk tuft-2, což vedlo k vyšší míře kontaminace T-buňkami, které jsme odstranili v naší analýze jednotlivých buněk pomocí neřízeného shlukování.

Pro třídění scRNA plné délky byla 96jamková destička pevně uzavřena těsněním Microseal F a centrifugována při 800 g po dobu 1 min. Destička byla okamžitě zmrazena na suchém ledu a uchovávána při teplotě -80 °C, dokud nebyla připravena k vyčištění lyzátu. Buňky hromadné populace byly roztříděny do Eppendorfovy zkumavky obsahující 100 μl roztoku TCL s 1% 2-merkaptoethanolem a uloženy při -80 °C.

Pro scRNA-seq na bázi kapek byly buňky tříděny se stejnými parametry jako pro scRNA-seq na bázi destiček, ale byly tříděny do Eppendorfovy zkumavky obsahující 50 μl 0,4% BSA-PBS a uchovávány na ledu až do přechodu na platformu GemCode pro jednotlivé buňky.

ScRNA-seq na bázi destiček

Jednotlivé buňky. Knihovny byly připraveny pomocí upraveného protokolu SMART-Seq216. Stručně řečeno, vyčištění lyzátu RNA bylo provedeno pomocí kuliček RNAClean XP (Agencourt), následovala reverzní transkripce pomocí reverzní transkriptázy Maxima (Life Technologies) a celopřepisová amplifikace (WTA) pomocí směsi KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) po dobu 21 cyklů. Produkty WTA byly přečištěny pomocí kuliček Ampure XP (Beckman Coulter), kvantifikovány pomocí Qubit dsDNA HS Assay Kit (ThermoFisher) a vyhodnoceny pomocí vysoce citlivého DNA čipu (Agilent). Knihovny RNA-seq byly sestaveny z purifikovaných produktů WTA pomocí sady Nextera XT DNA Library Preperation Kit (Illumina). Na každé destičce byly populační kontroly a kontroly bez buněk zpracovány stejnou metodou jako u jednotlivých buněk. Knihovny byly sekvenovány na přístroji Illumina NextSeq 500.

Skupinové vzorky. Hromadné populační vzorky byly zpracovány extrakcí RNA pomocí sady RNeasy Plus Micro Kit (Qiagen) podle doporučení výrobce a následným postupem podle upraveného protokolu SMART-Seq2 po vyčištění lyzátu, jak je popsáno výše.

Sekvence scRNA na bázi kapek

Jednotlivé buňky byly zpracovány prostřednictvím platformy GemCode Single Cell Platform pomocí sad GemCode Gel Bead, Chip a Library Kits (10X Genomics, Pleasanton) podle protokolu výrobce. Stručně řečeno, jednotlivé buňky byly tříděny do 0,4% BSA-PBS. Do každého kanálu bylo přidáno 6 000 buněk s průměrnou mírou obnovy 1 500 buněk. Buňky byly poté rozděleny do gelových kuliček v emulzi v přístroji GemCode, kde proběhla lýza buněk a reverzní transkripce RNA s čárovým kódem, následovaná amplifikací, střihem a připojením 5′ adaptoru a indexu vzorku. Knihovny byly sekvenovány na přístroji Illumina NextSeq 500.

Imnofluorescence a smFISH

Imnofluorescence. Barvení tkání tenkého střeva bylo provedeno podle předchozího popisu34. Stručně řečeno, tkáně byly fixovány po dobu 14 hodin ve formalínu, vloženy do parafínu a nařezány na řezy o tloušťce 5 μm. Řezy byly deparafinovány standardními technikami, inkubovány s primárními protilátkami přes noc při 4 °C a poté se sekundárními protilátkami při pokojové teplotě po dobu 30 min. Sklíčka byla upevněna pomocí přípravku Slowfade Mountant + DAPI (Life Technologies, S36964) a zapečetěna.

smFISH. Byla použita sada RNAScope Multiplex Flourescent Kit (Advanced Cell Diagnostics) podle doporučení výrobce s následujícími změnami. Doba varu pro vyhledání cíle byla upravena na 12 min a inkubace s proteázou IV při 40 °C byla upravena na 8 min. Sklíčka byla upevněna pomocí přípravku Slowfade Mountant+DAPI (Life Technologies, S36964) a zapečetěna.

Kombinovaná imunofluorescence a smFISH. Tato metoda byla provedena tak, že se nejprve provedla smFISH, jak je popsáno výše, s následujícími změnami. Po Amp 4 byly řezy tkání promyty v promývacím pufru, inkubovány s primárními protilátkami přes noc při 4 °C, třikrát promyty v 1× TBST a poté inkubovány se sekundárními protilátkami po dobu 30 min při pokojové teplotě. Sklíčka byla upevněna přípravkem Slowfade Mountant + DAPI (Life Technologies, S36964) a zapečetěna.

Image analysis

Snímky tkáňových řezů byly pořízeny konfokálním mikroskopem Fluorview FV1200 s použitím Kalmanova a sekvenčního laserového záření pro snížení šumu a překrývání signálů. Měřítka byla ke každému snímku přidána pomocí konfokálního softwaru FV10-ASW 3.1 Viewer. Obrázky byly překryty a vizualizovány pomocí softwaru Image J45.

Protilátky a sondy

Střevní organoidní kultury

Po izolaci krypt byla jednobuněčná suspenze resuspendována v Matrigelu (BD Bioscience) s 1 μM peptidu Jagged-1 (Ana-Spec). Do každé jamky 24jamkové destičky bylo vyseto zhruba 300 krypt zalitých v 25 μl Matrigela. Po ztuhnutí byl Matrigel inkubován v 600 μl kultivačního média (Advanced DMEM/F12, Invitrogen) se streptomycinem/penicilinem a glutamátem a doplněn o EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), Y-276432 dihydrochlorid monohydrát (10 μM, Tochris), N-acetyl-1-cystein (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) a Wnt3A (25 ng ml-1, R&D Systems). Čerstvá média byla vyměněna 3. den a organoidy byly pasážovány disociací pomocí TrypLE a 6. den znovu suspendovány v novém Matrigelu s poměrem rozdělení 1:3. U vybraných experimentů byly organoidy navíc ošetřeny RANKL (100 ng ml-1, Biolegends). Ošetřené organoidy byly disociovány a podrobeny scRNA-seq oběma metodami.

Kvantitativní PCR

Pro relativní qPCR byla použita cDNA 16 jednotlivých buněk tuft-1, tuft-2 a náhodných EpCam+ z destiček založených na scRNA-seq v plné délce. Exprese genů byla analyzována pomocí kvantitativní PCR v reálném čase na přístroji LightCycler 480 Instrument II (Roche) s použitím směsi LightCycler 480 SYBR green (Roche) s následujícími sadami primerů: HPRT1-F, GTTAAGCAGTACAGCCCCAAA; HPRT1-R, AGGGCATATCCAACAACAAACTT; UBC-F, CAGCCGTATCTTCCCAGACT; UBC-R, CTCAGAGGGATGCCAGTAATCTA; tslp-F, TACTCTCAATCCTATCCCTGCTG; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC; Alpi-F, TCCTACCTCCATTCTCTATGG, Alpi-R, CCGCCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCTACCATC; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. Primery qPCR byly u všech transkriptů navrženy pro exon-exonovou hranici.

Výpočetní analýza

Předběžné zpracování dat scRNA-seq z kapek. De-multiplexování, zarovnání k transkriptomu mm10 a přiřazení jedinečného molekulárního identifikátoru (UMI) bylo provedeno pomocí sady nástrojů Cellranger (verze 1.0.1) poskytnuté společností 10X Genomics. U každé buňky jsme kvantifikovali počet genů, pro které byl mapován alespoň jeden čtení, a poté jsme vyloučili všechny buňky s méně než 800 detekovanými geny. Hodnoty exprese Ei,j pro gen i v buňce j byly vypočteny vydělením počtu UMI pro gen i součtem počtu UMI v buňce j, aby se normalizovaly rozdíly v pokrytí, a poté vynásobením 10 000 pro vytvoření hodnot podobných TPM a nakonec výpočtem log2(TPM + 1). Dávková korekce byla provedena pomocí programu ComBat46 implementovaného v balíčku R sva47 s použitím výchozího parametrického režimu úpravy. Výstupem byla opravená expresní matice, která byla použita jako vstup pro další analýzu.

Výběr variabilních genů byl proveden fitováním zobecněného lineárního modelu na vztah mezi čtvercovým variačním koeficientem a střední úrovní exprese v logaritmickém prostoru a výběrem genů, které se významně odchýlily (P < 0,05) od fitované křivky48.

Předběžné zpracování dat scRNA-seq SMART-Seq2. Soubory BAM byly převedeny na sloučené, de-multiplexované soubory FASTQ pomocí softwarového balíčku Bcl2Fastq v2.17.1.14 dodaného společností Illumina. Párová čtení byla mapována na myší transkriptom UCSC mm10 pomocí Bowtie49 s parametry „-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6“, který umožňuje zarovnání sekvencí s jednou neshodou. Úrovně exprese genů byly kvantifikovány pomocí hodnot TPM vypočtených programem RSEM50 v1.2.3 v režimu párových konců. Pro každou buňku jsme kvantifikovali počet genů, pro které byl mapován alespoň jeden čtení, a poté jsme vyloučili všechny buňky s méně než 3 000 detekovanými geny nebo s mapováním transkriptomu menším než 40 %. Poté jsme identifikovali vysoce variabilní geny, jak je popsáno výše.

Snížení dimenzionality pomocí PCA a t-SNE. Expresní matici jsme omezili na výše uvedené podmnožiny variabilních genů a vysoce kvalitních buněk a poté jsme hodnoty centrovali a škálovali před jejich zadáním do analýzy hlavních komponent (PCA), která byla implementována pomocí funkce R prcomp z balíčku stats pro datovou sadu SMART-seq2. Pro soubor dat založený na kapičkách jsme použili náhodnou aproximaci PCA, implementovanou pomocí funkce rpca z balíčku rsvd R, s parametrem k nastaveným na 100. Tato aproximace s nízkou rangí byla použita, protože její výpočet je pro velmi široké matice o několik řádů rychlejší. Vzhledem k tomu, že mnoho hlavních komponent vysvětluje jen velmi malou část rozptylu, lze poměr signálu k šumu podstatně zlepšit výběrem podmnožiny n „významných“ hlavních komponent. Po PCA byly významné hlavní komponenty identifikovány pomocí permutačního testu51 , implementovaného pomocí funkce permutationPA z balíčku jackstraw R. Tento test identifikoval 13 a 15 významných hlavních komponent v souborech dat 10X na obr. 1b a SMART-Seq2 na obr. 2a s rozšířenými daty. Jako vstup pro další analýzu bylo použito skóre pouze z těchto významných hlavních komponent.

Pro vizualizaci byla dimenzionalita datových souborů dále redukována pomocí přibližné verze „Barnes-hut“ programu t-SNE52,53 . Tato metoda byla implementována pomocí funkce Rtsne z balíčku Rtsne R s použitím 20 000 iterací a nastavení perplexity, které se pohybovalo od 10 do 30 v závislosti na velikosti souboru dat.

Identifikace trajektorií buněčné diferenciace pomocí difuzních map

Před spuštěním redukce dimenzionality pomocí difuzních map jsme vybrali vysoce variabilní geny v datech následujícím způsobem. Nejprve jsme v datech dosadili nulový model pro základní variabilitu exprese genů mezi buňkami pomocí power-law vztahu mezi variačním koeficientem a průměrem počtu UMI všech exprimovaných genů, podobně jako v předchozí práci54. Poté jsme pro každý gen vypočítali rozdíl mezi hodnotou jeho pozorovaného variačního koeficientu a hodnotou očekávanou podle nulového modelu (CVdiff). Histogram CVdiff vykazoval „tlustý“ chvost. Vypočítali jsme průměr μ a směrodatnou odchylku σ tohoto rozdělení a vybrali jsme všechny geny, pro které CVdiff > μ + 1,67σ, čímž jsme získali 761 genů pro další analýzu.

Snížení dimenzionality jsme provedli pomocí metody difuzní mapy22. Stručně řečeno, matice přechodu mezi buňkami byla vypočtena pomocí Gaussova jádra, přičemž šířka jádra byla upravena podle lokálního okolí každé buňky55. Tato matice byla po normalizaci převedena na markovskou matici. Pravé vlastní vektory vi (i = 0, 1, 2, …) této matice byly vypočteny a seřazeny v pořadí klesající vlastní hodnoty λi (i = 0, 1, 2, …) po vyloučení „horního“ vlastního vektoru v0, který odpovídá λ0 = 1 (což odráží normalizační omezení markovské matice). Zbývající vlastní vektory vi (i = 1, 2, …) definují vložení difuzní mapy a označují se jako difuzní komponenty (DCk, k = 1, 2, …). Zaznamenali jsme spektrální mezeru mezi λ4 a λ5, a proto jsme zachovali DC1-DC4 jak pro původní soubor dat (Obr. 4 rozšířených dat), tak pro data získaná z odlišných oblastí střeva (Obr. 2c).

Odstranění kontaminujících imunitních buněk a dubletů

Ačkoli byly buňky před sekvenováním tříděny pomocí EpCAM, v souboru dat 10X bylo pozorováno malé množství kontaminujících imunitních buněk. Těchto 264 buněk bylo odstraněno úvodním kolem neřízeného shlukování (shlukování na základě hustoty mapy t-SNE pomocí dbscan56 z balíčku R fpc), protože tvořily extrémně výrazný shluk. U datové sady SMART-Seq2 bylo několik buněk odlehlých z hlediska složitosti knihovny, což by mohlo odpovídat více než jedné jednotlivé buňce na sekvenační knihovnu („dublety“). Tyto buňky byly poté odstraněny výpočtem horního 1% kvantilu distribuce genů detekovaných na buňku a odstraněním všech buněk v tomto kvantilu.

Klastrová analýza

Pro shlukování jednotlivých buněk podle jejich exprese jsme použili přístup neřízeného shlukování založený na algoritmu shlukování grafů Infomap9 podle přístupů pro jednobuněčná data CyTOF57 a scRNA-seq10. Stručně řečeno, na datech jsme vytvořili graf k nejbližších sousedů, přičemž jsme pro každou dvojici buněk použili euklidovskou vzdálenost mezi skóre významných hlavních komponent k určení k nejbližších sousedů. Parametr k byl zvolen tak, aby odpovídal velikosti souboru dat. Konkrétně byla hodnota k nastavena na 200 a 80 pro soubor dat založený na kapičkách o 7 216 buňkách (obr. 1b) a pro soubor dat SMART-Seq2 o 1 522 buňkách (rozšířená data, obr. 2a). Organoidy ošetřené RANKL obsahovaly 5 434 buněk a k bylo nastaveno na 200; datový soubor Salmonella a H. polygyrus obsahoval 9 842 buněk a k bylo nastaveno na 500. Pro shlukovou analýzu v rámci buněčných typů, konkrétně podskupin enteroendokrinních a chomáčových buněk, jsme místo euklidovské vzdálenosti použili Pearsonovu korelační vzdálenost a nastavili k = 15, k = 30 a k = 40 pro enteroendokrinní podtypy (533 buněk) a pro 166 a 102 chomáčových buněk v souborech dat 10X a SMART-Seq2. Graf nejbližšího souseda byl vypočten pomocí funkce nng z balíčku R cccd. Graf k-nejbližších sousedů byl poté použit jako vstup do programu Infomap9, implementovaného pomocí funkce infomap.community z balíčku igraph R.

Zjištěné shluky byly mapovány na typy buněk nebo mezistavy pomocí známých markerů pro subtypy střevních epiteliálních buněk. (Obr. 1g rozšířených dat, obr. 2a rozšířených dat). Pro subanalýzu enteroendokrinních buněk (EEC) (obr. 3) byla sloučena jakákoli skupina klastrů progenitorů EEC s průměrnými párovými korelacemi mezi významnými skóre hlavních komponent r > 0,85, čímž vznikly čtyři klastry. Tyto čtyři shluky jsme označili jako progenitorové „A“ na základě vysokých hladin Ghrl nebo progenitorové (časné), (střední) nebo (pozdní) (v tomto pořadí) na základě klesajících hladin kmenových genů (Slc12a2, Ascl2, Axin2) a genů buněčného cyklu a rostoucích hladin známých regulačních faktorů EEC (Neurod1, Neurod2 a Neurog3) (Extended Data Fig. 5c, Supplementary Table 6). U datové sady SMART-Seq2 byly dva shluky exprimující vysoké hladiny markerových genů kmenových buněk (Extended Data Obr. 2a) sloučeny do „kmenového“ shluku a dva další shluky byly sloučeny do shluku „TA“.

Pro shlukovou analýzu datové sady 4 700 buněk folikulárního epitelu byly mikroskopické buňky velmi vzácné (0,38 %), a proto byla k jejich identifikaci použita metoda ClusterDP58 , protože na této datové sadě vykazovala empiricky lepší výsledky než algoritmus k-nejbližšího souseda grafu. Stejně jako u metody k-nejbližšího souseda byla metoda ClusterDP spuštěna s použitím významných (P < 0,05) skóre hlavních komponent (v tomto případě 19) jako vstupu a byla implementována pomocí funkcí findClusters a densityClusters z balíčku densityClust R s použitím parametrů rho = 1. V tomto případě byly jako vstupy použity skóre hlavních komponent (19).1 a delta = 0,25.

Vyčlenění vzácných typů buněk pro další analýzu

Počáteční shlukování souboru dat z celého střeva (7 216 buněk; obr. 1b) ukázalo shluk 310 buněk EEC a 166 buněk chomáčků. Buňky chomáčků byly pro dílčí analýzu převzaty „tak, jak jsou“ (obr. 4a, b), zatímco buňky EEC byly spojeny s druhým shlukem 239 buněk EEC, které byly identifikovány v regionálním souboru dat (obr. 2a, vpravo), celkem tedy 549 buněk EEC. Skupina 16 buněk koexprimovala markery EEC Chga a Chgb s markery Panethových buněk, včetně Lyz1, Defa5 a Defa22, a byla proto interpretována jako dublety a z analýzy odstraněna, takže zůstalo 533 buněk EEC, které byly základem pro analýzu na obr. 3. Pro porovnání expresních profilů enterocytů z proximálního a distálního tenkého střeva (obr. 2b) bylo použito 1 041 enterocytů identifikovaných z 11 665 buněk v regionálním souboru dat (obr. 2a).

Definice signatur buněčných typů

Pro identifikaci maximálně specifických genů pro buněčné typy jsme provedli testy rozdílné exprese mezi každou dvojicí klastrů pro všechna možná párová srovnání. Poté byly pro daný klastr domnělé signaturní geny filtrovány pomocí maximální hodnoty Q FDR a seřazeny podle minimální log2(fold change). Minimální fold change a maximální hodnota Q představují nejslabší velikost účinku ve všech párových srovnáních; jedná se tedy o přísné kritérium. Geny signatury buněčného typu zobrazené na obr. 1c, v rozšířených datech na obr. 2b, v rozšířených datech na obr. 8e a v doplňkových tabulkách 2-4 a 8 byly získány s použitím maximální hodnoty FDR 0,05 a minimální hodnoty log2(fold change) 0,5. V případě signatur postmitotických buněčných typů prošly tímto prahem všechny geny jak v souborech dat 3′ (obr. 1c), tak v souborech dat plné délky (Extended Data obr. 2b).

V případě signaturních genů pro podtypy v rámci buněčných typů (obr. 3b, obr. 4b, Extended Data obr. 8b) byly všechny geny označeny jako signaturní. 7b) byla vypočtena kombinovaná hodnota P (napříč párovými testy) pro obohacení pomocí Fisherovy metody – což je mírnější kritérium než pouhé vzetí maximální hodnoty P – a byla použita maximální hodnota FDR Q 0,01 spolu s mezní hodnotou minimálního log2(fold change) 0,25 pro podtypy tuftových buněk (obr. 4b, Extended Data obr. 7b, Doplňková tabulka 7) a 0,1 pro podtypy EEC (obr. 3b, Doplňková tabulka 6). Všechny geny v signatuře tuft cell prošly touto hranicí jak v souborech dat 3′ (obr. 4b), tak v souborech dat plné délky (Extended Data obr. 7b), zatímco signatury podtypů EEC byly definovány pouze pomocí 3′. Vzhledem k nízkému počtu buněk (n = 18) byla Fisherova kombinovaná hodnota P použita také pro in vivo mikroskopickou buněčnou signaturu, přičemž mezní hodnota FDR byla 0,001 (obr. 5d, doplňková tabulka 8). Markerové geny byly seřazeny podle minimální log2(fold change). Testy rozdílné exprese byly provedeny pomocí Mannova-Whitneyho U testu (známého také jako Wilcoxonův rank-sum test) implementovaného pomocí funkce R wilcox.test. Pro experimenty s infekcí (obr. 6) jsme použili dvoudílný „překážkový“ model pro kontrolu technické kvality i variability mezi jednotlivými myšmi. Ten byl implementován pomocí balíčku R MAST59 a hodnoty P pro diferenciální expresi byly vypočteny pomocí testu poměru pravděpodobnosti. Korekce testování vícenásobných hypotéz byla provedena kontrolou FDR60 pomocí funkce R p.adjust.

Skórování buněk pomocí signaturních genových sad

Pro získání skóre pro specifickou sadu n genů v dané buňce byla definována „backgroundová“ genová sada pro kontrolu rozdílů v pokrytí sekvenování a složitosti knihovny mezi buňkami podobně jako v ref. 12. Soubor genů pozadí byl vybrán tak, aby byl podobný genům zájmu z hlediska úrovně exprese. Konkrétně bylo vybráno 10n nejbližších sousedů ve dvourozměrném prostoru definovaném průměrnou expresí a detekční frekvencí ve všech buňkách. Skóre signatury pro danou buňku pak bylo definováno jako průměrná exprese n signaturních genů v této buňce minus průměrná exprese 10n genů pozadí v této buňce.

Odhady četností vzorkování buněčných typů

Pro každý buněčný typ je pravděpodobnost pozorování alespoň n buněk ve vzorku o velikosti k modelována pomocí kumulativní distribuční funkce negativního binomu NBcdf(k, n, p), kde p je relativní četnost tohoto buněčného typu. Pro m typů buněk se stejným parametrem p je celková pravděpodobnost výskytu každého typu alespoň n-krát NBcdf(k; n, p)m. Takovou analýzu lze provést s uživatelem zadanými parametry na http://satijalab.org/howmanycells.

Dendrogram EEC

Průměrné vektory exprese byly vypočteny pro všech 12 podskupin EEC s použitím hodnot log2(TPM + 1) a omezeny na podskupinu 1 361 genů, které byly identifikovány jako významně variabilní mezi podskupinami EEC (P < 0,05), jak je popsáno výše. Průměrné expresní vektory zahrnující tyto geny byly hierarchicky shlukovány pomocí balíčku R pvclust (Spearmanova vzdálenost, metoda shlukování ward.D2), který poskytuje bootstrapový odhad spolehlivosti na každém uzlu dendrogramu jako empirickou hodnotu P na 100 000 pokusů (Extended Data Fig. 6a).

Buněčně specifické transkripční faktory, GPCR a proteiny bohaté na leucinové repetice

Seznam všech genů identifikovaných jako působící jako transkripční faktory u myší byl získán z databáze AnimalTFDB61. Soubor GPCR byl získán z databáze 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). Funkční anotace pro jednotlivé proteiny (rozšířená data obr. 2d) byly získány od Britské farmakologické společnosti (BPS) a Mezinárodní unie pro základní a klinickou farmakologii (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). Seznam proteinů bohatých na leucinové repetice byl převzat z ref. 62. Pro mapování z lidských na myší názvy genů byly lidské a myší ortology staženy z databáze Ensembl (poslední verze 86; http://www.ensembl.org/biomart/martview) a synonyma lidských a myších genů z NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). Pro každý lidský gen s opakováním bohatým na leucin byla všechna lidská synonyma mapována na ortologický gen u myši pomocí seznamu ortologů a názvy myších genů byly mapovány na názvy v datech o jedné buňce pomocí seznamu synonym.

Poté byly identifikovány transkripční faktory, GPCR a proteiny s opakováním bohatým na leucin obohacené na buněčný typ protnutím seznamu genů obohacených v každém buněčném typu se seznamy transkripčních faktorů, GPCR a proteinů s opakováním bohatým na leucin definovanými výše. Geny obohacené na buněčný typ byly definovány pomocí datové sady SMART-Seq2 jako geny s minimální log2(fold change) 0 a maximální FDR 0,5, přičemž bylo zachováno maximálně 10 genů na buněčný typ v rozšířených datech obr. 2e, f (kompletní seznamy jsou uvedeny v doplňkové tabulce 5). Kromě toho byl výběrem mírnějšího prahu identifikován rozsáhlejší panel GPCR specifických pro buněčný typ (Extended Data Obr. 2d). Toho bylo dosaženo porovnáním každého typu buněk se všemi ostatními buňkami namísto párových porovnání popsaných v předchozí části a výběrem všech genů GPCR, které byly diferenciálně exprimovány (FDR < 0,001).

Testování změn v zastoupení buněčných typů

Modelovali jsme zjištěný počet každého typu buněk u každé analyzované myši jako náhodnou počítanou proměnnou pomocí Poissonova procesu. Míra detekce je pak modelována uvedením celkového počtu profilovaných buněk u dané myši jako offsetové proměnné, přičemž stav každé myši (léčba nebo kontrola) je uveden jako kovariát. Model byl sestaven pomocí příkazu glm z balíčku stats. Hodnota P pro významnost účinku vyvolaného léčbou byla posouzena pomocí Waldova testu na regresním koeficientu.

Pro posouzení významnosti prostorového rozložení podskupin EEC (obr. 3e) zahrnovalo srovnání více než dvě skupiny. Naší nulovou hypotézou bylo zejména to, že podíl jednotlivých podskupin EEC zjištěných ve třech oblastech střeva (duodenum, jejunum a ileum) je stejný. K ověření této hypotézy jsme použili analýzu rozptylu (ANOVA) s χ2 testem na výše popsanou shodu s Poissonovým modelem, implementovanou pomocí funkce anova z balíčku stats.

Obohacení genových souborů a analýza genové ontologie

Analýza genové ontologie byla provedena pomocí balíčku goseq R63 s použitím významně diferencovaně exprimovaných genů (FDR < 0.05) jako cílové geny a všechny geny exprimované s log2(TPM + 1) > 3 v nejméně deseti buňkách jako pozadí.

Dostupnost dat

Dostupnost kódu

.