A single-cell survey of the small intestinal epithelium

Mice

Wszystkie prace na myszach przeprowadzono zgodnie z Institutional Animal Care and Use Committees (IACUC) oraz z odpowiednimi wytycznymi w Broad Institute i Massachusetts Institute of Technology, odpowiednio z protokołami 0055-05-15 i 0612-058-18. Do wszystkich eksperymentów myszy zostały losowo przydzielone do grup leczenia po dopasowaniu do płci i wieku 7-10-tygodniowych samic lub samców typu dzikiego C57BL/6J lub Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP), uzyskanych z Laboratorium Jacksona (Bar Harbour) lub myszy Gfi1beGFP/+ (Gfi1b-GFP)43. Myszy trzymano w warunkach wolnych od specyficznych patogenów w obiektach dla zwierząt w Broad Institute, Massachusetts Institute of Technology lub Harvard T. H. Chan School of Public Health.

Zakażenie Salmonella enterica i H. polygyrus. Myszy C57BL/6J (Jackson Laboratory) zakażono 200 larwami trzeciego stadium H. polygyrus lub 108 Salmonella enterica, utrzymywanymi w warunkach wolnych od specyficznych patogenów w Massachusetts General Hospital (Charlestown), z protokołem 2003N000158. H. polygyrus rozmnażano zgodnie z wcześniejszym opisem44. Myszy poddano eutanazji 3 i 10 dni po zakażeniu H. polygyrus. W przypadku Salmonella enterica, myszy zakażono naturalnie opornym na streptomycynę szczepem SL1344 S. Typhimurium (108 komórek), jak opisano wcześniej44, i poddano eutanazji 48 h po zakażeniu.

Dysocjacja komórek i izolacja krypt

Izolacja krypt. Jelito cienkie myszy C57BL/6J wild-type, Lgr5-GFP lub Gfi1b-GFP izolowano i płukano w zimnym PBS. Tkankę otwierano podłużnie i krojono na małe fragmenty o długości około 2 mm. Tkankę inkubowano w 20 mM EDTA-PBS na lodzie przez 90 min, wstrząsając co 30 min. Tkankę następnie energicznie wytrząsano, a supernatant zbierano jako frakcję 1 do nowej probówki stożkowej. Tkankę inkubowano w świeżym EDTA-PBS i co 30 min zbierano nową frakcję. Frakcje zbierano do momentu, gdy supernatant składał się prawie wyłącznie z krypt. Ostatnią frakcję (wzbogaconą o krypty) przemywano dwukrotnie w PBS, odwirowywano przy 300g przez 3 min i dysocjowano za pomocą TrypLE express (Invitrogen) przez 1 min w temp. 37 °C. Zawiesina jednokomórkowa była następnie przepuszczana przez filtr 40 μm i wybarwiana do FACS dla scRNA-seq (poniżej) lub używana do hodowli organoidów. Potwierdziliśmy solidność tej metody testując dodatkowe metody izolacji pojedynczych komórek – albo „całe” (zeskrobywanie nabłonka) albo „wzbogacone w kosmki” (frakcja 1; patrz wyżej) – i stwierdziliśmy, że ze względu na wysoki wskaźnik śmiertelności (przez anoikis) postmitotycznych zróżnicowanych komórek (których głównym składnikiem są dojrzałe enterocyty), zawiesina pojedynczych komórek wzbogacona w krypty wiernie reprezentuje skład typów komórek jelita cienkiego (dane nie pokazane).

Izolacja nabłonka związanego z pęcherzykami. Komórki nabłonkowe z nabłonka związanego z pęcherzykami izolowano poprzez pobranie małych wycinków (0,2-0,5 cm) zawierających łatki Peyera z jelita cienkiego myszy C57Bl/6J lub Gfi1beGFP/+.

Sortowanie komórek

