單細(xì)胞測序--Seurat（上）--創(chuàng)建對象和數(shù)據(jù)預(yù)處理

參考：https://zhuanlan.zhihu.com/p/134124866

需要這三個(gè)文件：

1）barcodes.tsv：樣本中所有的細(xì)胞的條碼庸诱。這文件會(huì)按照 matrix.tsv 的順序來排列

[圖片上傳失敗...(image-b39475-1621137785277)]

2）genes.tsv：參照 (reference name)會(huì)按照Ensembl、NCBI、UCSC網(wǎng)站而有所不同涧衙， gene symbol會(huì)與 matrix.tsv 的順序一致蜘渣。

image.png

3）matrix.mtx：包含 count matrix衷旅，行名與gene.tsv上的行名對應(yīng)，列名與barcodes.tsv相對應(yīng)缕减。

image.png

Seurat 內(nèi)的【Read10X()】函數(shù)可以將上述三文件 raw data整合為一個(gè)稀疏矩陣 (sparse matrix) 庭再，

1捞奕、加載數(shù)據(jù) Load the PBMC dataset

pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")

dim(pbmc.data)

[1] 32738 2700

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

參數(shù) ：

counts: 未標(biāo)準(zhǔn)化的數(shù)據(jù)，如原始計(jì)數(shù)或TPMs

project: 設(shè)置Seurat對象的項(xiàng)目名稱拄轻；默認(rèn)為"SeuratProject"

assay： 與初始輸入數(shù)據(jù)對應(yīng)的分析名稱颅围。

meta.data: 添加到Seurat對象的其他細(xì)胞水平（cell-level）數(shù)據(jù)。一個(gè)matrix恨搓，其中行是cell name谷浅，列是附加的元數(shù)據(jù)字段。

min.cells --該feature至少在n個(gè)細(xì)胞內(nèi)被覆蓋奶卓； 該基因（feature）至少在3個(gè)細(xì)胞中被檢測到

min.features--規(guī)定了至少檢測到這些feature的細(xì)胞一疯。即檢測到的基因至少有200個(gè)細(xì)胞才被用于分析

看看 pbmc 里面有啥

pbmc

An object of class Seurat

13714 features across 2700 samples within 1 assay

Active assay: RNA (13714 features, 0 variable features)

dim(pbmc) #行，列

[1] 13714 2700

head(pbmc@assays))

$RNA

Assay data with 13714 features for 2700 cells

First 10 features:

AL627309.1, AP006222.2, RP11-206L10.2, RP11-206L10.9,

LINC00115, NOC2L, KLHL17, PLEKHN1, RP11-54O7.17, HES4

head(pbmc@meta.data)) QC指標(biāo)存儲

                                orig.ident  nCount_RNA nFeature_RNA

AAACATACAACCAC-1 pbmc3k 2419 779

AAACATTGAGCTAC-1 pbmc3k 4903 1352

AAACATTGATCAGC-1 pbmc3k 3147 1129

AAACCGTGCTTCCG-1 pbmc3k 2639 960

AAACCGTGTATGCG-1 pbmc3k 980 521

AAACGCACTGGTAC-1 pbmc3k 2163 781

head(pbmc@active.assay))

[1] "RNA"

head(pbmc@active.ident))

AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1

pbmc3k pbmc3k pbmc3k pbmc3k

AAACCGTGTATGCG-1 AAACGCACTGGTAC-1

pbmc3k pbmc3k

Levels: pbmc3k

pbmc@project.name

[1] "pbmc3k"

head(rownames(pbmc))

[1] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9"

[5] "LINC00115" "NOC2L"

head(colnames(pbmc))

[1] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1"

[4] "AAACCGTGCTTCCG-1" "AAACCGTGTATGCG-1" "AAACGCACTGGTAC-1"

2夺姑、標(biāo)準(zhǔn)的預(yù)處理流程 Standard pre-processing workflow

1）基于質(zhì)量控制指標(biāo)（QC metrics）和刪除不需要的細(xì)胞作進(jìn)一步分析：

Seurat允許根據(jù)任何用戶定義的標(biāo)準(zhǔn)過濾單元格（cells）墩邀。

通常使用的QC指標(biāo)：

1. 每個(gè)細(xì)胞內(nèi)被檢測到特有的基因（****unique genes****）的數(shù)目，unique feature會(huì)因?yàn)閿?shù)據(jù)質(zhì)量而調(diào)整盏浙。

   ▲低質(zhì)量的細(xì)胞或空的液滴一般含有較少的基因眉睹；

   ▲細(xì)胞雙重態(tài)或多重態(tài)可能呈現(xiàn)異常高的基因count值

