A vékonybélhám egysejtes felmérése

Egerek

Minden egérrel végzett munka a Broad Institute és a Massachusetts Institute of Technology intézményi állatgondozási és -használati bizottságainak (IACUC) és a vonatkozó irányelveknek megfelelően történt, a 0055-05-15 és a 0612-058-18 protokollok szerint. Minden kísérlethez az egereket véletlenszerűen osztották be a kezelési csoportokba, miután a Jackson Laboratory (Bar Harbour) vagy a Gfi1beGFP/+ (Gfi1b-GFP) egerekből43 származó 7-10 hetes vad típusú C57BL/6J vagy Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP) egerek nemének és korának megfeleltetése után 7-10 hetes nőstény vagy hím vad típusú egereket kaptak. Az egereket a Broad Institute, Massachusetts Institute of Technology vagy a Harvard T. H. Chan School of Public Health állattartó létesítményeiben, specifikus patogénmentes körülmények között tartották.

Salmonella enterica és H. polygyrus fertőzés. C57BL/6J egereket (Jackson Laboratory) 200 H. polygyrus harmadik stádiumú lárvával vagy 108 Salmonella enterica lárvával fertőztek meg, amelyeket a Massachusetts General Hospitalban (Charlestown) specifikus patogénmentes körülmények között tartottak, a 2003N000158 protokoll szerint. A H. polygyrus szaporítása a korábban leírtak szerint történt44. Az egereket 3 és 10 nappal a H. polygyrus fertőzés után elaltatták. A Salmonella enterica esetében az egereket a S. Typhimurium természetesen sztreptomicin-rezisztens SL1344 törzsével (108 sejt) fertőztük meg a korábban leírtak szerint44 , és 48 órával a fertőzés után elaltattuk őket.

Celldissociáció és kripták izolálása

Kripták izolálása. C57BL/6J vad típusú, Lgr5-GFP vagy Gfi1b-GFP egerek vékonybelét izoláltuk és hideg PBS-ben öblítettük. A szövetet hosszanti irányban felnyitottuk, és nagyjából 2 mm hosszúságú kis darabokra szeleteltük. A szövetet 20 mM EDTA-PBS-ben inkubáltuk jégen 90 percig, 30 percenként rázva. Ezután a szövetet erőteljesen felráztuk, és a felülúszót 1. frakcióként egy új kúpos csőbe gyűjtöttük. A szövetet friss EDTA-PBS-ben inkubáltuk, és 30 percenként új frakciót gyűjtöttünk. A frakciókat addig gyűjtöttük, amíg a felülúszó majdnem teljesen kriptákból állt. Az utolsó frakciót (a kripták számára dúsított) kétszer mostuk PBS-ben, 3 percig 300 g-nél centrifugáltuk, majd 1 percig 37 °C-on TrypLE express (Invitrogen) segítségével disszociáltuk. Az egysejtes szuszpenziót ezután 40 μm-es szűrőn átpasszíroztuk, és FACS festést végeztünk az scRNS-seq (alább) számára, vagy organoid tenyésztéshez használtuk. A módszer robusztusságát további egysejtes izolációs módszerek – akár “teljes” (a hámbélés lekaparása), akár “villusokkal dúsított” (1. frakció; lásd fentebb) – tesztelésével igazoltuk, és megállapítottuk, hogy a posztmitotikus differenciált sejtek (amelyek elsődleges összetevője az érett enterociták) magas mortalitási aránya miatt (anoikisz révén) a kriptával dúsított egysejtes szuszpenzió hűen reprezentálja a vékonybél sejttípusainak összetételét (az adatok nem láthatóak).

A tüszőkhöz kapcsolódó epithéliumok izolálása. A tüsző-asszociált epitéliák epitélsejtjeit C57Bl/6J vagy Gfi1beGFP/+ egerek vékonybeléből Peyer-foltokat tartalmazó kis metszetek (0,2-0,5 cm) kivonásával izoláltuk.

Cellaválogatás