Dla eksperymentów opartych na płytkach full-length scRNA-seq, maszyna FACS (Astrios) została użyta do sortowania pojedynczych komórek do każdej studzienki 96-dołkowej płytki PCR zawierającej 5 μl buforu TCL z 1% 2-merkaptoetanolem. Dla izolacji EpCAM+, komórki barwiono na 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience), EpCAM+ (eBioscience); dla specyficznych komórek nabłonkowych, barwiliśmy również na CD24+/- (eBioscience) i c-Kit+/- (eBioscience). W celu wzbogacenia specyficznych populacji komórek nabłonka jelitowego, komórki zostały wyizolowane z myszy Lgr5-GFP, wybarwione przeciwciałami wymienionymi powyżej i bramkowane na GFP-high (komórki macierzyste), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (linie wydzielnicze) lub GFP-/CD24-/EpCAM+ (komórki nabłonkowe). W celu lepszego odzyskania komórek Paneta, dopuściliśmy wyższe parametry rozproszenia bocznego i rozproszenia do przodu w połączeniu z CD24+/c-Kit+, aby zweryfikować odzyskanie komórek Paneta w komórkach EpCAM+. W celu wyizolowania tuft-2, komórki nabłonkowe od trzech różnych myszy barwiono jak powyżej, ale używając EpCAM+/CD45+ do sortowania 2000 pojedynczych komórek. Użyliśmy łagodnej bramki sortowania, aby zapewnić, że uzyskaliśmy wystarczającą liczbę tych rzadkich komórek tuft-2, co doprowadziło do wyższego wskaźnika zanieczyszczenia komórek T, które usunęliśmy w naszej analizie pojedynczych komórek przy użyciu nienadzorowanego grupowania.

Dla sortowania scRNA-seq o pełnej długości, 96-dołkowa płytka została szczelnie zamknięta za pomocą Microseal F i odwirowana przy 800g przez 1 min. Płytka została natychmiast zamrożona na suchym lodzie i przechowywana w temperaturze -80 °C do czasu przygotowania do oczyszczenia lizatu. Komórki populacji masowej były sortowane do probówki Eppendorfa zawierającej 100 μl roztworu TCL z 1% 2-merkaptoetanolem i przechowywane w temperaturze -80 °C.

Dla scRNA-seq opartego na kropelkach, komórki były sortowane z takimi samymi parametrami jak dla scRNA-seq opartego na płytkach, ale były sortowane do probówki Eppendorfa zawierającej 50 μl 0,4% BSA-PBS i przechowywane w lodzie do czasu przejścia na platformę GemCode single-cell.

ScRNA-seq oparty na płytkach

Komórki pojedyncze. Biblioteki zostały przygotowane przy użyciu zmodyfikowanego protokołu SMART-Seq216. W skrócie, oczyszczono lizat RNA przy użyciu kulek RNAClean XP (Agencourt), a następnie przeprowadzono odwrotną transkrypcję przy użyciu odwrotnej transkryptazy Maxima (Life Technologies) i amplifikację całej transkrypcji (WTA) przy użyciu KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) przez 21 cykli. Produkty WTA były oczyszczane za pomocą kulek Ampure XP (Beckman Coulter), kwantyfikowane za pomocą Qubit dsDNA HS Assay Kit (ThermoFisher) i oceniane za pomocą wysokoczułego chipa DNA (Agilent). Biblioteki RNA-seq zostały skonstruowane z oczyszczonych produktów WTA przy użyciu Nextera XT DNA Library Preperation Kit (Illumina). Na każdej płytce, populacja i kontrole bezkomórkowe były przetwarzane przy użyciu tej samej metody, co dla pojedynczych komórek. Biblioteki były sekwencjonowane na urządzeniu Illumina NextSeq 500.

Próbki zbiorcze. Zbiorcze próbki populacji były przetwarzane przez ekstrakcję RNA za pomocą RNeasy Plus Micro Kit (Qiagen), zgodnie z zaleceniami producenta, a następnie przez zmodyfikowany protokół SMART-Seq2 po oczyszczeniu lizatu, jak opisano powyżej.

Droplet-based scRNA-seq

Pojedyncze komórki były przetwarzane przez GemCode Single Cell Platform przy użyciu GemCode Gel Bead, Chip i Library Kits (10X Genomics, Pleasanton), zgodnie z protokołem producenta. W skrócie, pojedyncze komórki były sortowane do 0,4% BSA-PBS. Do każdego kanału dodawano 6 000 komórek, przy czym średni wskaźnik odzysku wynosił 1 500 komórek. Następnie komórki były rozdzielane na kulki żelowe w emulsji w urządzeniu GemCode, gdzie następowała liza komórek i kodowana barkodami odwrotna transkrypcja RNA, a następnie amplifikacja, ścinanie oraz dołączanie adaptera 5′ i indeksu próbki. Biblioteki były sekwencjonowane na Illumina NextSeq 500.

Immunofluorescencja i smFISH