2. 細(xì)胞內(nèi)被監(jiān)測到的分子的總數(shù)目（與unique genes高度相關(guān)）

3. 匹配到線粒體基因組的read的百分比

   ▲低質(zhì)量/將要死去的細(xì)胞經(jīng)常呈現(xiàn)過度的線粒體污染情況荔茬；

   ▲使用function計(jì)算線粒體QC metrics, 此function可以計(jì)算源于feature集的count值的百分比

   ▲使用所有以【MT-】為起始的基因集合，作為線粒體基因集

數(shù)據(jù)質(zhì)控

質(zhì)控的參數(shù)主要有兩個(gè)：

1.每個(gè)細(xì)胞測到的unique feature數(shù)目（unique feature代表一個(gè)細(xì)胞檢測到的基因的數(shù)目竹海，可以根據(jù)數(shù)據(jù)的質(zhì)量進(jìn)行調(diào)整）

2.每個(gè)細(xì)胞檢測到的線粒體基因的比例慕蔚，理論上線粒體基因組與核基因組相比，只占很小一部分斋配。所以線粒體基因表達(dá)比例過高的細(xì)胞會(huì)被過濾孔飒。

（參考：https://zhuanlan.zhihu.com/p/145991506）

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

"[[" 操作符可以向?qū)ο笤獢?shù)據(jù)添加列

QC metrics 在 Seurat ：head(pbmc@meta.data)) QC指標(biāo)存儲

▲過濾掉擁有 feature count>2500, or < 200 的細(xì)胞；

▲過濾掉含有> 5% 的線粒體 count值的細(xì)胞

使用 VlnPlot（）進(jìn)行可視化****(基因表達(dá)艰争、指標(biāo)坏瞄、PC分?jǐn)?shù)等)

nFeature RNA: 每個(gè)細(xì)胞所檢測到的基因數(shù)目，也就是以前版本的nGene甩卓；

nCount RNA: 每個(gè)細(xì)胞測到所有基因的表達(dá)量之和鸠匀，即這些基因數(shù)目一共測到的count數(shù)目，也就是以前版本的UMI數(shù)目逾柿；

percent.mt: 每個(gè)細(xì)胞所檢測到的線粒體基因缀棍，即測到的線粒體基因的比例。

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

image.png

FeatureScatter 一般用于可視化 feature-feature 關(guān)系机错，使用點(diǎn)圖查看兩個(gè)數(shù)據(jù)之間的相關(guān)性爬范，也可以用于計(jì)算對象的任何東西，i.e. 對象數(shù)據(jù)中的列毡熏，PC分?jǐn)?shù)等坦敌。

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")

plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1 + plot2

image.png

選擇 200< gene 數(shù)目 <2500（根據(jù)violin圖1） & 線粒體數(shù)目 <5%的細(xì)胞（根據(jù)violin圖2）

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

pbmc

An object of class Seurat

13714 features across 2638 samples within 1 assay #genes:13714, cells：2638

Active assay: RNA (13714 features, 0 variable features)

2）數(shù)據(jù)標(biāo)準(zhǔn)化Normalizing the data

在去除了數(shù)據(jù)集中不需要的細(xì)胞后侣诵，下一步就是標(biāo)準(zhǔn)化數(shù)據(jù)痢法。默認(rèn)設(shè)置下，使用全局-規(guī)模（global-scaling）的標(biāo)準(zhǔn)化方法 "LogNormalize"杜顺，這方法按總體表達(dá)量财搁，對每一個(gè)細(xì)胞的特征表達(dá)量（feature expression）進(jìn)行標(biāo)準(zhǔn)化，乘以一個(gè)比例因子（****scale factor****）（默認(rèn)為10,000）躬络，然后log轉(zhuǎn)換此結(jié)果尖奔。

標(biāo)準(zhǔn)化后的值儲存在【pbmc[["RNA"]]@data】。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization

0% 10 20 30 40 50 60 70 80 90 100%

