Seurat 指導(dǎo)聚類教程

參照官網(wǎng)教程 用了自己的一批真實的數(shù)據(jù)跟压，總共有7038個細胞。以下是cellranger count跑出來的標準結(jié)果。

我們從讀取數(shù)據(jù)開始。Read10X函數(shù)從10X讀取cellranger流程的輸出臭杰，返回UMI計數(shù)矩陣。矩陣中的值表示每個特征(即基因;在每個細胞(列)中檢測到的谚中。

接下來渴杆，我們使用count矩陣創(chuàng)建一個Seurat對象。該對象充當(dāng)一個容器藏杖，其中包含單細胞數(shù)據(jù)集的數(shù)據(jù)(如計數(shù)矩陣)和分析(如PCA或聚類結(jié)果)将塑。

library(Seurat)
library(dplyr)

讀取數(shù)據(jù)：

real_10x_data = Read10X(data.dir = "\\example_G48E2L1\\filtered_feature_bc_matrix")
Read10X(data.dir = NULL, gene.column = 2, unique.features = TRUE)

data.dir 參數(shù)包含矩陣的目錄。包含matrix.mtx蝌麸，gene.tsv（或features.tsv）和barcodes.tsv。為了加載多個數(shù)據(jù)目錄艾疟，可以給出一個向量或命名向量来吩。如果給定了命名向量，則cell barcode 名稱將以該名稱為前綴蔽莱。

gene.column 參數(shù)指定基因在哪一列弟疆。features.tsv或gene.tsv用于基因名稱的tsv;默認是2，表示第二列是基因名盗冷，我們來看一下features.tsv怠苔，包含3列：

unique.features 參數(shù)默認為TRUE，表示 使features name unique仪糖。

如果features.csv表明數(shù)據(jù)具有多個數(shù)據(jù)類型柑司，則返回一個包含每種類型數(shù)據(jù)的稀疏矩陣的列表。否則將返回一個包含表達式數(shù)據(jù)的稀疏矩陣锅劝。

使用原始數(shù)據(jù)(非規(guī)范化數(shù)據(jù))初始化Seurat對象攒驰。

CreateSeuratObject(counts, project = "SeuratProject", assay = "RNA",
  min.cells = 0, min.features = 0, names.field = 1,
  names.delim = "_", meta.data = NULL)

注意:在以前的版本(<3.0)中近迁，該函數(shù)還接受一個參數(shù)來設(shè)置“檢測到的”特征(基因)的表達閾值。為了簡化初始化過程/假設(shè)簸州，刪除了此功能鉴竭。如果您仍然希望為特定的數(shù)據(jù)集設(shè)置這個閾值，那么只需在調(diào)用此函數(shù)之前對輸入表達式矩陣進行篩選即可岸浑。

> G48E2L1 <-CreateSeuratObject(counts = real_10x_data, project = "G48E2L1",)
> G48E2L1
An object of class Seurat
33538 features across 7038 samples within 1 assay
Active assay: RNA (33538 features)

可以發(fā)現(xiàn)7038samples 和網(wǎng)頁版報告一致搏存，33538 features也和features.csv數(shù)量一致。

count matrix 數(shù)據(jù)長什么矢洲？

例如璧眠，計數(shù)矩陣存儲在G48E2L1[["RNA"]]@counts中读虏。

只看幾個基因：

>real_10x_data[c("CD3D", "TCL1A","MS4A1"),1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGGAATGGA’, ‘AAACCTGAGGCTCAGA’, ‘AAACCTGCAGCCTATA’ ... ]]
CD3D  6 . 7 3 3 12 5 1 6 1 . 5 2 5 7 3 4 14 4 . 1 4 2 1 7 2 3 5 13 7
TCL1A . . . . .  . . . . . . . . . . . .  . . . . . . . . . . .  . .
MS4A1 . . . . .  . . . . . . . . . . . .  . . . . . . . . . . .  . .

點. 矩陣中的值表示0(未檢測到分子)责静。由于scRNA-seq矩陣中的大多數(shù)值都是0，所以Seurat盡可能使用稀疏矩陣表示掘譬。這為Drop-seq/inDrop/10x數(shù)據(jù)節(jié)省了大量內(nèi)存和速度泰演。

如果需要查看稀疏矩陣的空間大小（個人理解）葱轩，這些可以忽略睦焕。