Immunofluorescencja. Barwienie tkanek jelita cienkiego przeprowadzono zgodnie z wcześniejszym opisem34. W skrócie, tkanki były utrwalane przez 14 h w formalinie, osadzane w parafinie i cięte na odcinki o grubości 5 μm. Sekcje były deparafinowane przy użyciu standardowych technik, inkubowane z przeciwciałami pierwszorzędowymi przez noc w 4 °C, a następnie z przeciwciałami drugorzędowymi w temperaturze pokojowej przez 30 min. Szkiełka były utrwalane za pomocą Slowfade Mountant + DAPI (Life Technologies, S36964) i zamykane.

smFISH. Zestaw RNAScope Multiplex Flourescent Kit (Advanced Cell Diagnostics) był używany zgodnie z zaleceniami producenta z następującymi zmianami. Czas wrzenia w celu odzyskania celu został dostosowany do 12 min, a inkubacja z proteazą IV w temperaturze 40 °C została dostosowana do 8 min. Szkiełka montowano przy użyciu Slowfade Mountant+DAPI (Life Technologies, S36964) i uszczelniano.

Połączona immunofluorescencja i smFISH. W tym celu najpierw wykonywano smFISH, jak opisano powyżej, z następującymi zmianami. Po Amp 4, sekcje tkanek płukano w buforze płuczącym, inkubowano z przeciwciałami pierwotnymi przez noc w 4 °C, płukano w 1× TBST trzykrotnie, a następnie inkubowano z przeciwciałami wtórnymi przez 30 min w temperaturze pokojowej. Szkiełka montowano za pomocą Slowfade Mountant + DAPI (Life Technologies, S36964) i uszczelniano.

Analiza obrazu

Obrazy wycinków tkankowych wykonano za pomocą mikroskopu konfokalnego Fluorview FV1200 z wykorzystaniem Kalmana i sekwencyjnej emisji lasera w celu redukcji szumu i nakładania się sygnałów. Paski skali zostały dodane do każdego obrazu przy użyciu oprogramowania konfokalnego FV10-ASW 3.1 Viewer. Obrazy zostały nałożone na siebie i zwizualizowane przy użyciu oprogramowania Image J45.

Przeciwciała i sondy

Hodowle organoidów jelitowych

Po izolacji krypt, zawiesina jednokomórkowa została ponownie zawieszona w Matrigelu (BD Bioscience) z 1 μM peptydu Jagged-1 (Ana-Spec). Około 300 krypt osadzonych w 25 μl Matrigelu posiano do każdego dołka płytki 24 dołkowej. Po zestaleniu, Matrigel inkubowano w 600 μl podłoża hodowlanego (Advanced DMEM/F12, Invitrogen) ze streptomycyną/penicyliną i glutamataksem oraz uzupełnionego EGF (100 ng ml-1, Peprotech), R-spondin-1 (600 ng ml-1, R&D), noggin (100 ng ml-1, Prepotech), monohydrat dwuhydrochlorku Y-276432 (10 μM, Tochris), N-acetylo-1-cysteina (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) i Wnt3A (25 ng ml-1, R&D Systems). Świeża pożywka była wymieniana w dniu 3, a organoidy były pasażowane przez dysocjację z TrypLE i ponownie zawieszane w nowym Matrigelu w dniu 6 w stosunku 1:3. W wybranych eksperymentach organoidy były dodatkowo traktowane RANKL (100 ng ml-1, Biolegends). Traktowane organoidy zostały zdysocjowane i poddane scRNA-seq przy użyciu obu metod.

Jakościowy PCR

CDNA z 16 amplifikowanych całotranskryptomowo pojedynczych komórek tuft-1, tuft-2 i losowych EpCam+ z płytek scRNA-seq opartych na pełnej długości zostało użyte do względnego qPCR. Ekspresję genów analizowano metodą ilościowego PCR w czasie rzeczywistym na aparacie LightCycler 480 Instrument II (Roche) przy użyciu LightCycler 480 SYBR green mix (Roche) z następującymi zestawami primerów: HPRT1-F, GTTAAGCAGTACAGCCCCAAA; HPRT1-R, AGGGCATCCAACAACAAACTT; UBC-F, CAGCCGTATCTTCCCAGACT; UBC-R, CTCAGAGGATGCCAGTAATCTA; tslp-F, TACTCTCAATCCTATCCCTGGCTG; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC; Alpi-F, TCCTACCTCCATTCTATGG, Alpi-R, CCGCCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCATCTACACCATC; Dclk1-R, CCAGCTTCTTAAAGGCTCGAT. Startery qPCR zostały zaprojektowane dla granicy ekson-ekson we wszystkich transkryptach.