A lemez alapú teljes hosszúságú scRNS-seq kísérletekhez egy FACS gépet (Astrios) használtunk, hogy egy-egy sejtet válogassunk egy 96 lyukú PCR lemez minden egyes lyukába, amely 5 μl TCL puffert tartalmazott 1% 2-merkaptoetanollal. Az EpCAM+ izoláláshoz a sejteket 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience) festéssel festettük; a specifikus epithelsejtek esetében CD24+/- (eBioscience) és c-Kit+/- (eBioscience) festést is végeztünk. A specifikus bélhámsejt-populációk dúsításához Lgr5-GFP egerekből izoláltunk sejteket, amelyeket a fent említett antitestekkel festettünk, és GFP-high (őssejtek), GFP-low (TA-k), GFP-/CD24+/c-Kit+/- (szekréciós vonalak) vagy GFP-/CD24-/EpCAM+ (epitélsejtek) gátlással jelöltünk. A Paneth-sejtek jobb visszanyerése érdekében a CD24+/c-Kit+ kombinációval magasabb oldalsó szórás és előre szórás paramétereket engedélyeztünk, hogy az EpCAM+ sejtek Paneth-sejt visszanyerését igazolni tudjuk. A tuft-2 izoláláshoz három különböző egérből származó epithelsejteket festettünk a fentiek szerint, de EpCAM+/CD45+ használatával 2000 egyedi sejtet válogattunk. Enyhe válogatási kaput alkalmaztunk annak érdekében, hogy elegendő számú ilyen ritka tuft-2 sejtet kapjunk, ami a T-sejtek magasabb kontaminációs arányához vezetett, amelyeket az egysejtes elemzés során a felügyelet nélküli klaszterezéssel eltávolítottunk.

A teljes hosszúságú scRNS-seq válogatáshoz a 96 lyukú lemezt szorosan lezártuk egy Microseal F-vel, és 1 percig 800 g-nél centrifugáltuk. A lemezt azonnal lefagyasztottuk szárazjégen, és -80 °C-on tartottuk a lizátum tisztításáig. Az ömlesztett populáció sejtjeit egy Eppendorf-csőbe válogattuk, amely 100 μl TCL-oldatot tartalmazott 1%-os 2-merkaptoetanollal, és -80 °C-on tároltuk.

A csepp alapú scRNS-seq esetében a sejteket ugyanazokkal a paraméterekkel válogattuk, mint a tányér alapú scRNS-seq esetében, de 50 μl 0,4%-os BSA-PBS-t tartalmazó Eppendorf-csőbe válogattuk, és jégen tároltuk, amíg a GemCode egysejtes platformra nem kerültek.

Tányér alapú scRNS-seq

Egysejtes sejtek. A könyvtárakat a módosított SMART-Seq2 protokoll16 segítségével állítottuk elő. Röviden, az RNS-lizátum tisztítását RNAClean XP gyöngyökkel (Agencourt) végeztük, amelyet reverz transzkripció követett Maxima reverz transzkriptázzal (Life Technologies) és teljes transzkripciós amplifikáció (WTA) KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) segítségével 21 cikluson keresztül. A WTA-termékeket Ampure XP gyöngyökkel (Beckman Coulter) tisztítottuk, a Qubit dsDNA HS Assay Kit (ThermoFisher) segítségével számszerűsítettük és nagy érzékenységű DNS-chippel (Agilent) értékeltük. A tisztított WTA-termékekből RNS-seq könyvtárakat készítettünk Nextera XT DNA Library Preperation Kit (Illumina) segítségével. Minden egyes lemezen a populációs és a sejt nélküli kontrollokat ugyanazzal a módszerrel dolgoztuk fel, mint az egyes sejtek esetében. A könyvtárak szekvenálása Illumina NextSeq 500-on történt.

Ömlesztett minták. Az ömlesztett populációs mintákat a gyártó ajánlásainak megfelelően az RNeasy Plus Micro Kit (Qiagen) segítségével történő RNS extrakcióval dolgoztuk fel, majd a lizátum tisztítását követően a fent leírtak szerint a módosított SMART-Seq2 protokollal folytattuk.

Cseppecske alapú scRNS-seq

Az egyes sejteket a GemCode Single Cell Platformon keresztül dolgoztuk fel a GemCode Gel Bead, Chip és Library Kits (10X Genomics, Pleasanton) segítségével a gyártó protokolljának megfelelően. Röviden, az egyes sejteket 0,4%-os BSA-PBS-be válogattuk. Minden csatornába 6000 sejtet adtunk, az átlagos visszanyerési arány 1500 sejt volt. A sejteket ezután a GemCode műszerben emulziós gélgyöngyökre osztottuk, ahol a sejtek lízise és az RNS vonalkódos reverz transzkripciója történt, majd az amplifikáció, a nyírás és az 5′ adaptor és a mintaindex rögzítése következett. A könyvtárak szekvenálása Illumina NextSeq 500 készülékkel történt.

Immunofluoreszcencia és smFISH