> dense.size = object.size(as.matrix(real_10x_data))
> dense.size
1891208224 bytes
> spare.size <- object.size(real_10x_data)
> spare.size
179737888 bytes
> dense.size/spare.size
10.5 bytes

標準的預(yù)處理流程：

以下步驟包含Seurat中scRNA-seq數(shù)據(jù)的標準預(yù)處理工作流藐握。這些代表細胞的選擇和過濾基于QC指標，數(shù)據(jù)歸一化和縮放垃喊，并檢測高度可變的特征猾普。

QC和選擇細胞進行進一步分析
Seurat允許您輕松地探索QC指標，并根據(jù)任何用戶定義的標準過濾單元格本谜。大家通常使用的一些QC指標包括

• 在每個細胞中檢測到的獨特基因的數(shù)量
  * 低質(zhì)量的細胞或空液滴通常只有很少的基因
  * 細胞雙重或多胞胎可能表現(xiàn)出異常高的基因計數(shù)
• 同樣初家，在一個細胞內(nèi)檢測到的分子總數(shù)(與獨特的基因密切相關(guān))
• reads map到線粒體基因組的比例
  * 低質(zhì)量/瀕死細胞常表現(xiàn)出廣泛的線粒體污染
  * 我們使用PercentageFeatureSet函數(shù)計算線粒體QC指標，該函數(shù)計算來自一組特征的計數(shù)的百分比
  * 我們使用所有以MT-開頭的基因作為一組線粒體基因

"[[" 操作符可以向?qū)ο笤獢?shù)據(jù)添加列乌助。新增一列儲存QC統(tǒng)計最好不過了溜在。
G48E2L1[["percent.mt"]] <- PercentageFeatureSet(G48E2L1, pattern = "^MT-")

QC指標存儲在哪里？

> head(G48E2L1@meta.data,5)
                 orig.ident nCount_RNA nFeature_RNA percent.mt
AAACCTGAGGAATGGA    G48E2L1       8276         2922   6.911551
AAACCTGAGGCTCAGA    G48E2L1       2765         1253  10.886076
AAACCTGCAGCCTATA    G48E2L1       8529         2841   8.136945
AAACCTGGTCAGAATA    G48E2L1       9102         3003   5.207647
AAACCTGGTTAAAGAC    G48E2L1       3539         1721  13.139305

在官方的樣例中他托，可視化QC指標掖肋，并進行cell過濾。

• 過濾unique feature counts 大于2500或少于200的細胞
• 過濾線粒體比例大于5%的細胞

先來可視化看一下：

VlnPlot()是Seurat中用于繪制單細胞數(shù)據(jù)的小提琴圖函數(shù)(基因表達赏参、指標志笼、PC分數(shù)等)，小提琴圖用于顯示數(shù)據(jù)分布及其概率密度把篓。

