celem proponowanej metody jest identyfikacja typów komórek obecnych w mieszaninie pojedynczych komórek. Wkład metody jest matryca ekspresji genów pojedynczych komórek (Mgene×cell), w której wiersze reprezentują geny, a kolumny reprezentują komórki. Poniżej podajemy więcej szczegółów na temat danych wejściowych i różnych etapów proponowanych ram. Ogólne podejście przedstawiono na Rys. 1.
ogólny obieg pracy proponowanej metody. Biorąc pod uwagę matrycę ekspresji genów jednokomórkowych, moduł (a) eliminuje geny, które nie ulegają ekspresji w żadnej komórce. Korzystając z macierzy wynikowej, moduł (B) oblicza odległość euklidesową między komórkami. Wyjście tego modułu jest macierzą odległości, w której wiersze i kolumny są komórkami (dcell×cell). Moduł (C) zmniejsza wymiarowość macierzy odległości za pomocą techniki t-rozproszonego stochastycznego osadzania sąsiada (t-SNE). W tym module stosuje się metodę średniej sylwetki, aby wybrać optymalną liczbę klastrów K. wreszcie w module (D), matryca odległości o niższym wymiarze i optymalna liczba klastrów K uzyskana z modułu (C) są używane jako dane wejściowe do identyfikacji najbardziej stabilnego klastrowania komórek. Rysunek 2 przedstawia szczegóły modułu D.
źródło danych
osiem publicznie dostępnych zestawów danych scRNA-seq oraz pięć zestawów danych symulacyjnych, które wykorzystaliśmy w naszej analizie, znajdują się w materiałach uzupełniających. Wśród ośmiu rzeczywistych zestawów danych, wszystkie z wyjątkiem trzech (Klein51, Patel52, Treutlein53) są uważane za „złoty standard”, ponieważ etykiety komórek są znane w sposób definitywny. Patel52 i Treutlein53 są określane jako „silver standard” przez Kiselev et al.28 ponieważ ich etykiety komórkowe są określane w oparciu o metody obliczeniowe i wiedzę autorów na temat podstawowej biologii.
pozyskaliśmy przetworzone dane ze strony internetowej Hemberg lab (https://hemberg-lab.github.io/scRNA.seq.datasets). Hemberg i in.54 do przechowywania danych wykorzystaj jednoskładnikowy Biokonduktor S4 class55, a pakiet scater56 do kontroli jakości i tworzenia wykresów. Znormalizowane dane są zdeponowane jako pojedynczy obiekt (.Plik RData) oraz informacja o typie komórki jest dostępna w kolumnie cell_type1 gniazda „colData” tego obiektu. Wartości ekspresji genów komórek są zorganizowane jako matryca, w której wiersze są komórkami, a kolumny są genami. W naszej analizie usuwane są geny (cechy), które nie ulegają ekspresji w żadnych komórkach. Nie filtrowaliśmy żadnej komórki w tej analizie.
filtrowanie genów
jak pokazano na Fig. 1A, usuwamy geny / transkrypty, które nie są wyrażone w żadnej komórce (wartość ekspresji wynosi zero we wszystkich komórkach). Takie geny nie mogą dostarczyć użytecznych informacji, które mogą odróżniać typy komórkowe57. Wynik wykonania metody filtrowania na matrycy ekspresji genów jednokomórkowych (Mgene×cell) jest wykorzystywany jako wejście do drugiego modułu proponowanego frameworka.
pomiar różnicy między komórkami
odległość między komórkami oblicza się za pomocą metryki euklidesowej (rys. 1B). Wyjście ten krok być odległość (odmienność) matryca Dcell×komórka. Zmniejszamy Wymiar D, wykonując technikę t-distributed Stochastic neighbor embedding (T-SNE)34,58, nieliniową redukcję wymiarowości/wizualizację (rys. 1C). Będziemy odnosić się do wyjścia jako D ’ cell×l, gdzie 2 ≤ l ≤ cell. W tym badaniu liczba wymiarów wynosi 2.
grupowanie
Identyfikacja optymalnej liczby klastrów
w tej sekcji opisano trzeci moduł proponowanej metody (rys. 1C). W tej analizie T-SNE jest wielokrotnie (N = 50) nakładany na matrycę odległości dcell×cell w celu uzyskania matrycy odległości d’cell×l o zmniejszonej wymiarowości. za każdym razem optymalną liczbę klastrów oblicza się na podstawie metody średniej sylwetki przy użyciu matrycy odległości d’o zmniejszonej wymiarowości. Aby znaleźć optymalną liczbę klastrów k, klastry k-means są stosowane na macierzy D ’ przy użyciu wartości zakresu (domyślnie = 2:20) i wybiera się k, które maksymalizuje średnią miarę sylwetki. Wreszcie, średnia wybranych liczb K w różnych powtórzeniach (N = 50) (zaokrąglona do najbliższej liczby całkowitej) jest uważana za ostateczną optymalną liczbę klastrów.
sylwetka ocenia jakość tego klastrowania w oparciu o to, jak dobrze jego punkty danych są klastrowane. Miara sylwetki jest przypisana do każdego punktu danych, reprezentując, jak blisko punktu danych jest do własnego klastra w porównaniu do innych klastrów. Dla każdego punktu danych i miara ta jest obliczana w następujący sposób:
gdzie A(i) jest średnią odległością między punktem danych i A wszystkimi innymi punktami danych w tym samym klastrze. b (i) jest najmniejszą średnią odległością i do wszystkich punktów w jakimkolwiek innym skupisku, którego i nie jest członkiem. s (i) przyjmuje wartości od -1 do 1, gdzie wysoki wynik dodatni pokazuje, że dany punkt danych jest dobrze zgrupowany (blisko innych punktów we własnym klastrze i daleko od punktów w innych klastrach). Z drugiej strony, wysoki wynik ujemny pokazuje, że punkt danych jest słabo zgrupowany.
K-oznacza grupowanie oparte na metodzie resamplingu
w tej sekcji opisano szczegółowo ostatni moduł proponowanej metody. Jak pokazano na Fig. 2, korzystając z matrycy odległości zmniejszonej wymiarowo D ’ i wybranej liczby klastrów k z poprzedniego kroku, identyfikujemy najbardziej stabilne klastry, generując różne rozwiązania klastrowe (clusteringi (i∈)) i mierzymy stabilność każdego rozwiązania klastrowego w oparciu o metodę resamplingu. Miara stabilności przypisana do każdego konkretnego klastrowania (oznaczona jako clusteringi) określa, jak często klastry K należące do tego klastra są zachowywane, gdy dane wejściowe (D’) są wielokrotnie próbkowane. Resamplowane zbiory danych są generowane z D ’ przez losowe zastąpienie 5% punktów danych (komórek) szumem. Te hałaśliwe zbiory danych są następnie używane jako wejście do algorytmu k-means. W związku z tym,kilka klastrów (clusteringi, j, j ∈ ) są generowane z resamplowanych danych (resamplowane wersje clusteringi).
W tej analizie, biorąc pod uwagę macierz odległości d ’ cell×l oraz optymalną liczbę klastrów k, obliczamy N różnych klastrów (clustering1, …, clusteringn) za pomocą algorytmu klastrowania K-oznacza. Następnie stabilność każdego klastrowania oceniana jest na podstawie metody resamplingu (szara ramka). Wynik stabilności jest przypisywany do każdego klastra na podstawie tego, jak często jego klastry są odzyskiwane, gdy dane wejściowe są zakłócane (ponownie próbkowane). Klastry z maksymalnym wynikiem stabilności są wybierane jako ostateczne rozwiązanie.
aby ocenić stabilność każdego klastra c w klastrze, klaster C jest porównywany ze wszystkimi klastrami w klastrze,które są uzyskiwane z danych resamplingu (clusteringi, j) na podstawie odległości Jaccard. Współczynnik Jaccard59, miara podobieństwa między zbiorami, służy do obliczenia podobieństwa między dwoma zbiorami w następujący sposób:
gdzie wyrażenie A i B są dwoma klastrami, składającymi się z niektórych punktów danych w X = {X1, …, xN}.
Jeśli podobieństwo Jaccard pomiędzy klastrem c (z pierwotnego klastra) a najbardziej podobnym klastrem w klastrze resamplowanym jest równe lub większe niż 0,75, klaster ten jest uważany za stabilny (zachowany). W ten sposób stabilność każdego klastra w klastrach jest obliczana jako procent czasów zachowania klastra (Współczynnik Jaccarda ≥ 0.75) przez m różnych resamplingów.
następnie uśredniamy miary stabilności klastrów K należących do clusteringi i uważamy je za ogólną miarę stabilności clusteringi. Spośród n różnych rozwiązań klastrowania (clusteringi (i∈)) wybieramy rozwiązanie klastrowania z miarą maksymalnej stabilności jako ostateczne rozwiązanie klastrowania.
Rysunek 3 pokazuje szczegóły metody resamplingu, którą wykonaliśmy w celu obliczenia miary stabilności dla każdego klastrowania. Klastry uzyskane przez zastosowanie średniej k na próbkowanym zbiorze danych są porównywane z klastrami z oryginalnych danych wejściowych wyłącznie w oparciu o punkty nieszumowe (punkty danych szumu są wykluczone, gdy dwa klastry są porównywane na podstawie metryki podobieństwa Jaccard.
struktura resamplingu do obliczania miary stabilności dla każdego klastra. Wejście zawiera N punktów danych X = {X1, …, xN}, liczba klastrów k, Liczba resamplingów m i klastrowanie C, które uzyskuje się przez zastosowanie k-means na X. Ta analiza generuje dane m resamplingu przez przypadkowe zastąpienie 5% punktów danych szumem i oblicza m resamplowane klastry na podstawie k-means klastrowania. Każdy klaster c w C jest porównywany z najbardziej podobnym klastrem w klastrze resamplingu i obliczany jest współczynnik Jaccard między dwoma klastrami, podczas gdy punkty szumu są wykluczone. Procent czasu, w którym współczynniki Jaccard są większe od 0.75 jest uważany za miarę stabilności dla klastra c. średnia miar stabilności dla wszystkich klastrów należących do klastra C jest obliczana i uważana za ogólną miarę stabilności dla klastra C.
metody walidacji
używamy 13 różnych zbiorów danych, w których znane są typy komórek (etykiety). Aby zmierzyć poziom podobieństwa między etykietami referencyjnymi i wnioskowanymi etykietami, które są uzyskiwane za pomocą każdej metody klastrowania, używamy trzech różnych wskaźników: skorygowany wskaźnik rand (ARI), skorygowana informacja wzajemna (AMI) i miara V, jak wyjaśniono poniżej.
skorygowany indeks rand
biorąc pod uwagę etykiety komórek, skorygowany indeks Rand (ARI)47 służy do oceny podobieństwa między wnioskowanym klastrowaniem a prawdziwym klastrowaniem. ARI waha się od 0, dla słabego dopasowania (losowe grupowanie), do 1 dla idealnego porozumienia z prawdziwym grupowaniem. Dla zestawu n punktów danych tabela awaryjności jest skonstruowana na podstawie wspólnej liczby punktów danych między dwoma klastrami. Załóżmy, Że X = {X1, X2,…, XR} i Y = {Y1, Y2, …, YC} reprezentują dwa różne klastry z klastrami R I c, odpowiednio. Nakładanie się X i Y można podsumować jako tabelę liczbową MR×C = , gdzie i = 1…R, j = 1…C. Xi i Yj oznaczają klaster w klastrach X i Y, a i I j odnoszą się odpowiednio do numeru wiersza i numeru kolumny tabeli awaryjnej. ARI definiuje się następująco:
skorygowane wzajemne informacje
H(X) jest nieujemne i przyjmuje wartość 0 tylko wtedy, gdy nie ma niepewności określającej przynależność punktu danych do klastra (jest tylko jeden gromada). Wzajemną informację (MI) między dwoma klastrami X i Y oblicza się w następujący sposób:
gdzie P(I, J) oznacza prawdopodobieństwo,że punkt danych należy zarówno do klastra XI w X, jak i do klastra yj w y:
mi jest nieujemną wielkością górną ograniczoną przez entropie H(x) I H(y). Określa on ilościowo informacje dzielone przez oba klastry i dlatego może być uważany za miarę podobieństwa klastrów. Skorygowana miara wzajemnej informacji jest zdefiniowana w następujący sposób:
gdzie oczekiwana wzajemna informacja pomiędzy dwoma losowymi klastrami jest:
gdzie AI i bj są sumami częściowymi tabeli awaryjnej: \({a}_{i}={\sum }_{J\mathrm{=1}}^{C}{n}_{IJ}\) i \({B}_{J}={\sum }_{i\mathrm{=1}}^{R}{n}_{IJ}\).
skorygowana informacja wzajemna (AMI) przyjmuje wartość 1, gdy dwa klastry są identyczne i 0, gdy MI między dwiema partycjami jest równa wartości oczekiwanej z powodu samego przypadku.
miara V
miara V50 jest średnią harmoniczną między dwoma miarami: jednorodnością i kompletnością. Kryteria jednorodności są spełnione, jeśli klastrowanie przypisuje tylko te punkty danych, które są członkami jednej klasy (true cluster) do jednego klastra. Tak więc rozkład klas w każdym klastrze powinien być przekrzywiony do jednej klasy (Entropia zerowa). Aby określić, jak blisko danego klastra jest do tego ideału, Entropia warunkowa rozkładu klasowego danego klastra jest obliczana jako H (C / K), gdzie C = {C1, C2,…, Cl} jest zbiorem klas, A K jest zbiorem K = {K1, K2, …, Km}. W idealnie jednorodnym przypadku wartość ta wynosi 0. Jednak wartość ta zależy od wielkości zbioru danych i rozkładu wielkości klas. Tak więc, ta Entropia warunkowa jest znormalizowana przez maksymalną redukcję entropii, jaką może dostarczyć informacja grupująca, H(C). W związku z tym jednorodność definiuje się następująco:
kompletność jest symetryczna do jednorodności50. W celu spełnienia kryteriów kompletności, klastry muszą przypisać wszystkie te punkty danych, które są członkami jednej klasy do jednego klastra. Aby zmierzyć kompletność, ocenia się rozkład przypisań klastrów w każdej klasie. W doskonale kompletnym rozwiązaniu klastrowania, każdy z tych rozkładów będzie całkowicie przekrzywiony do jednego klastra.
biorąc pod uwagę jednorodność h i kompletność C, miara V jest obliczana jako ważona średnia harmoniczna jednorodności i kompletności:
jeśli β jest większe niż 1, kompletność jest ważona silniej w obliczeniach. Jeśli β jest mniejsze niż 1, jednorodność jest ważona silniej. Ponieważ obliczenia jednorodności, kompletności i miary V są całkowicie niezależne od liczby klas, liczby klastrów, wielkości zbioru danych i algorytmu klastrowania, środki te mogą być stosowane do oceny dowolnego rozwiązania klastrowania.