Immunofluoreszcencia. A vékonybélszövetek festését a korábban leírtak szerint végeztük34. Röviden, a szöveteket 14 órán át formalinban rögzítettük, paraffinba ágyaztuk és 5 μm vastagságú metszetekre vágtuk. A metszeteket standard technikákkal deparaffinizáltuk, majd 4 °C-on egy éjszakán át primer antitestekkel, majd 30 percig szobahőmérsékleten másodlagos antitestekkel inkubáltuk. A tárgylemezeket Slowfade Mountant + DAPI-val (Life Technologies, S36964) montíroztuk és lezártuk.

smFISH. RNAScope Multiplex Flourescent Kit-et (Advanced Cell Diagnostics) használtunk a gyártó ajánlásai szerint, a következő módosításokkal. A target retrieval forralási időt 12 percre állítottuk be, a proteáz IV-vel 40 °C-on történő inkubálást pedig 8 percre. A tárgylemezeket Slowfade Mountant+DAPI-val (Life Technologies, S36964) montíroztuk és lezártuk.

Kombinált immunfluoreszcencia és smFISH. Ezt úgy valósítottuk meg, hogy először smFISH-t végeztünk a fent leírtak szerint, a következő változtatásokkal. Az Amp 4 után a szövetszelvényeket mosópufferrel mostuk, primer antitestekkel inkubáltuk egy éjszakán át 4 °C-on, háromszor mostuk 1× TBST-ben, majd másodlagos antitestekkel inkubáltuk 30 percig szobahőmérsékleten. A tárgylemezeket Slowfade Mountant + DAPI-val (Life Technologies, S36964) montíroztuk és lezártuk.

Képelemzés

A szöveti metszetek képeit konfokális mikroszkóp Fluorview FV1200 készülékkel készítettük, a zaj és a jelek átfedésének csökkentése érdekében Kalman és szekvenciális lézeremisszió alkalmazásával. A méretarányos sávokat az FV10-ASW 3.1 Viewer konfokális szoftver segítségével adtuk hozzá az egyes képekhez. A képeket az Image J szoftver45 segítségével fedtük egymásra és vizualizáltuk.

Az antitestek és szondák

Bélorganoid kultúrák

A kripták izolálását követően az egysejtű szuszpenziót 1 μM Jagged-1 peptiddel (Ana-Spec) ellátott Matrigelben (BD Bioscience) reszuszpendáltuk. Nagyjából 300, 25 μl Matrigelbe ágyazott kriptát ültettünk egy 24 lyukú lemez minden egyes lyukába. Miután megszilárdult, a Matrigelt 600 μl táptalajban (Advanced DMEM/F12, Invitrogen) inkubáltuk sztreptomicinnel/penicillinnel és glutamataxszal, valamint EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), Y-276432 dihidroklorid-monohidrát (10 μM, Tochris), N-acetil-1-cisztein (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) és Wnt3A (25 ng ml-1, R&D Systems). A 3. napon friss médiumot cseréltünk, és az organoidokat TrypLE-vel történő disszociációval passziváltuk, majd a 6. napon 1:3 osztási arányban új Matrigelbe szuszpendáltuk. A kiválasztott kísérletekhez az organoidokat RANKL-lal (100 ng ml-1, Biolegends) is kezeltük. A kezelt organoidokat disszociáltuk és mindkét módszerrel scRNS-seq-nek vetettük alá.

Kvantitatív PCR

A teljes hosszúságon alapuló scRNS-seq lemezekről származó tuft-1, tuft-2 és random EpCam+ 16 teljes transzkriptom-amplifikált egyedi sejtjének cDNS-ét használtuk a relatív qPCR-hez. A génexpressziót kvantitatív valós idejű PCR-rel elemeztük LightCycler 480 Instrument II (Roche) készülékkel, LightCycler 480 SYBR green mix (Roche) használatával, a következő primer készletekkel: HPRT1-F, GTTAAGCAGCAGTACAGCCCCCCAAA; HPRT1-R, AGGGCATATATCCAACAACAAACAAACTT; UBC-F, CAGCCGTATATCTTCTCCCCCAGACT; UBC-R, CTCAGAGAGGGATGCCAGTAATCTA; tslp-F, TACTCTCTCAATATCCTATCCCCCTGGCTG; Tlsp-R, CCATTTCCTCTGAGTACCGTCATTTC; Alpi-F, TCCTACACCTCCTCCATTCTCTCTATGG, Alpi-R, CCGCCTGCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCATCTACACCATC; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. A qPCR-primereket minden transzkriptumban exon-exon határra terveztük.

