## Analyse von scRNA-Seq-Daten mit Seurat ### Seurat Das R-Paket *Seurat* dient der Analyse von single-cell-RNA-Seq-Daten. Es wird ausführlich auf dieser Webseite beschrieben: https://satijalab.org/seurat Sie sog. "Vignetten" demonstieren verschiedene Anwendungen mit konkreten beispiel-Daten. Wir folgen der ersten Vignette, dem "Guided Clustering Tutorial", [hier](https://satijalab.org/seurat/articles/pbmc3k_tutorial.html). ### Schnell-Durchlauf Hier ist ein "Schnell-Durchlauf" durch die Haupt-Schritte des Tutorials: ```r library( tidyverse ) library( Seurat ) Read10X( "/var/data/pbmc_3k/filtered_gene_bc_matrices/hg19/" ) -> count_matrix CreateSeuratObject( count_matrix, min.cells=3 ) -> seu NormalizeData( seu ) -> seu ScaleData( seu ) -> seu FindVariableFeatures( seu ) -> seu RunPCA( seu ) -> seu FindNeighbors( seu, dims = 1:20 ) -> seu FindClusters( seu ) -> seu RunUMAP( seu, dims=1:20 ) -> seu DimPlot( seu, reduction = "umap" ) + coord_fixed() ``` ### Datensatz Die verwendeten Beispieldaten sind der PBMC3K-Datensatz, der von 10X (dem Hersteller von "Chromiun", dem mikrofluidischen DropSeq-Gerät für scRNASeq) zur Vefügung gestellt wird. Es enthält 3000 ("3K") PBMCs (peripheral-blood mononuclear cells). Den Datensatz findet man [hier](https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz). Die Datei muss entpackt werden; man erhält dann ein Verzeichnis mit 3 einzelnen Daten ("barcodes.tsv", "genes.tsv" und "matrix.mtx"). Den Pfad zu diesem Verzeichnis muss man in `Read10X` angeben. (Der im Beispiel-Code oben angegebene Pfad funktioniert auf dem papagei-Server.) ### Erläuterungen zum Schnelldurchlauf Die o.g. Befehle führen dann folgende Schritte aus: - EInlesen der Count-Matrix mit Read10X - Erstellen eines Seurat-Objekt (ein Daten-Objekt, in dessen viele "Schubfächer" die Ergebnisse der nachfolgenden Schritte einsortiert werden) - Normalisieren der Daten: $q = \log( k/s * 10000 + 1 )$, wobei $k$ die Anzahl der Read-Counts für ein Gen in einer Zelle ist und $s$ die Gesamtzahl der Reads pro Zelle. Wie wir es schon zuvor getan haben wird ein kleiner Wert ($10^{-4}$ auf der Skala der Read-Anzahlen) hinzu addiert, um zu vermeiden, dass man den Logarithmus von 0 zieht. - Skalieren der log-normalisierten Daten auf Mittelwert 0 und Standardabweichung 1. Die skalierten Daten werden *nur* für die PCA (s.u.) verwendet - Finden von 1000 Genen, die am stärksten über die Zellen variieren. Diese sind besonders geeignet, um zu erkennen, ob sich zwei Zellen ähneln. - Durchführung einer Hauptkomponentenanalyse (principal components analysis, PCA) für die hoch-variablen Gene. Wir erhalten für jede Zelle einen Vektor von einigen wenigen Zahlen, die komprimiert die nützlichen Informationen aus den hochvariablen-Genen so zusammen fassen, das sich Poisson-Rauschen bestmöglich herausmittelt - Nachbarschafts-Suche: Für jede Zelle werden 20 nächste Nachbarn (nearest neighbors, NN) gesucht, d.h. Zellen, die der jeweiligen Zelle möglichst ähnlich sind. Ähnlichkeit wird anhand der ersten 20 Hauptkomponenten (PCA-Dimensionen) bewertet. - Cluster-Suche: Zellen, deren Nachbarschaften sich stark überlappen, werden zu "Clustern" zusammengefasst. Jede Zelle wird so einem Cluster zugeordnet. - UMAP: Um einen zwei-dimensionalen Übersichtsplot zu erhalten, wird jeder Zelle eine Position in einem 2D-Koordinatensystem zugewiesen, und zwar so, dass Nachbarn (siehe vorheriger Schritt) möglichst nah beieinander liegen, und Zellen, die nicht benachtbarts sind, möglichst nicht nah beieinander liegen. Hierzu gibt es verschiedene Verfahren; wir verwenden die "UMAP"-Methode. - Mit DimPlot wird das Ergebnis der UMAP-Funktion nun angezeigt. Die Zellen werden dabei nach Cluster-Zugehörigkeit eingefärbt. ### Plot Wir erhalten folgenden Plot: ![](https://i.imgur.com/dlVCxLL.png) ### Cluster-Zuordnung Mit dem folgenden Seurat-Befehl können wir die Genexpression in Cluster 3 (dunkelgrün, oben) mit der Genexpression in allen anderen Zellen vergleichen: ```r FindMarkers( seu, ident.1 = 3 ) ``` Wir erhalten folgende Ausgabe ``` p_val avg_log2FC pct.1 pct.2 p_val_adj MS4A1 0.000000e+00 3.3862713 0.847 0.053 0.000000e+00 CD79A 0.000000e+00 4.2811486 0.929 0.042 0.000000e+00 TCL1A 3.086810e-274 3.5971860 0.616 0.022 4.233251e-270 CD79B 5.802988e-271 3.4412214 0.903 0.143 7.958218e-267 LINC00926 1.138457e-270 2.8242350 0.554 0.010 1.561280e-266 HLA-DQA1 1.171759e-269 3.0605946 0.884 0.117 1.606951e-265 VPREB3 8.133906e-239 2.4219333 0.483 0.007 1.115484e-234 HLA-DQB1 9.686520e-229 3.0601910 0.852 0.147 1.328409e-224 CD74 6.132410e-190 2.9510962 1.000 0.818 8.409987e-186 HLA-DRA 2.246742e-187 2.7667371 1.000 0.490 3.081182e-183 HLA-DPB1 3.488394e-170 2.3645303 0.986 0.443 4.783983e-166 ... ``` Das Gen MS4A1 (ganz oben) wurde also in 84,7% aller Zellen in Cluster 3 detektiert (`pct.1`) aber nur in 5,3% der anderen Zellen (`pct.2`). Die durchnicttliche Expression ist im Cluster um $2^{3.386}$ als anderswo (`avg_log2FC`, average log2 fold change). Das ist hoch-signifikant (`p_val` is extrem klein). Eine kurzes Googlen zeigt, das MS4A1 ein Protein ist, das sich nur auf der Oberfläche von B-Zellen findet. Cluster 3 enthält also B-Zellen. Wir können auch alle Zellen in dem Block links unten mit dem Rest vergleichen: ```r FindMarkers( seu, ident.1 = c( 0, 1, 4, 6 ) ) ``` Hier finden wir CD3D, die Delta-Kette des T-Zell-Korezeptors. Wo genau wird dieses Gen exprimiert? ```r FeaturePlot( seu, "CD3D" ) ``` ![](https://i.imgur.com/0aphwN8.png) Offebnsichtlich nicht überall in dem "Blob" links unten, sondern nur in Clustern 0, 1, und 4, aber nicht in 6. WIr vergleichen also die Zellen in Cluster 0, 1, 4 mit den Zellen in Cluster 6: ```r FindMarkers( seu, ident.1 = c( 0, 1, 4 ), ident.2 = 6 ) ``` ### Hausaufgabe Was kommt da heraus? Versuchen Sie es selbst, als Hausaufgabe. Versuchen Sie auch, herauszufinden, ob es Unterschiede innerhalb der Cluster-Gruppe 0, 1, und 4 gibt. Wie interpretieren Sie sie? Und was ist mit den Zellen rechts unten im Bild?