VlnPlot(G48E2L1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

FeatureScatter通常用于可視化 feature-feature 關(guān)系纫溃，也可以用于計算對象的任何東西，i.e. 對象數(shù)據(jù)中的列韧掩，PC分數(shù)等紊浩。 個人理解：就是用點圖看兩個數(shù)據(jù)之間的相關(guān)性。

plot1 <- FeatureScatter(G48E2L1, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(G48E2L1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

image.png

官方教程中在這里過濾掉 2500 > nFeature_RNA >200 和percent.mt < 5的數(shù)據(jù)：

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

但是我不想過濾揍很，本文數(shù)據(jù)沒有做過濾處理郎楼。哈哈哈！

數(shù)據(jù)標準化

從數(shù)據(jù)集中刪除不需要的細胞后窒悔，下一步是數(shù)據(jù)標準化。默認情況下敌买，我們使用全局縮放歸一化方法“LogNormalize”简珠，它將每個細胞的特征表達式測量值歸一化為總的表達式，再乘以一個縮放因子(默認為10,000)虹钮，對結(jié)果對數(shù)化處理聋庵。標準化的數(shù)值存儲在pbmc[["RNA"]]@data中。

G48E2L1 <- NormalizeData(G48E2L1, normalization.method = "LogNormalize", scale.factor = 10000)
G48E2L1[["RNA"]]@data

高度可變特征的識別(feature selection)

接下來芙粱，我們將計算數(shù)據(jù)集中顯示高細胞間差異的特征子集(i祭玉。e，它們在一些細胞中高表達春畔，在另一些細胞中低表達)脱货。我們和其他人發(fā)現(xiàn)岛都，在下游分析中關(guān)注這些基因有助于在單細胞數(shù)據(jù)集中突出生物信號。

這里詳細描述了Seurat3中的過程振峻，并通過直接建模單細胞數(shù)據(jù)中固有的均值-方差關(guān)系改進了先前的版本臼疫，并在FindVariableFeatures函數(shù)中實現(xiàn)。默認情況下扣孟，我們?yōu)槊總€數(shù)據(jù)集返回2,000個特性烫堤。這些將用于下游分析，如PCA凤价。

Find variable features

識別“平均變異性圖”上的異常點鸽斟。

FindVariableFeatures(object, ...)

如何選擇****selection.method:

vst:首先，用局部多項式回歸(loess)擬合對數(shù)(方差)與對數(shù)(均值)的關(guān)系利诺。然后使用觀察到的平均值和期望的方差(由擬合線給出)對特征值進行標準化富蓄。然后，在裁剪到最大值之后立轧，根據(jù)標準化的值計算特征方差(參見clip)格粪。max參數(shù))。

mean.var.plot (mvp):首先氛改，使用一個函數(shù)計算每個特征的平均表達式(mean.function)和離散度(diffusion .function)帐萎。接下來，根據(jù)每個bin的平均表達式將特征劃分為number .bin (默認 20)胜卤，并計算每個bin內(nèi)的離散度z-score疆导。這樣做的目的是識別變量特征，同時控制可變性和平均表達之間的強烈關(guān)系葛躏。

dispersion(disp):選擇分散值最高的基因

G48E2L1 <- FindVariableFeatures(G48E2L1, selection.method = "vst", nfeatures = 2000)

找出10個差異最大的基因：

top10 <- head(VariableFeatures(G48E2L1),10)
# 繪制帶有和不帶有標簽的變量特性
plot1 <- VariableFeaturePlot(G48E2L1)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))

數(shù)據(jù)的縮放（Scaling the data）

接下來澈段，我們應(yīng)用一個線性變換(“scaling”)，這是一個標準的預(yù)處理步驟舰攒，比PCA等降維技術(shù)更重要败富。

ScaleData函數(shù)功能:

• 改變每個基因的表達，使細胞間的平均表達為0
• 測量每個基因的表達摩窃，使細胞間的差異為1
  * 這一步給予下游分析同等的權(quán)重兽叮，這樣高表達基因就不會占主導(dǎo)地位
• 其結(jié)果存儲在G48E2L1[["RNA"]]@scale.data中

all.genes <- rownames(G48E2L1)
G48E2L1 <- ScaleData(G48E2L1, features = all.genes)

線性降維

接下來，我們對縮放的數(shù)據(jù)執(zhí)行PCA猾愿。默認情況下鹦聪，只使用前面確定的變量特性作為輸入，但是如果您希望選擇不同的子集蒂秘，可以使用features參數(shù)來定義泽本。

G48E2L1 <- RunPCA(G48E2L1, features = VariableFeatures(object = G48E2L1))

Seurat提供了幾種有用的方法來可視化細胞和定義PCA的特性，包括VizDimReduction姻僧、DimPlot和DimHeatmap

檢查和可視化PCA結(jié)果的幾種不同的方法

> print(G48E2L1[["pca"]], dims=1:5, nfeatures = 5)
PC_ 1
Positive:  MALAT1, NEAT1, MT-ND6, XIST, AC004687.1
Negative:  RAN, H2AFZ, PRDX1, SLC25A5, PFN1
PC_ 2
Positive:  EEF1A1, BTF3, C1QBP, UNG, RPL5
Negative:  MKI67, CENPF, TOP2A, ASPM, UBE2C
PC_ 3
Positive:  MIR155HG, PTPRK, LIF, XIST, ANK3
Negative:  S100A4, LTB, S100A11, TMSB4X, FTL
PC_ 4
Positive:  CLDND1, BTG1, ARL6IP1, SELL, STAT1
Negative:  MIF, RPS18, RPS2, GZMB, SRM
PC_ 5
Positive:  GZMB, CCL5, CD8A, CD8B, CTSW
Negative:  HIST1H4C, NME4, CORO1B, STAT1, CDK1

VizDimLoadings(G48E2L1, dims = 1:2, reduction = "pca")

DimPlot(G48E2L1, reduction = "pca")