Analiza obliczeniowa

Preprocessing of droplet-based scRNA-seq data. De-multipleksowanie, wyrównanie do transkryptomu mm10 i unikalny identyfikator molekularny (UMI) zostały wykonane przy użyciu zestawu narzędzi Cellranger (wersja 1.0.1) dostarczonego przez 10X Genomics. Dla każdej komórki określiliśmy ilościowo liczbę genów, dla których zmapowano przynajmniej jeden odczyt, a następnie wykluczyliśmy wszystkie komórki z mniej niż 800 wykrytymi genami. Wartości ekspresji Ei,j dla genu i w komórce j zostały obliczone przez podzielenie liczby UMI dla genu i przez sumę liczby UMI w komórce j, aby znormalizować różnice w pokryciu, a następnie pomnożenie przez 10,000, aby uzyskać wartości podobne do TPM, a na koniec obliczenie log2(TPM + 1). Korekcję wsadową przeprowadzono przy użyciu programu ComBat46 zaimplementowanego w pakiecie R sva47, stosując domyślny tryb korekcji parametrycznej. Wyjściem była skorygowana macierz ekspresji, którą wykorzystano jako dane wejściowe do dalszej analizy.

Wybór zmiennych genów przeprowadzono poprzez dopasowanie uogólnionego modelu liniowego do zależności między kwadratem współczynnika zmienności a średnim poziomem ekspresji w przestrzeni logarytmicznej i wybranie genów, które znacząco odbiegały (P < 0.05) od dopasowanej krzywej48.

Preprocessing of SMART-Seq2 scRNA-seq data. Pliki BAM zostały przekonwertowane do połączonych, zdemultipleksowanych FASTQs przy użyciu dostarczonego przez Illumina pakietu oprogramowania Bcl2Fastq v2.17.1.14. Sparowane odczyty mapowano do transkryptomu myszy UCSC mm10 za pomocą Bowtie49 z parametrami ’-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6′, co pozwala na wyrównanie sekwencji z jednym niedopasowaniem. Poziom ekspresji genów określano ilościowo przy użyciu wartości TPM obliczanych przez RSEM50 v1.2.3 w trybie paired-end. Dla każdej komórki określiliśmy ilościowo liczbę genów, dla których zmapowano przynajmniej jeden odczyt, a następnie wykluczyliśmy wszystkie komórki, w których wykryto mniej niż 3000 genów lub mapowanie transkryptomu wynosiło mniej niż 40%. Następnie zidentyfikowaliśmy geny o wysokiej zmienności w sposób opisany powyżej.

Redukcja wymiarowości za pomocą PCA i t-SNE. Ograniczyliśmy macierz ekspresji do podzbiorów zmiennych genów i komórek wysokiej jakości opisanych powyżej, a następnie wyśrodkowaliśmy i przeskalowaliśmy wartości przed wprowadzeniem ich do analizy składowych głównych (PCA), która została zaimplementowana przy użyciu funkcji R prcomp z pakietu stats dla zbioru danych SMART-seq2. Dla zbioru danych opartego na kropelkach użyliśmy randomizowanego przybliżenia PCA, zaimplementowanego przy użyciu funkcji rpca z pakietu rsvd R, z parametrem k ustawionym na 100. To przybliżenie o niskiej randze zostało użyte, ponieważ jest ono o kilka rzędów wielkości szybsze do obliczenia dla bardzo szerokich macierzy. Biorąc pod uwagę, że wiele głównych składowych wyjaśnia bardzo małą część wariancji, stosunek sygnału do szumu może być znacznie poprawiony poprzez wybór podzbioru n „znaczących” głównych składowych. Po PCA zidentyfikowano znaczące składowe główne za pomocą testu permutacji51 , zaimplementowanego przy użyciu funkcji permutationPA z pakietu jackstraw R. Test ten zidentyfikował 13 i 15 znaczących składowych głównych. Test ten zidentyfikował 13 i 15 znaczących składowych głównych odpowiednio w zbiorach danych 10X i SMART-Seq2 z Rys. 1b i Extended Data Rys. 2a. Wyniki tylko z tych istotnych składowych głównych zostały użyte jako dane wejściowe do dalszej analizy.

Dla celów wizualizacji, wymiarowość zbiorów danych została zredukowana przy użyciu przybliżonej wersji t-SNE52,53 „Barnes-hut”. Zostało to zaimplementowane przy użyciu funkcji Rtsne z pakietu Rtsne R przy użyciu 20 000 iteracji i ustawieniu perplexity, które wahało się od 10 do 30 w zależności od wielkości zbioru danych.