Computációs elemzés

A csepp alapú scRNS-seq adatok előfeldolgozása. A de-multiplexelést, az mm10 transzkriptomhoz való igazítást és az egyedi molekuláris azonosító (UMI) kollapszálást a 10X Genomics által biztosított Cellranger eszközkészlet (1.0.1 verzió) segítségével végeztük. Minden egyes sejt esetében számszerűsítettük azon gének számát, amelyekhez legalább egy leolvasást leképeztünk, majd kizártuk azokat a sejteket, amelyekben 800-nál kevesebb detektált gén volt. A j sejtben az i génre vonatkozó Ei,j expressziós értékeket úgy számoltuk ki, hogy az i génre vonatkozó UMI-számokat elosztottuk a j sejtben lévő UMI-számok összegével, hogy normalizáljuk a lefedettségi különbségeket, majd megszoroztuk 10 000-rel, hogy TPM-szerű értékeket kapjunk, végül pedig kiszámítottuk a log2(TPM + 1) értéket. A kötegelt korrekciót az sva47 R-csomagban implementált ComBat46 segítségével végeztük, az alapértelmezett paraméteres korrekciós módot használva. A kimenet egy korrigált expressziós mátrix volt, amelyet a további elemzés bemeneteként használtunk.

A változó gének kiválasztása úgy történt, hogy egy általánosított lineáris modellt illesztettünk a variációs együttható négyzete és az átlagos expressziós szint közötti kapcsolatra logaritmikus térben, és kiválasztottuk azokat a géneket, amelyek jelentősen (P < 0,05) eltértek az illesztett görbétől48.

A SMART-Seq2 scRNA-seq adatok előfeldolgozása. A BAM-fájlokat az Illumina által biztosított Bcl2Fastq szoftvercsomag v2.17.1.14 segítségével egyesített, de-multiplexált FASTQ-kká alakítottuk át. A párosított végű leolvasásokat a UCSC mm10 egér transzkriptomra térképeztük le a Bowtie49 segítségével, a “-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6” paraméterekkel, amely lehetővé teszi az egy eltérést tartalmazó szekvenciák összehangolását. A gének expressziós szintjét az RSEM50 v1.2.3 által párosított végű módban számított TPM-értékek segítségével számszerűsítettük. Minden egyes sejt esetében számszerűsítettük azoknak a géneknek a számát, amelyekhez legalább egy leolvasás leképeződött, majd kizártuk azokat a sejteket, amelyekben vagy kevesebb mint 3000 detektált gén volt, vagy a transzkriptom leképeződése kevesebb mint 40%-os volt. Ezután a fent leírtak szerint azonosítottuk a nagymértékben változó géneket.

Dimenziócsökkentés PCA és t-SNE segítségével. Az expressziós mátrixot a fent említett változó gének és jó minőségű sejtek alcsoportjaira korlátoztuk, majd az értékeket központosítottuk és skáláztuk, mielőtt bevittük volna őket a főkomponens-elemzésbe (PCA), amelyet a SMART-seq2-adatkészlethez a stats csomagból származó prcomp R-funkcióval hajtottunk végre. A cseppalapú adathalmaz esetében a PCA randomizált közelítését alkalmaztuk, amelyet az rsvd R csomag rpca függvényével valósítottunk meg, a k paramétert 100-ra állítva. Ezt az alacsony rangú közelítést azért használtuk, mert nagyon széles mátrixok esetén több nagyságrenddel gyorsabb a számítása. Tekintettel arra, hogy sok főkomponens a variancia nagyon kis részét magyarázza, a jel-zaj arány jelentősen javítható az n “jelentős” főkomponensből álló részhalmaz kiválasztásával. A PCA után a szignifikáns főkomponenseket a permutációs teszt51 segítségével azonosítottuk, amelyet a jackstraw R csomag permutationPA függvényével hajtottunk végre. Ez a teszt 13, illetve 15 szignifikáns főkomponenst azonosított az 1b. ábra 10X és a SMART-Seq2 adatkészletében, illetve a 2a. ábra kiterjesztett adataiban. Csak ezeknek a szignifikáns főkomponenseknek a pontszámait használtuk a további elemzéshez.