[----|----|----|----|----|----|----|----|----|----|

高可變度基因的鑒定（基因的選擇） Identification of highly variable features (feature selection)

下一步穷当，計(jì)算數(shù)據(jù)集里面基因子集提茁，該子集呈現(xiàn)出高度的細(xì)胞間變異（在有些細(xì)胞內(nèi)高表達(dá)，在其他細(xì)胞內(nèi)低表達(dá)）馁菜。還發(fā)現(xiàn)了專注在下游分析內(nèi)的基因可幫助突出單細(xì)胞數(shù)據(jù)內(nèi)的生物信號茴扁。

Seurat3 ，根據(jù)舊有版本升級了一下汪疮，通過使用【FindVariableFeatures 】函數(shù)峭火，直接對單細(xì)胞數(shù)據(jù)中固有內(nèi)的均值方差（mean-variance）的關(guān)系建模毁习。在默認(rèn)設(shè)置下，每個(gè)數(shù)據(jù)集返回2,000個(gè)feature卖丸。這些數(shù)據(jù)會(huì)被用于下游數(shù)據(jù)分析纺且，如PCA。

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

#鑒定前10個(gè)高度可變的基因 Identify the 10 most highly variable genes

top10 <- head(VariableFeatures(pbmc), 10)

top10

[1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1"

[9] "GNG11" "S100A8"

#在有l(wèi)abel或無label下稍浆，畫出2000個(gè)高變異的基因

plot1 <- VariableFeaturePlot(pbmc)

plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

plot1 + plot2

或

CombinePlots(plots = list(plot1, plot2))

image

對數(shù)據(jù)去中心化 Scaling the data

使用線性轉(zhuǎn)換(縮放'scaling')载碌，它是一個(gè)降維技巧（如PCA）的標(biāo)準(zhǔn)預(yù)處理步驟，

目的:所謂去中心化粹湃，就是將樣本X中的每個(gè)觀測值都減掉樣本均值恐仑，這樣做的好處是能夠使得求解協(xié)方差矩陣變得更容易。（參考：http://www.reibang.com/p/c2850302b644）

【ScaleData】函數(shù)：

▲轉(zhuǎn)換每個(gè)基因的表達(dá)量为鳄，使得細(xì)胞之間的平均表達(dá)量為0裳仆；

▲按比例縮小每個(gè)基因的表達(dá)量，使得細(xì)胞之間的方差（variance）為1孤钦；

   ?此步驟給予了下游分析的同等權(quán)重（equal weight）歧斟，使得高表達(dá)基因不會(huì)占主要地位；

▲改結(jié)果會(huì)被儲存在【pbmc[["RNA"]]@scale.data】

all.genes <- rownames(pbmc)

pbmc <- ScaleData****(pbmc, features = all.genes)

? 在【Seurat v2】里面偏形，可以使用【ScaleData】函數(shù)來去除單細(xì)胞數(shù)據(jù)集里面一些不需要的變量來源静袖。例如，我們可以“退回（regress out）”與例如細(xì)胞周期階段相關(guān)的異質(zhì)性俊扭，或與線粒體污染相關(guān)的異質(zhì)性队橙。

pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
<meta name="source" content="lake">

單細(xì)胞測序--Seurat（下）--線性降維、識別差異基因萨惑、分配細(xì)胞類型標(biāo)識

1捐康、先跑 PCA 進(jìn)行線性降維 Perform linear dimensional reduction

Seurat 在識別細(xì)胞亞型時(shí)，先用 PCA 挑出幾個(gè)貢獻(xiàn)率最大的主成分庸蔼，再用已選出的主成分分值來將進(jìn)一步聚類分析解总，而非全部細(xì)胞一次性進(jìn)行聚類分析。

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

多種可視化 PCA 結(jié)果的方法 Examine and visualize PCA results a few different ways

（1）

print(pbmc[["pca"]], dims = 1:5, nfeatures = 5) #將其分為5個(gè)主成分

PC_ 1

Positive: CST3, TYROBP, LST1, AIF1, FTL

Negative: MALAT1, LTB, IL32, IL7R, CD2

PC_ 2

Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1

Negative: NKG7, PRF1, CST7, GZMB, GZMA

PC_ 3

Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1

Negative: PPBP, PF4, SDPR, SPARC, GNG11

PC_ 4

Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1

Negative: VIM, IL7R, S100A6, IL32, S100A8

PC_ 5

Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY

Negative: LTB, IL7R, CKB, VIM, MS4A7

（2）

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca") #將其分成兩個(gè)主成分

image

（3）

DimPlot(pbmc, reduction = "pca")

image

（4）

可使用【DimHeatmap】對數(shù)據(jù)集內(nèi)的異質(zhì)性的主要數(shù)據(jù)作簡單探索姐仅，在決定哪個(gè)PC可包含在內(nèi)被用作下有分析時(shí)花枫，此函數(shù)顯得特別有用。細(xì)胞和feature根據(jù)PCA的打分值來設(shè)置掏膏。設(shè)置【cells】為一數(shù)值繪制出范圍兩端的“極端”細(xì)胞劳翰，這可顯著地加速大數(shù)據(jù)集圖表的繪制。盡管這是一監(jiān)督算法分析馒疹，發(fā)現(xiàn)這仍是用于探索相關(guān)feature特征的有價(jià)值的工具佳簸。

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

image

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

image

2、評估最顯著的主成分行冰，決定數(shù)據(jù)集的維數(shù) Determine the 'dimensionality' of the dataset

使用了受JackStraw 程序啟發(fā)的重新取樣測試（ resampling test ）溺蕉。鑒定了“顯著的”PC作為那些擁有低P-value特征的強(qiáng)有力的富集伶丐。

bmc <- JackStraw(pbmc, num.replicate = 100) #此步耗時(shí)較長

pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

【JackStrawPlot】函數(shù)為每一個(gè)PC的P-value的分布與一均勻分布（ uniform distribution）的對比提供了可視化的工具（虛線）。“顯著的”PC會(huì)顯示帶有較低的P-value的特征的強(qiáng)富集（實(shí)線上的虛線）疯特。在此情況下哗魂，在第一個(gè)10~12 PC之后，有一個(gè)顯著的斷崖式下降漓雅。

JackStrawPlot(pbmc, dims = 1:15)

可見下圖前5個(gè)PC 與后5個(gè)PC顯著分開录别，挑選 p-value < 0.05的主成分。

image

3邻吞、細(xì)胞的聚類分析 Cluster the cells

【FindNeighbors】函數(shù)將細(xì)胞距離矩陣分區(qū)為小聚群组题，將一圖劃分為高度互通的“quasi-cliques”或“communities”內(nèi)。該函數(shù)將先前定義好的數(shù)據(jù)集（起初10個(gè)PC）的維度作為輸入抱冷。

pbmc <- FindNeighbors(pbmc, dims = 1:10)

pbmc <- F****indClusters(pbmc, resolution = 0.5)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Ec

Number of nodes: 2638

Number of edges: 96033

Running Louvain algorithm...

Maximum modularity in 10 random starts: 0.8720

Number of communities: 9

Elapsed time: 0 seconds

查看前五個(gè)聚類ID結(jié)果:

head(Idents(pbmc), 5)

AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1

1 3 1 2

AAACCGTGTATGCG-1

6

Levels: 0 1 2 3 4 5 6 7 8

4崔列、UMAP/tSNE的非線性降維分析 Run non-linear dimensional reduction (UMAP/tSNE)

非線性多降維的技巧，如 tSNE 和 UMAP旺遮，用于可視化和探索這些數(shù)據(jù)集赵讯。它們將相似的細(xì)胞放置在低維度的空間，使用相同的 PC 作為聚類分群分析的它們輸入數(shù)據(jù)耿眉”咭恚可設(shè)置 【label = TRUE】或【LabelClusters】函數(shù)幫助標(biāo)記細(xì)胞群。

pbmc <- RunUMAP(pbmc, dims = 1:10)

DimPlot(pbmc, reduction = "umap")

image

5鸣剪、識別差異表達(dá)的 marker 基因 Finding differentially expressed features (cluster biomarkers)

Seurat可尋找通過差異表達(dá)決定分群的 markers组底。【FindAllMarkers】函數(shù)自動(dòng)化完成所有小集群的分組筐骇，你可以測試小集群的分組债鸡，或其他細(xì)胞與細(xì)胞之間的，或針對所有細(xì)胞的分組拥褂∧锼【ident.1】函數(shù)鑒定了單一群落里面的陽性和陰性細(xì)胞（positive and negative markers）牙寞。

【min.pct】需要檢測到兩組之間最小的差異饺鹃，這差異必須是存在一定量的差異表達(dá)。

找出 cluster1 中所有的marker

cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)

head(cluster1.markers, n = 5)

p_val avg_logFC pct.1 pct.2 p_val_adj

IL32 1.894810e-92 0.8373872 0.948 0.464 2.598542e-88

LTB 7.953303e-89 0.8921170 0.981 0.642 1.090716e-84

CD3D 1.655937e-70 0.6436286 0.919 0.431 2.270951e-66

IL7R 3.688893e-68 0.8147082 0.747 0.325 5.058947e-64

LDHB 2.292819e-67 0.6253110 0.950 0.613 3.144372e-63

找出區(qū)別于cluster 0 和cluster 3 與c luster 5的所有標(biāo)記

cluster5.markers <- FindMarkers(pbmc, 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

FCGR3A 7.583625e-209 2.963144 0.975 0.037 1.040018e-204

IFITM3 2.500844e-199 2.698187 0.975 0.046 3.429657e-195

CFD 1.763722e-195 2.362381 0.938 0.037 2.418768e-191

CD68 4.612171e-192 2.087366 0.926 0.036 6.325132e-188

RP11-290F20.3 1.846215e-188 1.886288 0.840 0.016 2.531900e-184

找出與所有剩余細(xì)胞對比的间雀，每個(gè)cluster的標(biāo)記悔详，僅匯報(bào)陽性細(xì)胞

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC) #管道符號

Registered S3 method overwritten by 'cli':

method from

print.boxx spatstat

A tibble: 18 x 7

Groups: cluster [9]

p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene

1 1.96e-107 0.730 0.901 0.594 2.69e-103 0 LDHB

2 1.61e- 82 0.922 0.436 0.11 2.20e- 78 0 CCR7

3 7.95e- 89 0.892 0.981 0.642 1.09e- 84 1 LTB

4 1.85e- 60 0.859 0.422 0.11 2.54e- 56 1 AQP3

5 0. 3.86 0.996 0.215 0. 2 S100A9

6 0. 3.80 0.975 0.121 0. 2 S100A8

7 0. 2.99 0.936 0.041 0. 3 CD79A

8 9.48e-271 2.49 0.622 0.022 1.30e-266 3 TCL1A

9 2.96e-189 2.12 0.985 0.24 4.06e-185 4 CCL5

10 2.57e-158 2.05 0.587 0.059 3.52e-154 4 GZMK

11 3.51e-184 2.30 0.975 0.134 4.82e-180 5 FCGR3A

12 2.03e-125 2.14 1 0.315 2.78e-121 5 LST1

13 7.95e-269 3.35 0.961 0.068 1.09e-264 6 GZMB

14 3.13e-191 3.69 0.961 0.131 4.30e-187 6 GNLY

15 1.48e-220 2.68 0.812 0.011 2.03e-216 7 FCER1A

16 1.67e- 21 1.99 1 0.513 2.28e- 17 7 HLA-DPB1

17 7.73e-200 5.02 1 0.01 1.06e-195 8 PF4

18 3.68e-110 5.94 1 0.024 5.05e-106 8 PPBP

Seurat有多種微分表達(dá)式測試（tests for differential expression），可以在【test.use】參數(shù)內(nèi)被設(shè)置惹挟。例如茄螃，ROC 測試會(huì)給任何獨(dú)立的maker（范圍由0-隨機(jī)，到1-完美）返回一個(gè)“分類能力”连锯。

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

【VlnPlot】顯示了集群里面表達(dá)量的概率分布归苍，【FeaturePlot】（可用于可視化 tSNE 或 PCA圖的特征表達(dá)量）是最經(jīng)常使用的可視化工具用狱。建議探索【 RidgePlot】、【CellScatter】和【DotPlot】 作為另一種方法來探索數(shù)據(jù)集拼弃。

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

image

使用 raw counts 畫圖

VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts****", log = TRUE)

image

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",

•                            "CD8A"))
  

image

【DoHeatmap】對給定的細(xì)胞和特征生成一熱圖夏伊。在本例子中，我們對每個(gè)集群畫出前20個(gè)marker（如果marker量少于20就全部的畫上）

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()

image

6吻氧、給cluster 分配細(xì)胞類型標(biāo)識 Assigning cell type identity to clusters

使用規(guī)范 marker 溺忧，以便匹配無偏倚的聚類到已知的細(xì)胞類型。

image

names(new.cluster.ids) <- levels(pbmc)

pbmc <- RenameIdents(pbmc, new.cluster.ids)

DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

image

saveRDS(pbmc, file = "../pbmc3k_final.rds")
參考:單細(xì)胞測序--Seurat（下）--線性降維盯孙、識別差異基因鲁森、分配細(xì)胞類型標(biāo)識 · 語雀 (yuque.com)