Identyfikacja trajektorii różnicowania komórek przy użyciu map dyfuzji

Przed uruchomieniem redukcji wymiarowości mapy dyfuzji wybraliśmy wysoce zmienne geny w danych w następujący sposób. Najpierw dopasowaliśmy model zerowy dla podstawowej zmienności ekspresji genów komórka-komórka w danych, wykorzystując zależność power-law pomiędzy współczynnikiem zmienności a średnią liczbą UMI wszystkich ekspresjonowanych genów, podobnie jak w poprzedniej pracy54. Następnie dla każdego genu obliczono różnicę między wartością obserwowanego współczynnika zmienności a wartością oczekiwaną przez model zerowy (CVdiff). Histogram CVdiff wykazywał „gruby” ogon. Obliczyliśmy średnią μ i odchylenie standardowe σ tego rozkładu i wybraliśmy wszystkie geny, dla których CVdiff > μ + 1.67σ, uzyskując 761 genów do dalszej analizy.

Zredukowaliśmy wymiarowość stosując metodę mapy dyfuzji22. W skrócie, macierz przejścia komórka-komórka została obliczona przy użyciu jądra gaussowskiego, z szerokością jądra dostosowaną do lokalnego sąsiedztwa każdej komórki55. Macierz ta po normalizacji została przekształcona w macierz Markowa. Prawe wektory własne vi (i = 0, 1, 2, …) tej macierzy zostały obliczone i posortowane w kolejności malejącej wartości własnej λi (i = 0, 1, 2, …), po wyłączeniu „górnego” wektora własnego v0, odpowiadającego λ0 = 1 (co odzwierciedla ograniczenie normalizacji macierzy markowskiej). Pozostałe wektory własne vi (i = 1, 2, …) definiują osadzenie mapy dyfuzyjnej i są określane jako składowe dyfuzji (DCk, k = 1, 2, …). Zauważyliśmy lukę spektralną pomiędzy λ4 i λ5, a zatem zachowaliśmy DC1-DC4 zarówno dla początkowego zbioru danych (Extended Data Fig. 4), jak i danych wyodrębnionych z różnych regionów jelita (Fig. 2c).

Usuwanie zanieczyszczających komórek odpornościowych i dubletów

Pomimo, że komórki były sortowane przed sekwencjonowaniem przy użyciu EpCAM, niewielka liczba zanieczyszczających komórek odpornościowych została zaobserwowana w zbiorze danych 10X. Te 264 komórki zostały usunięte we wstępnej rundzie nienadzorowanego grupowania (grupowanie oparte na gęstości mapy t-SNE przy użyciu dbscan56 z pakietu R fpc), ponieważ tworzyły one bardzo wyraźne skupisko. W przypadku zbioru danych SMART-Seq2 kilka komórek było wartościami odstającymi pod względem złożoności biblioteki, co mogło odpowiadać więcej niż jednej pojedynczej komórce na bibliotekę sekwencjonowania („dublety”). Komórki te zostały następnie usunięte przez obliczenie górnego 1% kwantyla dystrybucji genów wykrytych na komórkę i usunięcie wszystkich komórek w tym kwantylu.

Analiza skupień

Aby pogrupować pojedyncze komórki według ich ekspresji, zastosowaliśmy nienadzorowane podejście do grupowania, oparte na algorytmie grupowania grafów Infomap9, wzorując się na podejściach dla danych CyTOF z pojedynczych komórek57 i scRNA-seq10. W skrócie, skonstruowaliśmy graf k-najbliższych sąsiadów na danych wykorzystując, dla każdej pary komórek, odległość euklidesową pomiędzy wynikami znaczących głównych składowych w celu zidentyfikowania k najbliższych sąsiadów. Parametr k został dobrany tak, aby był zgodny z rozmiarem zbioru danych. W szczególności, k zostało ustawione na 200 i 80 odpowiednio dla zbioru danych opartego na kropelkach, składającego się z 7 216 komórek (Rys. 1b) i dla zbioru danych SMART-Seq2, składającego się z 1 522 komórek (Extended Data Rys. 2a). Organoidy poddane działaniu RANKL zawierały 5 434 komórek, a k zostało ustawione na 200; zestaw danych Salmonella i H. polygyrus zawierał 9 842 komórek, a k zostało ustawione na 500. Do analizy skupień w obrębie typów komórek, w szczególności podzbiorów komórek enteroendokrynnych i tuftowych, użyliśmy odległości korelacyjnej Pearsona zamiast odległości euklidesowej i ustawiliśmy k = 15, k = 30 i k = 40 dla podtypów enteroendokrynnych (533 komórki) oraz dla 166 i 102 komórek tuftowych odpowiednio w zbiorach danych 10X i SMART-Seq2. Wykres najbliższego sąsiedztwa został obliczony przy użyciu funkcji nng z pakietu R cccd. Graf k-najbliższych sąsiadów został następnie użyty jako wejście do Infomap9, zaimplementowanej przy użyciu funkcji infomap.community z pakietu igraph R.