特別是规丽，DimHeatmap可以方便地探索數(shù)據(jù)集中主要的異構(gòu)來源蒲牧，并且在決定哪些PCs可以用于進一步的下游分析時非常有用。細胞和特征都是根據(jù)它們的PCA分數(shù)排序的嘁捷。將cells設(shè)置為一個數(shù)字造成，可以繪制光譜兩端的“極端”細胞，這極大地加快了繪制大型數(shù)據(jù)集的速度雄嚣。雖然這顯然是一個監(jiān)督分析晒屎，但我們發(fā)現(xiàn)這是一個有價值的工具，用于探索相關(guān)的特征集缓升。

DimHeatmap(G48E2L1, dims = 1, cells = 500, balanced = TRUE)

DimHeatmap(G48E2L1, dims = 1:15, cells = 500, balanced = TRUE)

image.png

確定數(shù)據(jù)的“維數(shù)”

為了克服scRNA-seq數(shù)據(jù)中單個特征中大量的技術(shù)噪聲鼓鲁，Seurat根據(jù)他們的PCA評分將細胞分組，每個PC實質(zhì)上代表一個“元特征”港谊，它將跨相關(guān)特征集的信息組合在一起骇吭。因此，最主要的組件代表了數(shù)據(jù)集的健壯壓縮歧寺。但是燥狰，我們應(yīng)該選擇包含多少個主成分? 10 個？ 20個? 100個斜筐？

在Macosko et al文章中龙致，我們實現(xiàn)了一個重采樣測試的靈感來自JackStraw程序。我們隨機排列數(shù)據(jù)的一個子集(默認為1%)并重新運行PCA顷链，構(gòu)造一個特征得分的“null distribution”目代，然后重復(fù)這個過程。我們認為最“significant” 的PC是那些具有豐富的低p值特征的嗤练。

# Note：對于大型數(shù)據(jù)集榛了，此過程可能需要很長時間，為了方便起見煞抬，請注釋掉霜大。
# 可以使用ElbowPlot()中實現(xiàn)的更近似的技術(shù)來減少計算時間
# 計算時間
G48E2L1 <- JackStraw(G48E2L1, num.replicate = 100)
G48E2L1 <- ScoreJackStraw(G48E2L1, dims = 1:20)

JackStrawPlot函數(shù)提供了一個可視化工具叠穆，用于用均勻分布(虛線)比較每個PC的p-values分布绒窑。“顯著的”PCs將顯示出一個低p值(虛線以上的實線)的強富集特性锦援。在這種情況下蝗碎，在最初的10-12個PCs之后，重要性似乎急劇下降旗扑。

JackStrawPlot(G48E2L1, dims = 1:15)

另一種啟發(fā)式方法生成“Elbow plot”:根據(jù)各成分解釋的方差百分比對主要成分進行排序(ElbowPlot 函數(shù))蹦骑。在這個例子中，我們可以觀察到PC9-10周圍的一個拐點（“elbow”）臀防，這表明大部分真實信號是在前10個pc中捕獲的眠菇。

ElbowPlot(G48E2L1)

對用戶來說边败，確定數(shù)據(jù)集的真實維數(shù)是一項挑戰(zhàn)/不確定的工作。因此捎废，我們建議考慮這三種方法笑窜。第一個是更有監(jiān)督的，探索PCs以確定相關(guān)的異質(zhì)性來源登疗，并可與GSEA聯(lián)合使用排截。第二個實現(xiàn)了一個基于隨機空模型的統(tǒng)計測試，但是對于大型數(shù)據(jù)集來說非常耗時辐益，并且可能不會返回一個明確的PC截止時間断傲。第三種是一種常用的啟發(fā)式算法，可以立即計算出來智政。在這個例子中认罩，所有這三種方法都產(chǎn)生了相似的結(jié)果，但是我們可能有理由選擇PC 7-12之間的任何一個作為截止時間续捂。

我們在這里選擇了10個垦垂，但鼓勵用戶考慮以下幾點:

• 樹突狀細胞和NK細胞愛好者可能認識到，與PCs 12和PCs 13密切相關(guān)的基因定義了罕見的免疫亞群(即MZB1是漿細胞樣dc的標記)牙瓢。但是劫拗，這些組非常罕見，在沒有預(yù)先知識的情況下一罩，很難從這種大小的數(shù)據(jù)集的背景噪聲中區(qū)分出來杨幼。

• 我們鼓勵用戶使用不同數(shù)量的pc(10臺、15臺甚至50臺!)重復(fù)下游分析聂渊。正如你將觀察到的差购，結(jié)果往往沒有顯著的不同。