A vizualizáció érdekében az adatkészletek dimenzióját tovább csökkentettük a t-SNE52,53 “Barnes-hut” közelítő változatának segítségével. Ezt az Rtsne R csomagból származó Rtsne függvény segítségével hajtottuk végre, 20 000 iterációt és egy perplexitási beállítást használva, amely az adathalmaz méretétől függően 10 és 30 között változott.

A sejtdifferenciálódási pályák azonosítása diffúziós térképek segítségével

A diffúziós térkép dimenziócsökkentésének lefuttatása előtt az adatokban erősen változó géneket választottunk ki az alábbiak szerint. Először egy nullmodellt illesztettünk az alapszintű sejt-sejt génexpressziós variabilitásra az adatokban, a variációs együttható és az összes expresszált gén UMI-számának átlaga közötti power-law kapcsolatot használva, hasonlóan a korábbi munkákhoz54. Ezután minden gén esetében kiszámítottuk a megfigyelt variációs együttható értéke és a nullmodell által várt érték közötti különbséget (CVdiff). A CVdiff hisztogramja “kövér” farkat mutatott. Kiszámítottuk ennek az eloszlásnak az μ átlagát és σ szórását, és kiválasztottuk az összes olyan gént, amelynél a CVdiff > μ + 1,67σ, így 761 gént kaptunk a további elemzéshez.

A dimenzionalitáscsökkentést a diffúziós térképes megközelítéssel22 végeztük el. Röviden, a sejt-sejt átmenet mátrixát egy Gauss-kernel segítségével számoltuk ki, a kernel szélességét az egyes sejtek helyi szomszédságához igazítva55. Ezt a mátrixot normalizálás után Markov-mátrixszá alakítottuk át. A mátrix vi (i = 0, 1, 2, …) jobb oldali sajátvektorait kiszámítottuk és a λi (i = 0, 1, 2, …) csökkenő sajátérték szerinti sorrendbe rendeztük, miután kizártuk a λ0 = 1-nek megfelelő v0 sajátvektort (ami a Markov-mátrix normalizálási kényszerét tükrözi). A fennmaradó vi sajátvektorok (i = 1, 2, …) meghatározzák a diffúziós térkép beágyazását, és diffúziós komponenseknek (DCk, k = 1, 2, …) nevezzük őket. Észrevettük a λ4 és λ5 közötti spektrális szakadékot, ezért megtartottuk a DC1-DC4 komponenseket mind az eredeti adathalmaz (4. ábra, bővített adatok), mind a különböző bélrégiókból kinyert adatok esetében (2c. ábra).

Szennyező immunsejtek és doubletek eltávolítása

Noha a sejteket a szekvenálás előtt EpCAM segítségével szelektáltuk, a 10X adathalmazban kis számú szennyező immunsejtet figyeltünk meg. Ezeket a 264 sejtet egy első körös felügyelet nélküli klaszterezéssel (a t-SNE térkép sűrűségalapú klaszterezése az fpc R csomagból származó dbscan56 segítségével) eltávolítottuk, mivel rendkívül határozott klasztert alkottak. A SMART-Seq2-adatkészlet esetében több sejt a könyvtár komplexitása szempontjából kiugró értéket képviselt, ami valószínűleg szekvenáló könyvtáranként egynél több egyedi sejtnek felelhetett meg (“doubletek”). Ezeket a sejteket ezután úgy távolítottuk el, hogy kiszámítottuk a sejtenként detektált gének eloszlásának felső 1%-os kvantilisét, és eltávolítottuk az ebben a kvantilisben lévő sejteket.

Klaszterelemzés

Az egyes sejtek expressziójuk szerinti klaszterezéséhez az Infomap gráf-klaszterező algoritmuson9 alapuló, felügyelet nélküli klaszterezési megközelítést használtunk, az egysejtes CyTOF-adatok57 és az scRNA-seq10 esetében alkalmazott megközelítéseket követve. Röviden, egy k legközelebbi szomszédság gráfot építettünk az adatokra, minden egyes sejtpár esetében a szignifikáns főkomponensek pontszámai közötti euklideszi távolságot használva a k legközelebbi szomszédok azonosítására. A k paramétert úgy választottuk meg, hogy összhangban legyen az adathalmaz méretével. Konkrétan a k értéket 200-ra és 80-ra állítottuk be a 7 216 sejtből álló cseppalapú adathalmaz (1b. ábra) és az 1522 sejtből álló SMART-Seq2 adathalmaz (2a. ábra) esetében. A RANKL-kezelt organoidok 5 434 sejtet tartalmaztak, és k értéke 200 volt; a Salmonella és H. polygyrus adathalmaz 9 842 sejtet tartalmazott, és k értéke 500 volt. A sejttípusokon belüli klaszterelemzésekhez, különösen az enteroendokrin és a bojtos sejtek alcsoportjaihoz az euklideszi távolság helyett a Pearson-féle korrelációs távolságot használtuk, és k = 15, k = 30 és k = 40 értéket állítottunk be az enteroendokrin altípusok (533 sejt), valamint a 10X és SMART-Seq2 adatkészletek 166 és 102 bojtos sejtje esetében. A legközelebbi szomszédsági gráfot a cccd R-csomagból származó nng függvény segítségével számoltuk ki. A k-közelebbi szomszédsági gráfot ezután az Infomap9 bemenetként használtuk, amelyet az igraph R csomag infomap.community függvényével valósítottunk meg.