Wykryte skupiska zostały zmapowane do typów komórek lub stanów pośrednich przy użyciu znanych markerów dla podtypów komórek nabłonka jelitowego. (Extended Data Fig. 1g, Extended Data Fig. 2a). Dla subanalizy komórek enteroendokrynnych (EEC) (ryc. 3), każda grupa klastrów progenitorowych EEC ze średnimi korelacjami parami pomiędzy znaczącymi wynikami składowych głównych r > 0.85 została połączona, co dało cztery klastry. Oznaczyliśmy te cztery klastry jako progenitor 'A’ na podstawie wysokiego poziomu Ghrl, lub progenitor (wczesny), (środkowy) lub (późny) (w tej kolejności) na podstawie zmniejszających się poziomów genów pnia (Slc12a2, Ascl2, Axin2) i cyklu komórkowego oraz zwiększających się poziomów znanych czynników regulacyjnych EEC (Neurod1, Neurod2 i Neurog3) (Extended Data Fig. 5c, Supplementary Table 6). Dla zbioru danych SMART-Seq2, dwa klastry wyrażające wysoki poziom genów markerów komórek macierzystych (Extended Data Fig. 2a) zostały połączone w klaster „stem”, a dwa inne klastry zostały połączone w klaster „TA”.

Do analizy skupień zestawu danych nabłonka związanego z mieszkiem włosowym, składającego się z 4700 komórek, komórki mikrofold były bardzo rzadkie (0,38%) i dlatego do ich identyfikacji zastosowano metodę ClusterDP58 , ponieważ empirycznie wypadła ona lepiej niż algorytm grafu k-najbliższych sąsiadów na tym zestawie danych. Podobnie jak w przypadku metod k-nearest-neighbour, ClusterDP została uruchomiona przy użyciu znaczących (P < 0.05) wyników składowych głównych (w tym przypadku 19) jako danych wejściowych i została zaimplementowana przy użyciu funkcji findClusters i densityClust z pakietu densityClust R przy użyciu parametrów rho = 1.1 i delta = 0.25.

Wyodrębnienie rzadkich typów komórek do dalszej analizy

Wstępne grupowanie zbioru danych z całego jelita (7,216 komórek; Fig. 1b) wykazało skupisko 310 komórek EEC i 166 komórek kępkowych. Komórki kępki zostały wzięte „jako takie” do analizy cząstkowej (Ryc. 4a, b), podczas gdy komórki EEC zostały połączone z drugim skupiskiem 239 komórek EEC, które zostały zidentyfikowane w regionalnym zbiorze danych (Ryc. 2a, po prawej) dla uzyskania całkowitej liczby 549 komórek EEC. Grupa 16 komórek wykazywała współekspresję markerów EEC Chga i Chgb z markerami komórek Panetha, w tym Lyz1, Defa5 i Defa22, dlatego też zostały one zinterpretowane jako dublety i usunięte z analizy, pozostawiając 533 komórki EEC, które stanowiły podstawę do analizy na rycinie 3. Do porównania profili ekspresji enterocytów z proksymalnego i dystalnego jelita cienkiego (ryc. 2b) użyto 1 041 enterocytów zidentyfikowanych z 11 665 komórek w regionalnym zestawie danych (ryc. 2a).

Definiowanie sygnatur typu komórkowego