• 我們建議用戶在選擇這個參數(shù)時汉嗽，不要選得太高欲逃。例如，僅使用5個PCs執(zhí)行下游分析會顯著影響結(jié)果饼暑。

細胞聚類（Cluster the cells）

Seurat v3應(yīng)用了一種基于圖的集群方法稳析，建立在(Macosko等人)的初始策略之上。重要的是弓叛，驅(qū)動聚類分析的距離度量(基于先前確定的PCs)保持不變彰居。然而，我們將細胞距離矩陣劃分成集群的方法已經(jīng)得到了極大的改進撰筷。我們的方法受到最近手稿的很大啟發(fā)陈惰，這些手稿將基于圖的聚類方法應(yīng)用于scRNA-seq數(shù)據(jù) [SNN-Cliq, Xu and Su, Bioinformatics, 2015]和CyTOF數(shù)據(jù) [PhenoGraph, Levine et al., Cell, 2015]。簡單地說毕籽，這些方法將單元格嵌入到一個圖結(jié)構(gòu)中——例如k -最近鄰(KNN)圖抬闯，在具有相似特征表達模式的單元格之間繪制邊緣井辆，然后嘗試將這個圖劃分為高度互連的準團或社區(qū)。

和表現(xiàn)型一樣溶握，我們首先在PCA空間中構(gòu)造一個基于歐氏距離的KNN圖杯缺，然后根據(jù)任意兩個細胞在局部區(qū)域的共享重疊(Jaccard相似性)來細化它們之間的邊權(quán)值。此步驟使用FindNeighbors函數(shù)執(zhí)行睡榆，并將之前定義的數(shù)據(jù)集維度(前10個pc)作為輸入萍肆。

為了對單元進行聚類，我們接下來應(yīng)用模塊化優(yōu)化技術(shù)肉微，如Louvain算法(default)或SLM [SLM, Blondel et al., Journal of Statistical Mechanics]匾鸥，以迭代方式將單元分組在一起，目標是優(yōu)化標準模塊化函數(shù)碉纳。FindClusters函數(shù)實現(xiàn)這個過程勿负，并包含一個分辨率參數(shù)，該參數(shù)設(shè)置下游集群的粒度劳曹，增加的值將導(dǎo)致更多的集群奴愉。我們發(fā)現(xiàn)，將該參數(shù)設(shè)置在0.4-1.2之間铁孵，對于3K左右的單細胞數(shù)據(jù)集通常會得到良好的結(jié)果锭硼。對于較大的數(shù)據(jù)集，最佳分辨率通常會增加蜕劝√赐罚可以使用Idents函數(shù)找到集群。