A detektált klasztereket a bélhámsejtek altípusainak ismert markerei segítségével sejttípusokhoz vagy köztes állapotokhoz képeztük le. (Kibővített adatok 1g. ábra, kibővített adatok 2a. ábra). Az enteroendokrin (EEC) sejtek alanalízise (3. ábra) esetében az EEC progenitor klaszterek bármely olyan csoportját, amelynek szignifikáns főkomponens-pontszámok közötti átlagos páronkénti korrelációja r > 0,85 volt, összevonták, így négy klaszter jött létre. Ezt a négy klasztert a Ghrl magas szintje alapján “A” progenitornak, illetve (korai), (középső) vagy (késői) progenitornak (ebben a sorrendben) jelöltük a törzs (Slc12a2, Ascl2, Axin2) és sejtciklus gének csökkenő szintje és az ismert EEC szabályozó faktorok (Neurod1, Neurod2 és Neurog3) növekvő szintje alapján (bővített adatok 5c. ábra, 6. kiegészítő táblázat). A SMART-Seq2 adathalmaz esetében két, magas szintű őssejt-markergéneket expresszáló klasztert (Kiegészített adatok 2a. ábra) összevontunk egy “stem” klaszterré, két másik klasztert pedig egy “TA” klaszterré.

A 4700 sejtből álló tüsző-asszociált epitélium-adatkészlet klaszterelemzésénél a mikrofoltos sejtek nagyon ritkák voltak (0,38%), ezért a ClusterDP módszert58 használtuk az azonosításukra, mivel az empirikusan jobban teljesített, mint a k-nearest-neighbour gráf algoritmus ezen az adatkészleten. A k-nearest-neighbour módszerekhez hasonlóan a ClusterDP-t is szignifikáns (P < 0,05) főkomponens-pontszámok (ebben az esetben 19) felhasználásával futtattuk, és a densityClust R csomagból származó findClusters és densityClust függvények segítségével valósítottuk meg, rho = 1 paraméterekkel.1 és delta = 0,25.

A ritka sejttípusok kivonása a további elemzéshez

A teljes bélállomány (7 216 sejt; 1b. ábra) kezdeti klaszterezése 310 EEC-sejtből és 166 bojtos sejtből álló klasztert mutatott. A tuftsejteket “változatlanul” vettük a részelemzéshez (4a. ábra, b), míg az EGK-sejteket egy második, 239 EGK-sejtből álló klaszterrel kombináltuk, amelyet a regionális adathalmazban azonosítottunk (2a. ábra, jobbra), így összesen 549 EGK-sejtet kaptunk. Egy 16 sejtből álló csoport a Chga és Chgb EEC-markereket a Paneth-sejtek markereivel, köztük a Lyz1, Defa5 és Defa22 markerekkel együttesen fejezte ki, ezért ezeket kettős sejtként értelmeztük és eltávolítottuk az elemzésből, így 533 EEC-sejt maradt, amelyek a 3. ábrán látható elemzés alapját képezték. A proximális és distális vékonybélből származó enterociták expressziós profiljainak összehasonlításához (2b. ábra) a regionális adathalmaz (2a. ábra) 11 665 sejtjéből azonosított 1041 enterocitát használtuk.

A sejttípus-szignatúrák meghatározása