Aby zidentyfikować maksymalnie specyficzne geny dla typów komórek, przeprowadziliśmy testy ekspresji różnicowej pomiędzy każdą parą klastrów dla wszystkich możliwych porównań parami. Następnie, dla danego klastra, potencjalne geny sygnaturowe były filtrowane przy użyciu maksymalnej wartości FDR Q i uszeregowane według minimalnej zmiany log2(fold change). Minimalna zmiana fałdu i maksymalna wartość Q reprezentują najsłabszą wielkość efektu we wszystkich porównaniach parami; jest to więc kryterium rygorystyczne. Geny sygnaturowe typu komórkowego przedstawione na Rys. 1c, Extended Data Rys. 2b, Extended Data Rys. 8e oraz w Tabelach 2-4 i 8 zostały uzyskane przy zastosowaniu maksymalnego FDR równego 0,05 i minimalnej log2(fold change) równej 0,5. W przypadku sygnatur typów komórek postmitotycznych, wszystkie geny przekroczyły ten próg zarówno w 3′ (Ryc. 1c), jak i w pełnej długości (Extended Data Fig. 2b).

W przypadku genów sygnatur dla podtypów w obrębie typów komórek (Ryc. 3b, Ryc. 4b, Extended Data Fig. 7b), połączona wartość P (dla wszystkich testów parami) dla wzbogacenia została obliczona przy użyciu metody Fishera – bardziej łagodnego kryterium niż po prostu biorąc maksymalną wartość P – i zastosowano maksymalną wartość FDR Q 0,01, wraz z odcięciem minimalnej wartości log2(fold change) 0,25 dla podtypów komórek tuft (Fig. 4b, Extended Data Fig. 7b, Supplementary Table 7) i 0,1 dla podtypów EEC (Fig. 3b, Supplementary Table 6). Wszystkie geny w sygnaturze komórek kępkowych przeszły ten próg odcięcia zarówno w zestawach danych 3′ (Ryc. 4b), jak i pełnej długości (Extended Data Fig. 7b), podczas gdy sygnatury podtypów EEC zostały zdefiniowane tylko przy użyciu 3′. Ze względu na niską liczbę komórek (n = 18), połączona wartość P Fishera została również użyta dla sygnatury komórek mikroflory in vivo, z punktem odcięcia FDR równym 0,001 (Fig. 5d, Tabela uzupełniająca 8). Geny markerowe uszeregowano według minimalnej log2(fold change). Testy ekspresji różnicowej przeprowadzono przy użyciu testu U Manna-Whitneya (znanego również jako test Wilcoxona rank-sum) zaimplementowanego przy użyciu funkcji R wilcox.test. Do eksperymentów infekcyjnych (Rys. 6) użyliśmy dwuczęściowego modelu „hurdle”, aby kontrolować zarówno jakość techniczną, jak i zmienność między myszami. Model ten został zaimplementowany przy użyciu pakietu R MAST59, a wartości P dla ekspresji różnicowej zostały obliczone przy użyciu testu likelihood-ratio. Korektę testowania wielokrotnych hipotez przeprowadzono poprzez kontrolę FDR60 przy użyciu funkcji R p.adjust.

Punktowanie komórek przy użyciu zestawów genów sygnaturowych

Aby uzyskać wynik dla określonego zestawu n genów w danej komórce, zdefiniowano zestaw genów „tła”, aby kontrolować różnice w pokryciu sekwencjonowania i złożoności biblioteki między komórkami w sposób podobny do ref. 12. Zestaw genów tła został wybrany tak, aby był podobny do genów będących przedmiotem zainteresowania pod względem poziomu ekspresji. W szczególności, wybrano 10n najbliższych sąsiadów w dwuwymiarowej przestrzeni zdefiniowanej przez średnią ekspresję i częstotliwość wykrywania we wszystkich komórkach. Wynik sygnatury dla tej komórki został następnie zdefiniowany jako średnia ekspresja n genów sygnatury w tej komórce, pomniejszona o średnią ekspresję 10n genów tła w tej komórce.

Oszacowanie częstotliwości pobierania próbek typu komórek

Dla każdego typu komórek prawdopodobieństwo zaobserwowania co najmniej n komórek w próbce o rozmiarze k jest modelowane przy użyciu funkcji rozkładu kumulatywnego ujemnego dwumianu NBcdf(k, n, p), gdzie p jest względną obfitością tego typu komórek. Dla m typów komórek z tym samym parametrem p, całkowite prawdopodobieństwo zaobserwowania każdego typu co najmniej n razy wynosi NBcdf(k; n, p)m. Taka analiza może być przeprowadzona z parametrami określonymi przez użytkownika w http://satijalab.org/howmanycells.

Dendrogram EWG