>G48E2L1 <- FindNeighbors(G48E2L1, dims = 1:10)
>G48E2L1 <- FindClusters(G48E2L1, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7038
Number of edges: 226547

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8435
Number of communities: 9
Elapsed time: 0 seconds

查看前5個細胞的cluster id

> head(Idents(G48E2L1),5)
AAACCTGAGGAATGGA AAACCTGAGGCTCAGA AAACCTGCAGCCTATA AAACCTGGTCAGAATA AAACCTGGTTAAAGAC
               6                1                0                0                4
Levels: 0 1 2 3 4 5 6 7 8

進行非線性降維（UMAP/tSNE）**

Seurat提供了幾種非線性的降維技術(shù)岖沛，如tSNE和UMAP暑始，以可視化和探索這些數(shù)據(jù)集。這些算法的目標是學(xué)習(xí)數(shù)據(jù)的底層流形婴削，以便在低維空間中將相似的單元放在一起廊镜。上面所確定的基于圖的集群中的單元應(yīng)該在這些降維圖上共同定位。作為UMAP和tSNE的輸入唉俗，我們建議使用相同的PCs作為聚類分析的輸入嗤朴。

G48E2L1 <- RunUMAP(G48E2L1, dims = 1:10)
# 注意，您可以設(shè)置'label = TRUE'虫溜，或者使用LabelClusters函數(shù)來幫助標記各個clusters.

DimPlot(G48E2L1, reduction = "umap")

此時可以保存對象雹姊，這樣就可以輕松地將其加載回來，而不必重新運行上面執(zhí)行的計算密集型步驟衡楞，或者輕松地與協(xié)作者共享容为。

saveRDS(G48E2L1, file = "../output/G48E2L1_tutorial.rds")

找差異表達features（cluster biomarkers）

Seurat可以幫助您找到通過差異表達式定義集群的標記。默認情況下，它識別單個簇的陽性和陰性標記(在ident.1中指定)坎背，與所有其他細胞相比較。Findallmarkers為所有集群自動化這個過程寄雀，但是您也可以測試集群組之間的相互關(guān)系得滤，或者測試所有細胞。

min.pct參數(shù)要求至少在兩組細胞中的任何一組中檢測一個特性盒犹，以及thresh.test參數(shù)要求一個特性在兩組之間有一定的差異(平均)懂更。您可以將這兩個值都設(shè)置為0，但是時間上有很大的增加——因為這將測試大量不太可能具有高度歧視性的特性急膀。作為加速這些計算的另一個選項沮协，max.cells.per.ident可以設(shè)置。這將對每個標識類進行采樣卓嫂，使其不具有比設(shè)置的細胞更多的細胞慷暂。雖然通常會有功率的損失，速度的增長可能是顯著的晨雳，最高度差異表達的特征可能仍然會上升到頂部行瑞。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(G48E2L1, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n=5)
               p_val  avg_logFC pct.1 pct.2     p_val_adj
MT-CO1  0.000000e+00  0.8509450 0.995 0.993  0.000000e+00
ACTB    0.000000e+00 -0.8800403 1.000 0.999  0.000000e+00
EEF1A1  0.000000e+00 -0.9152317 0.974 0.992  0.000000e+00
B2M     0.000000e+00 -1.2646953 0.900 0.991  0.000000e+00
RPL28  1.044112e-307  0.6582244 0.993 0.962 3.501744e-303

找出區(qū)分cluster 5與cluster 0和cluster 3的所有標記

cluster5.markers <- FindMarkers(G48E2L1, ident.1 = 5, ident.2 = c(0,3),min.pct = 0.25)
head(cluster5.markers, n=5)
               p_val  avg_logFC pct.1 pct.2     p_val_adj
B2M    6.780927e-144 -0.6074824 0.998 0.998 2.274187e-139
ACTB   3.018783e-139 -0.4931734 1.000 1.000 1.012439e-134
PRDX1  1.243987e-134 -0.6805939 0.797 0.994 4.172084e-130
EEF1A1 9.412366e-126 -0.4098785 1.000 1.000 3.156719e-121
MT-CO1 2.831898e-125  0.5641633 0.992 0.990 9.497620e-121

找出每個cluster的標記與所有剩余的細胞相比較，只報告陽性細胞

G48E2L1.markers <- FindAllMarkers(G48E2L1, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
G48E2L1.markers %>% group_by(cluster) %>% top_n(n=2, wt = avg_logFC)

Seurat有幾個關(guān)于微分表達式的測試餐禁，可以通過該測試設(shè)置血久。使用參數(shù)(詳情請參閱我們的DE vignette)。例如帮非，ROC測試返回任何單個標記(從0 - random到1 - perfect)的分類能力氧吐。

cluster1.markers <- FindMarkers(G48E2L1, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

我們包括一些可視化標記表達的工具。VlnPlot(顯示跨集群的表達式概率分布)和FeaturePlot(在tSNE或PCA圖上可視化特性表達式)是我們最常用的可視化方法末盔。我們還建議使用RidgePlot筑舅、CellScatter和DotPlot作為查看數(shù)據(jù)集的額外方法。

VlnPlot(G48E2L1, features = c("MS4A1", "CD79A"))

# you can plot raw counts as well
VlnPlot(G48E2L1, features = c("NKG7", "PF4"),slot = "counts", log = TRUE)

FeaturePlot(G48E2L1, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

DoHeatmap為給定的細胞和特征生成一個表達式heatmap庄岖。在本例中豁翎，我們繪制每個集群的前20個標記(如果小于20，則繪制所有標記)隅忿。

top10 <- G48E2L1.markers %>% group_by(cluster) %>% top_n(n=10, wt = avg_logFC)
DoHeatmap(G48E2L1, features = top10$gene) + NoLegend()

image.png

為clusters分配細胞類型標識

幸運的是心剥，在這個數(shù)據(jù)集的情況下，我們可以使用規(guī)范的標記背桐，以方便地匹配無偏聚類到已知的細胞類型:

image.png

new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC","Platelet")
names(new.cluster.ids) <- levels(G48E2L1)
G48E2L1 <- RenameIdents(G48E2L1, new.cluster.ids)
DimPlot(G48E2L1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

image.png

saveRDS(G48E2L1, file = "../output/G48E2L1_final.rds")