A sejttípusokra maximálisan specifikus gének azonosításához minden lehetséges páronkénti összehasonlítás esetén differenciális expressziós teszteket futtattunk minden egyes klaszterpár között. Ezután egy adott klaszter esetében a feltételezett szignatúra géneket a maximális FDR Q-érték segítségével szűrtük, és a minimális log2(fold változás) alapján rangsoroltuk. A minimális hajtásváltozás és a maximális Q-érték a leggyengébb hatásméretet képviseli az összes páronkénti összehasonlításban; ezért ez egy szigorú kritérium. Az 1c. ábrán, a bővített adatok 2b. ábráján, a bővített adatok 8e. ábráján és a 2-4. és 8. kiegészítő táblázatokban bemutatott sejttípus-szignatúra géneket 0,05 maximális FDR és 0,5 minimális log2(fold change) értékkel kaptuk. A posztmitotikus sejttípus szignatúrák esetében minden gén átlépte ezt a küszöbértéket mind a 3′ (1c. ábra), mind a teljes hosszúságú (kibővített adatok 2b. ábra) adatkészletekben.

A sejttípusokon belüli altípusok szignatúragénjei esetében (3b. ábra, 4b. ábra, kibővített adatok 2b. ábra). 7b. ábra) a feldúsulás kombinált P-értékét (az összes páros tesztre vonatkozóan) a Fisher-módszer segítségével számoltuk ki – ami enyhébb kritérium, mintha egyszerűen a maximális P-értéket vennénk -, és egy 0,01-es maximális FDR Q-értéket használtunk, valamint a minimális log2(fold change) 0,25-ös határértéket a tuftsejt-altípusok esetében (4b. ábra, 7b. ábra, kibővített adatok, 7. kiegészítő táblázat) és 0,1-et az EEC-altípusok esetében (3b. ábra, 6. kiegészítő táblázat). A tuftsejtes szignatúra minden génje megfelelt ezen a cut-off értéken mind a 3′ (4b. ábra), mind a teljes hosszúságú (kibővített adatok, 7b. ábra) adatkészletben, míg az EEC altípusok szignatúráját csak a 3′ alapján határoztuk meg. Az alacsony sejtszám (n = 18) miatt a Fisher-féle kombinált P-értéket használtuk az in vivo mikrofoltos sejtek szignatúrájára is, 0,001-es FDR cut-off értékkel (5d. ábra, 8. kiegészítő táblázat). A markergéneket a minimális log2(fold változás) szerint rangsoroltuk. A differenciális expressziós teszteket a wilcox.test R-funkcióval implementált Mann-Whitney U-teszt (más néven Wilcoxon rank-sum teszt) segítségével végeztük. A fertőzési kísérletekhez (6. ábra) egy kétrészes “akadály” modellt használtunk a technikai minőség és az egerek közötti eltérések ellenőrzésére. Ezt a MAST59 R-csomag segítségével valósítottuk meg, és a differenciális expresszió P-értékeit a likelihood-ratio teszt segítségével számoltuk ki. A többszörös hipotézisvizsgálat korrekcióját az FDR60 ellenőrzésével végeztük az R p.adjust függvény segítségével.

A sejtek pontozása szignált génkészletekkel

Az adott sejtben lévő n gén egy adott készletére vonatkozó pontszám meghatározásához egy “háttér” génkészletet definiáltunk, hogy a szekvenálási lefedettség és a könyvtár komplexitása közötti különbségeket a sejtek között a ref. 12. A háttérgénkészletet úgy választottuk ki, hogy az expressziós szint szempontjából hasonló legyen az érdeklődésre számot tartó génekhez. Konkrétan a 10n legközelebbi szomszédot választottuk ki a kétdimenziós térben, amelyet az átlagos expresszió és a detektálási gyakoriság határoz meg az összes sejtre vonatkozóan. Az adott sejtre vonatkozó szignatúra pontszámot ezután úgy határoztuk meg, hogy az n szignatúra gén átlagos expressziója az adott sejtben, csökkentve a 10n háttérgén átlagos expressziójával az adott sejtben.

A sejttípusok mintavételi gyakoriságának becslése

Minden sejttípus esetében a k méretű mintában legalább n sejt megfigyelésének valószínűségét a negatív binomiális NBcdf(k, n, p) kumulatív eloszlásfüggvényével modelleztük, ahol p az adott sejttípus relatív gyakorisága. Azonos p paraméterű m sejttípus esetén az egyes típusok legalább n alkalommal történő észlelésének teljes valószínűsége NBcdf(k; n, p)m. Ilyen elemzés végezhető a felhasználó által megadott paraméterekkel http://satijalab.org/howmanycells.

EEC dendrogram