Średnie wektory ekspresji obliczono dla wszystkich 12 skupisk podzbiorów EWG, używając wartości log2(TPM + 1), i ograniczono do podzbioru 1,361 genów zidentyfikowanych jako istotnie zmienne pomiędzy podzbiorami EWG (P < 0.05), jak opisano powyżej. Średnie wektory ekspresji obejmujące te geny zostały hierarchicznie zgrupowane przy użyciu pakietu R pvclust (odległość Spearmana, metoda grupowania ward.D2), który dostarcza bootstrapowe oszacowania ufności dla każdego węzła dendrogramu jako empiryczną wartość P na 100,000 prób (Extended Data Fig. 6a).

Cell-type-specific transcription factors, GPCRs and leucine-rich repeat proteins

Lista wszystkich genów zidentyfikowanych jako działające jako czynniki transkrypcyjne u myszy została uzyskana z AnimalTFDB61. Zestaw GPCRs uzyskano z bazy 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). Adnotacje funkcjonalne dla każdego białka (Extended Data Fig. 2d) uzyskano z British Pharmacological Society (BPS) oraz International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A). Lista białek o powtórzeniach bogatych w leucynę została zaczerpnięta z ref. 62. Aby odwzorować nazwy genów ludzkich na mysie, ortologi ludzkie i mysie zostały pobrane z Ensembl (najnowsze wydanie 86; http://www.ensembl.org/biomart/martview), a synonimy genów ludzkich i mysich z NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/). Dla każdego ludzkiego genu o powtórzeniu bogatym w leucynę wszystkie ludzkie synonimy zostały zmapowane do ortologicznego genu u myszy przy użyciu listy ortologów, a nazwy genów u myszy zostały zmapowane do nazw w danych jednokomórkowych przy użyciu listy synonimów.

Wzbogacone w typ komórki czynniki transkrypcyjne, GPCR i białka o powtórzeniu bogatym w leucynę zostały następnie zidentyfikowane poprzez przecięcie listy genów wzbogaconych w każdym typie komórki z listami czynników transkrypcyjnych, GPCR i białek o powtórzeniu bogatym w leucynę zdefiniowanymi powyżej. Geny wzbogacone w typ komórki zostały zdefiniowane przy użyciu zbioru danych SMART-Seq2 jako te z minimalną log2(fold change) równą 0 i maksymalnym FDR równym 0,5, zachowując maksymalnie 10 genów na typ komórki w Extended Data Fig. 2e, f (pełne listy znajdują się w Supplementary Table 5). Dodatkowo zidentyfikowano bardziej obszerny panel GPCR specyficznych dla danego typu komórki (Extended Data Fig. 2d) poprzez wybór bardziej łagodnego progu. Osiągnięto to poprzez porównanie każdego typu komórek do wszystkich innych komórek, zamiast porównań parami opisanych w poprzedniej sekcji, i wybierając wszystkie geny GPCR, które były różnie wyrażone (FDR < 0,001).

Testowanie zmian w proporcjach typów komórek

Modelujemy wykrytą liczbę każdego typu komórek u każdej analizowanej myszy jako losową zmienną liczoną za pomocą procesu Poissona. Szybkość wykrywania jest następnie modelowana przez dostarczenie całkowitej liczby komórek profilowanych u danej myszy jako zmiennej offsetowej, z warunkiem każdej myszy (leczenie lub kontrola) dostarczonym jako kowariant. Model został dopasowany przy użyciu polecenia R glm z pakietu stats. Wartość P dla istotności efektu wywołanego przez leczenie oceniano za pomocą testu Walda na współczynniku regresji.

Dla oceny istotności rozkładów przestrzennych podzbiorów EWG (ryc. 3e), porównanie obejmowało więcej niż dwie grupy. W szczególności, naszą hipotezą zerową było założenie, że udział każdego podzbioru EEC wykrytego w trzech regionach jelita (dwunastnica, jelito czcze i jelito kręte) jest równy. Aby przetestować tę hipotezę, użyliśmy analizy wariancji (ANOVA) z testem χ2 na opisanym powyżej dopasowaniu modelu Poissona, zaimplementowanym przy użyciu funkcji anova z pakietu stats.

Wzbogacenie zbioru genów i analiza ontologii genów

Analizę ontologii genów przeprowadzono przy użyciu pakietu goseq R63, używając istotnie różnie wyrażonych genów (FDR < 0.05) jako geny docelowe, a wszystkie geny ulegające ekspresji z log2(TPM + 1) > 3 w co najmniej dziesięciu komórkach jako tło.

Dostępność danych

Dostępność kodu

.