Az átlagos expressziós vektorokat mind a 12 EEC alcsoport klaszterére kiszámítottuk, log2(TPM + 1) értékek felhasználásával, és a fent leírtak szerint az EEC alcsoportok között szignifikánsan változónak (P < 0,05) azonosított 1361 génből álló alcsoportra korlátoztuk. Az ezeket a géneket tartalmazó átlagos expressziós vektorokat hierarchikusan klasztereztük a pvclust R-csomag segítségével (Spearman-távolság, ward.D2 klaszterezési módszer), amely minden dendrogramcsomópontra bootstrap megbízhatósági becslést ad empirikus P-értékként 100 000 vizsgálaton keresztül (Extended Data Fig. 6a).

Cellatípus-specifikus transzkripciós faktorok, GPCR-ek és leucinban gazdag ismétlődő fehérjék

Az egerekben transzkripciós faktorokként azonosított összes gén listáját az AnimalTFDB61-ből nyertük. A GPCR-ek készletét az UniProt adatbázisból nyertük (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). Az egyes fehérjék funkcionális annotációit (bővített adatok 2d. ábra) a British Pharmacological Society (BPS) és az International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A) adatbázisából szereztük be. A leucinban gazdag ismétlődő fehérjék listáját a ref. 62. A humán és egér génnevek közötti leképezéshez a humán és egér ortológokat az Ensembl-ből (legújabb, 86. kiadás; http://www.ensembl.org/biomart/martview), a humán és egér génszinonimákat pedig az NCBI-ból (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/) töltöttük le. Minden egyes humán leucinban gazdag ismétlődő gén esetében az összes humán szinonimát leképeztük az egér ortológ génre az ortológ lista segítségével, az egér génneveket pedig az egysejtes adatokban szereplő génekre a szinonimák listája segítségével.

A sejttípusban dúsított transzkripciós faktorokat, GPCR-eket és leucinban gazdag ismétlődő fehérjéket ezután az egyes sejttípusokban dúsított gének listájának és a fent meghatározott transzkripciós faktorok, GPCR-ek és leucinban gazdag ismétlődő fehérjék listájának keresztezésével azonosítottuk. A sejttípusokkal dúsított géneket a SMART-Seq2 adatkészlet felhasználásával úgy határoztuk meg, mint amelyeknél a log2(fold change) minimum 0 és maximum 0,5 FDR volt, megtartva sejttípusonként maximum 10 gént a 2e, f kiterjesztett adatok ábrán (a teljes listákat az 5. kiegészítő táblázat tartalmazza). Ezenkívül egy enyhébb küszöbérték kiválasztásával a sejttípus-specifikus GPCR-ek kiterjedtebb paneljét azonosítottuk (Extended Data 2d. ábra). Ezt úgy értük el, hogy az előző szakaszban leírt páros összehasonlítások helyett minden egyes sejttípust összehasonlítottunk az összes többi sejttel, és kiválasztottuk az összes olyan GPCR-gént, amely differenciálisan kifejeződött (FDR < 0,001).

A sejttípus-arányok változásának vizsgálata

Az egyes elemzett egerekben az egyes sejttípusok kimutatott számát Poisson-folyamat segítségével véletlenszerű számváltozóként modelleztük. A detektálási arányt ezután úgy modellezzük, hogy az adott egérben profilozott sejtek teljes számát adjuk meg offsetváltozóként, az egyes egerek állapotát (kezelés vagy kontroll) pedig kovariátorként adjuk meg. A modellt a stats csomagból származó glm R-paranccsal illesztettük. A kezelés által kiváltott hatás szignifikanciájának P-értékét a regressziós együtthatóra vonatkozó Wald-teszt segítségével értékeltük.

Az EGK-alcsoportok térbeli eloszlásainak szignifikanciájának értékeléséhez (3e. ábra) az összehasonlítás két csoportnál több csoportot érintett. Különösen az volt a nullhipotézisünk, hogy a három bélrégióban (duodenum, jejunum és ileum) kimutatott egyes EEC-alcsoportok aránya egyenlő. E hipotézis tesztelésére varianciaanalízist (ANOVA) alkalmaztunk a fent leírt Poisson-modell illesztésen végzett χ2-teszttel, amelyet a stats csomag anova függvényével valósítottunk meg.

Génkészlet-gazdagodás és génontológiai elemzés

A génontológiai elemzést a goseq R csomag63 segítségével végeztük el, a szignifikánsan differenciálisan kifejezett gének (FDR < 0.05) mint célgének, és háttérként a legalább tíz sejtben log2(TPM + 1) > 3 értékkel kifejezett összes gént.

Adatok elérhetősége

Kód elérhetősége