Seurat standard pipeline(10核心流程)

創(chuàng)建Seurat對象 Read10X CreateSeuratObjecet
質(zhì)控 PercentageFeatureSubject subset
標(biāo)準(zhǔn)化Normalization NormalizeData
高變基因選擇 FindVariableFeatures
數(shù)據(jù)縮放 ScaleData
線性降維 RunPCA
維數(shù)選擇 FindNeighbors
細胞聚類 FindClusters
非線性降維（UMAP/tSNE）RunUMAP DimPlot
鑒定差異表達特征（cluster markers）
細胞注釋 CellMarker---FeaturePlot---RenameIdents

下載示例數(shù)據(jù)：

Peripheral Blood Mononuclear Cells (PBMC) 是10X Genomics dataset page提供的一個數(shù)據(jù)，包含2700個單細胞鸟缕，出自Illumina NextSeq 500平臺粒梦。

PBMCs是來自健康供體具有相對少量RNA（around 1pg RNA/cell）的原代細胞阶牍。在Illumina NextSeq 500平臺晨另，檢測到2700個單細胞弦牡，每個細胞獲得69000 reads灭将。

tar -xvf pbmc3k_filtered_gene_bc_matrices.tar

文件夾下包含3個文件

barcodes.tsv
genes.tsv
matrix.mtx

一本慕、創(chuàng)建 Seurat 分析對象

1.1 Read10X() 函數(shù)可以自動讀入和解析數(shù)據(jù)

library(dplyr)
library(Seurat)
library(patchwork)

#讀取PBMC數(shù)據(jù)
pbmc.data <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")

#查看數(shù)據(jù)
dim(pbmc.data)
# 32738  2700

pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
# 3 x 30 sparse Matrix of class "dgCMatrix"

# CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
# TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
# MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .
#.表示0

summary(colSums(pbmc.data))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 548    1758    2197    2367    2763   15844

查看每個細胞有多少基因被檢測到

at_least_one <- apply(pbmc.data, 2, function(x) sum(x>0))


hist(at_least_one, breaks = 100,
     main = "Distribution of detected genes",
     xlab = "Genes with at least one tag")

hist(colSums(pbmc.data), breaks = 100, 
     main = "Expression sum per cell",
     xlab = "Sum expression")

1.2 pbmc 數(shù)據(jù)初始化 Seurat 對象

使用 CreateSeuratObject 創(chuàng)建對象揍庄。

sce<- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
sce
# An object of class Seurat 
# 13714 features across 2700 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

head(sce$RNA@data[,1:5])
6 x 5 sparse Matrix of class "dgCMatrix"
              AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
AL627309.1                   .                .                .                .                .
AP006222.2                   .                .                .                .                .
RP11-206L10.2                .                .                .                .                .
RP11-206L10.9                .                .                .                .                .
LINC00115                    .                .                .                .                .
NOC2L                        .                .                .                .                .

#過濾檢測 
min.features = 200：一個細胞最少要檢測到200個基因王污，
min.cells = 3：一個基因最少得在4個細胞中表達

其他參數(shù)解釋

• counts：未標(biāo)準(zhǔn)化的數(shù)據(jù)硝桩，如原始計數(shù)或TPMs
• project = "CreateSeuratObject"：設(shè)置Seurat對象的項目名稱
• assay = "RNA"：與初始輸入數(shù)據(jù)對應(yīng)的分析名稱
• names.field = 1：對于每個cell的初始標(biāo)識類路翻，從cell的名稱中選擇此字段椭豫。例如捌朴，如果cell在輸入矩陣中被命名為BARCODE_CLUSTER_CELLTYPE津函，則設(shè)置名稱蛉幸。字段設(shè)置為3以將初始標(biāo)識設(shè)置為 #CELLTYPE破讨。
• names.delim = "_"：對于每個cell的初始標(biāo)識類，從cell的列名中選擇此分隔符奕纫。例如提陶，如果cell命名為bar - cluster - celltype，則將此設(shè)置為“-”匹层，以便將cell名稱分離到其組成部分中隙笆，以選擇相關(guān)字段。
• meta.data = NULL：要添加到Seurat對象的其他單元級元數(shù)據(jù)升筏。應(yīng)該是data.frame撑柔，其中行是單元格名稱，列 #是附加的元數(shù)據(jù)字段您访。

這里還涉及到的知識點是R data中的S3類和S4類铅忿。
list一般情況下被認為是S3類，S4是指使用slots儲存數(shù)據(jù)的格式灵汪。

> sce
An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)

#1個數(shù)據(jù)集檀训，包含2700個細胞柑潦，13714個基因。

.在矩陣中的值表示0(未檢測到分子)峻凫。

由于scRNA-seq矩陣中的大多數(shù)值為0,因此Seurat在任何可能的情況下都使用稀疏矩陣表示妒茬。這為Drop-seq/inDrop/10x數(shù)據(jù)節(jié)省了大量內(nèi)存和速度。

所謂稀疏矩陣蔚晨，也就是在矩陣中乍钻，若數(shù)值為0的元素數(shù)目遠多于非0元素的數(shù)目，并且非0元素分布沒有規(guī)律時铭腕，則稱該矩陣為稀疏矩陣银择；與之相反，若非0元素數(shù)目占大多數(shù)時累舷，則稱該矩陣為稠密矩陣浩考。

二、數(shù)據(jù)預(yù)處理

這部分是基于數(shù)據(jù)質(zhì)控方法被盈，標(biāo)準(zhǔn)化和檢測到的變化基因?qū)?shù)據(jù)進行篩選析孽。

1. 對細胞的質(zhì)控 subset()

可以參考文章：Classification of low quality cells from single-cell RNA-seq data

常用的質(zhì)控指標(biāo)：
（1）每個細胞在檢測到的特異基因數(shù)

• 低質(zhì)量的細胞以及空泡油滴中一般檢測到很少的基因
• 包含多個細胞的油滴會檢測到異常多的基因

（2）每個細胞檢測到的分子總數(shù)（與基因密切相關(guān)）
（3）每個細胞的線粒體基因比例。低質(zhì)量/瀕死細胞常表現(xiàn)出廣泛的線粒體污染

使用PercentageFeatureSet()函數(shù)計算線粒體QC指標(biāo)只怎，使用所有以MT-開頭的基因作為一組線粒體基因

（1）線粒體基因占比計算

PercentageFeatureSet（）參數(shù)意義:
pattern = "^MT-"為正則表達式寫法袜瞬，意為匹配"MT-"開頭的基因，即線粒體基因身堡。區(qū)分MT大小寫邓尤，此處如果小寫則代表另一個物種。

為什么將線粒體基因作為一個過濾條件贴谎？“因為線粒體是參與細胞凋亡啟動和執(zhí)行的主要細胞器之一汞扎。線粒體基因高表達意為著樣品質(zhì)量差，導(dǎo)致大量細胞凋亡或裂解擅这。以及特定樣本的生物學(xué)澈魄，例如腫瘤活檢，可能由于代謝活動和/或壞死而增加線粒體基因表達仲翎。對于線粒體基因比例過高的細胞痹扇，會干擾細胞分群，在分析過程中為了避免這些細胞的影響谭确，通常會根據(jù)實際情況設(shè)計一個閾值進行細胞的過濾

###向 sce 新增一列 percent.mt 數(shù)據(jù)
sce[['percent.mt']] <- PercentageFeatureSet(sce,pattern='^MT-')

#展示前5個細胞的QC指標(biāo)
head(sce@meta.data, 5)

                orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898

QC指標(biāo)儲存在哪帘营？
每個細胞基因數(shù)和總分子數(shù)在建立seurat對象時就已經(jīng)自動計算好了。

可以使用sce@meta.data$percent.mt或者 sce['percent.mt']查看每個細胞的線粒體比例

（2）畫圖查看基因數(shù)目, UMI數(shù)目, 線粒體基因占比

# 使用小提琴圖可視化QC指標(biāo)
VlnPlot(sce, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3)

nFeature_RNA： 代表每個細胞測到的基因數(shù)目逐哈。
nCount_RNA ：代表每個細胞測到所有基因的表達量之和芬迄。
percent.mt：代表測到的線粒體基因的比例

（3）查看基因數(shù)目, 線粒體基因占比與UMI數(shù)目的關(guān)系
FeatureScatter通常用于可視化 feature-feature 相關(guān)性，

#nCount_RNA 與 percent.mt的相關(guān)性
plot1 <- FeatureScatter(sce, 
                        feature1 = "nCount_RNA", 
                        feature2 = "percent.mt")

#nCount_RNA 與 nFeature_RNA的相關(guān)性
plot2 <- FeatureScatter(sce, 
                        feature1 = "nCount_RNA", 
                        feature2 = "nFeature_RNA")

plot1 + plot2 #合并兩圖

拿 nCounts_RNA 與 nFeature_RNA 的散點圖（scatter）來說：

• 每個點代表一個細胞昂秃，斜率代表隨著count的增加gene的增加程度禀梳。
• count和gene一般呈現(xiàn)線性關(guān)系杜窄，斜率越大也就是較少的count就可以檢出較多的gene，說明這批細胞基因的豐度偏低（普遍貧窮）算途；反之塞耕，斜率較小，說明這批細胞基因豐度較高（少數(shù)富有）嘴瓤。
• 有的時候不是一條可擬合的線扫外，或者是兩條可擬合的線，也反映出細胞的異質(zhì)性廓脆。
• 總之筛谚，他就是一個散點圖，描述的是兩個變量的關(guān)系**停忿。

過濾線粒體基因表達比例過高的細胞驾讲，和一些極值細胞（可以根據(jù)小提琴圖判斷，查看兩端離群值）席赂。
（4）質(zhì)控 subset()

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

nCount_RNA:每個細胞的UMI數(shù)量
nFeature_RNA:每個細胞的基因數(shù)
選取 2500 > nFeature_RNA >200 和percent.mt < 5的數(shù)據(jù)

三吮铭、數(shù)據(jù)標(biāo)準(zhǔn)化 NormalizeData()

因為在測序之前會對抓取到的RNA進行PCR擴增，所以需要考慮文庫深度對測序的影響颅停。例：細胞A中總reads數(shù)10谓晌，基因a的reads數(shù)為2；細胞B中總reads數(shù)1000便监，基因a的reads數(shù)為20扎谎。并不能說基因a在細胞B中表達量大于A，故需要進行標(biāo)準(zhǔn)化烧董。所以需要對上一步得到的稀疏矩陣進行Normalize。

Normalize方式：每個細胞每個基因的特征計數(shù)除以該細胞（一列）的特征總計數(shù)胧奔，再乘以scale.factor(默認10,000)逊移，然后對結(jié)果進行l(wèi)og1p對數(shù)轉(zhuǎn)換（log1p=log(n+1)）。標(biāo)準(zhǔn)化的數(shù)值存儲在sce[["RNA"]]@data中龙填。

scale.factor = 10,000 的原因是

進行l(wèi)og轉(zhuǎn)換的原因是：
以bulk-rna-seq為例胳泉，用FPKM或者TPM繪制熱圖的時候，因為數(shù)值變化的范圍太過巨大岩遗，都需要進行一個log轉(zhuǎn)換扇商，讓數(shù)據(jù)壓縮在一個區(qū)間內(nèi)。其次宿礁，也是最重要的改變數(shù)據(jù)分布：測序數(shù)值本身不符合正態(tài)分布案铺，log轉(zhuǎn)換能讓數(shù)據(jù)趨近于正態(tài)分布，以便于后續(xù)進一步分析

默認使用數(shù)據(jù)標(biāo)準(zhǔn)化方法是LogNormalize, 每個細胞總的表達量都標(biāo)準(zhǔn)化到10000梆靖，然后log取對數(shù)控汉；標(biāo)準(zhǔn)化后的數(shù)據(jù)保存在結(jié)果存放于sce[["RNA"]]@data笔诵。，原始數(shù)保存在sce[["RNA"]]@counts中

標(biāo)準(zhǔn)化前姑子，每個細胞總的表達量

hist(colSums(sce$RNA@data),
     breaks = 100,
     main = "Total expression before normalisation",
     xlab = "Sum of expression")

標(biāo)準(zhǔn)化

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

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

若所有調(diào)用的參數(shù)都是默認值乎婿，則可省去

sce <- NormalizeData(sce)

標(biāo)準(zhǔn)化后，每個細胞總的表達量

hist(colSums(sce$RNA@data),
     breaks = 100,
     main = "Total expression after normalisation",
     xlab = "Sum of expression")  

標(biāo)準(zhǔn)化前和標(biāo)準(zhǔn)化后對比

四街佑、 鑒定高變基因 FindVariableFeatures()

接下來谢翎，鑒定在細胞間表達高度變化的基因 (即，它們在某些細胞中高表達沐旨，而在其他細胞中低表達)岳服。后續(xù)研究需要集中于這部分基因這些基因有助于突出單細胞數(shù)據(jù)集中的生物信號。

用FindVariableFeatures()函數(shù)實現(xiàn)希俩，首先計算每一個基因的均值和方差吊宋，并且直接模擬其關(guān)系。默認返回2000個features 颜武。這些將用于下游分析璃搜，如PCA。

高變基因方法選擇vst
1.使用loss(局部加權(quán)回歸)擬合平滑曲線模型
2.獲取模型計算的值作為y = var.exp值
3.var.standarlized = get variance after feature standardization: (每個基因-mena)/sd后取var()鳞上，注意sd=sqrt(var.exp)
4.按照var.standalized降序排列这吻，取前n(比如2000)個基因作為高變基因

sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|


# 查看最高變的10個基因
top10 <- head(VariableFeatures(sce), 10)
#head(sce$RNA@var.features,10)
# "PPBP"   "LYZ"    "S100A9" "IGLL5"  "GNLY"   "FTL"    "PF4"    "FTH1"   "GNG11"  "S100A8"


# 畫出不帶標(biāo)簽或帶標(biāo)簽基因點圖
plot1 <- VariableFeaturePlot(sce)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

五、數(shù)據(jù)縮放（歸一化）

線性變換(“縮放”)篙议，是在PCA降維之前的一個標(biāo)準(zhǔn)預(yù)處理步驟唾糯。最終每個基因均值為0，方差為1鬼贱。

ScaleData()函數(shù)功能:

• 轉(zhuǎn)換每個基因的表達值移怯，使每個細胞的平均表達為0
• 轉(zhuǎn)換每個基因的表達值，使細胞間的方差為1

此步驟在下游分析中具有相同的權(quán)重这难，因此高表達的基因不會占主導(dǎo)地位
結(jié)果存儲在sce[["RNA"]]@scale.data中

all.genes <- rownames(sce)
sce<- ScaleData(sce, features = all.genes)
Centering and scaling data matrix
  |================================================================| 100%

做了 Normalize 還做 scale 的原因是：
normalize改變的是數(shù)據(jù)的分布舟误，scale改變的是數(shù)據(jù)的范圍。
normalize是將偏態(tài)分布轉(zhuǎn)換成趨近于正態(tài)分布姻乓，sclae是將數(shù)據(jù)的分布限定在一個范圍內(nèi)嵌溢。
R語言的Z score計算是通過[scale()]函數(shù)求得，Seurat包中ScaleData()函數(shù)也基本參照了scale()函數(shù)的功能蹋岩。
scale方法中的兩個參數(shù)：center和scale

• center和scale默認為真赖草，即T或者TRUE
• center為真表示數(shù)據(jù)中心化：數(shù)據(jù)集中的各項數(shù)據(jù)減去數(shù)據(jù)集的均值。如有數(shù)據(jù)集1, 2, 3, 6, 3剪个，其均值為3 那么中心化之后的數(shù)據(jù)集為1-3,2-3,3-3,6-3,3-3秧骑，即：-2,-1,0,3,0 *scale為真表示數(shù)據(jù)標(biāo)準(zhǔn)化：指中心化之后的數(shù)據(jù)在除以數(shù)據(jù)集的標(biāo)準(zhǔn)差，即數(shù)據(jù)集中的各項數(shù)據(jù)減去數(shù)據(jù)集的均值再除以數(shù)據(jù)集的標(biāo)準(zhǔn)差。如有數(shù)據(jù)集1, 2, 3, 6, 3腿堤，其均值為3阀坏，其標(biāo)準(zhǔn)差為1.87 那么標(biāo)準(zhǔn)化之后的數(shù)據(jù)集為(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87，即：-1.069,-0.535,0,1.604,0

Z score的概念是指原始數(shù)據(jù)距離均值有多少個標(biāo)準(zhǔn)差笆檀。當(dāng)以標(biāo)準(zhǔn)差為單位進行測量時忌堂，Z score 衡量的是一個數(shù)值偏離總體均值以上或以下多少個標(biāo)準(zhǔn)差。如果原始數(shù)值高于均值酗洒，則 Z score得分為正士修，如果低于均值，則Z score為負樱衷。Z score其實是標(biāo)準(zhǔn)正態(tài)分布（Standard Normal Distribution）棋嘲，即平均值μ=0，標(biāo)準(zhǔn)差 σ=1 的正態(tài)分布矩桂。SND標(biāo)準(zhǔn)正態(tài)分布的直方圖如下所示：

六沸移、線性降維

1、PCA

接下來侄榴，對縮放的數(shù)據(jù)執(zhí)行PCA雹锣。默認情況下，只使用前面確定的變量特性作為輸入癞蚕，但是如果想選擇不同的子集蕊爵，可以使用features參數(shù)來定義(選擇高變基因)。

> sce <- RunPCA(sce,features = VariableFeatures(object = sce))
PC_ 1
Positive:  CST3, TYROBP, LST1, AIF1, FTL, FTH1, LYZ, FCN1, S100A9, TYMP
           FCER1G, CFD, LGALS1, S100A8, CTSS, LGALS2, SERPINA1, IFITM3, SPI1, CFP
           PSAP, IFI30, SAT1, COTL1, S100A11, NPC2, GRN, LGALS3, GSTP1, PYCARD
Negative:  MALAT1, LTB, IL32, IL7R, CD2, B2M, ACAP1, CD27, STK17A, CTSW
           CD247, GIMAP5, AQP3, CCL5, SELL, TRAF3IP3, GZMA, MAL, CST7, ITM2A
           MYC, GIMAP7, HOPX, BEX2, LDLRAP1, GZMK, ETS1, ZAP70, TNFAIP8, RIC3
PC_ 2
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1, HLA-DRA, LINC00926, CD79B, HLA-DRB1, CD74
           HLA-DMA, HLA-DPB1, HLA-DQA2, CD37, HLA-DRB5, HLA-DMB, HLA-DPA1, FCRLA, HVCN1, LTB
           BLNK, P2RX5, IGLL5, IRF8, SWAP70, ARHGAP24, FCGR2B, SMIM14, PPP1R14A, C16orf74
Negative:  NKG7, PRF1, CST7, GZMB, GZMA, FGFBP2, CTSW, GNLY, B2M, SPON2
           CCL4, GZMH, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX
           TTC38, APMAP, CTSC, S100A4, IGFBP7, ANXA1, ID2, IL32, XCL1, RHOC
PC_ 3
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1, HLA-DPA1, CD74, MS4A1, HLA-DRB1, HLA-DRA
           HLA-DRB5, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, CD37, HVCN1, FCRLA, IRF8
           PLAC8, BLNK, MALAT1, SMIM14, PLD4, LAT2, IGLL5, P2RX5, SWAP70, FCGR2B
Negative:  PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU
           HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, PTCRA, CA2, ACRBP, MMD, TREML1
           NGFRAP1, F13A1, SEPT5, RUFY1, TSC22D1, MPP1, CMTM5, RP11-367G6.3, MYL9, GP1BA
PC_ 4
Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1, CD74, HLA-DPB1, HIST1H2AC, PF4, TCL1A
           SDPR, HLA-DPA1, HLA-DRB1, HLA-DQA2, HLA-DRA, PPBP, LINC00926, GNG11, HLA-DRB5, SPARC
           GP9, AP001189.4, CA2, PTCRA, CD9, NRGN, RGS18, GZMB, CLU, TUBB1
Negative:  VIM, IL7R, S100A6, IL32, S100A8, S100A4, GIMAP7, S100A10, S100A9, MAL
           AQP3, CD2, CD14, FYB, LGALS2, GIMAP4, ANXA1, CD27, FCN1, RBP7
           LYZ, S100A11, GIMAP5, MS4A6A, S100A12, FOLR3, TRABD2A, AIF1, IL8, IFI6
PC_ 5
Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY, CCL4, CST7, PRF1, GZMA, SPON2
           GZMH, S100A9, LGALS2, CCL3, CTSW, XCL2, CD14, CLIC3, S100A12, CCL5
           RBP7, MS4A6A, GSTP1, FOLR3, IGFBP7, TYROBP, TTC38, AKR1C3, XCL1, HOPX
Negative:  LTB, IL7R, CKB, VIM, MS4A7, AQP3, CYTIP, RP11-290F20.3, SIGLEC10, HMOX1
           PTGES3, LILRB2, MAL, CD27, HN1, CD2, GDI2, ANXA5, CORO1B, TUBA1B
           FAM110A, ATP1A1, TRADD, PPA1, CCDC109B, ABRACL, CTD-2006K23.1, WARS, VMO1, FYB

查看對每個主成分影響比較大的基因集

> print(sce [["pca"]], dims = 1:5, nfeatures = 5)
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 

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

VizDimReduction, DimPlot, 和 DimHeatmap可以從基因或細胞角度可視化pca結(jié)果

可視化對每個主成分影響比較大的基因集

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

DimPlot(sce , reduction = "pca")

DimHeatmap繪制基于單個主成分的熱圖桦山，細胞和基因的排序都是基于他們的PCA分數(shù)攒射。對于數(shù)據(jù)異質(zhì)性的探索是很有幫助的，可以幫助用戶選擇哪些PC維度可以用于下一步的下游分析恒水。

DimHeatmap(sce, dims = 1, cells = 500, balanced = TRUE) #1個PC 500個細胞

展示多個主成分

DimHeatmap(sce,dims = 1:15,cells = 500,balanced = TRUE) #15個PC

2会放、數(shù)據(jù)維度

主成分分析的原理非常簡單，概括來說就是選擇包含信息量大的維度（features）寇窑，去除信息量少的“干擾”維度鸦概。所以這里會有個問題——如何知道保留幾個維度是最佳的呢？我們希望通過保留盡可能少的維度來留存盡可能多的信息甩骏。Seurat有兩種方法來確定維度。

JackStrawPlot()函數(shù)提供可視化方法先慷，用于比較每一個主成分的p-value的分布饮笛，虛線是均勻分布；顯著的主成分富集有小p-Value基因论熙，實線位于虛線左上方福青。下圖表明保留10個pca主成分用于后續(xù)分析是比較合理的。

JackStraw && Elboew plot

sce <- JackStraw(sce,num.replicate = 100)
sce <- ScoreJackStraw(sce,dims = 1:20)

> JackStrawPlot(sce,dims = 1:15)
Warning message:
Removed 23516 rows containing missing values (`geom_point()`).

主要看在哪個PC之后，顯著性大幅下降无午，也就是前多少個維度包含了大部分的樣本信息媒役。可以看出宪迟，在10-12個PC之后酣衷，顯著性大幅下降，也就是前10-12個維度包含了大部分的樣本信息次泽。

ElbowPlot(sce)

主要看PC在哪個附近有拐點（“elbow”）穿仪，拐點表明大部分真實信號是在前多少個pc中捕獲的。

可以看出意荤，PC9-10附近有一個拐點（“elbow”）啊片，這表明大部分真實信號是在前10個pc中捕獲的。

綜合以上方法玖像，選擇10個主成成分作為參數(shù)用于后續(xù)分析紫谷。

啟發(fā)式評估方法生成一個Elbow plot圖。在圖中展示了每個主成分對數(shù)據(jù)方差的解釋情況（百分比表示）捐寥，并進行排序笤昨。根據(jù)自己需要選擇主成分，圖中發(fā)現(xiàn)第9個主成分是一個拐點上真，后續(xù)的主成分(PC)變化都不大了咬腋。

注意：鑒別數(shù)據(jù)的真實維度不是件容易的事情；除了上面兩種方法睡互，Serat官當(dāng)文檔還建議將主成分(數(shù)據(jù)異質(zhì)性的相關(guān)來源有關(guān))與GSEA分析相結(jié)合根竿。Dendritic cell 和 NK aficionados可能識別的基因與主成分 12 和 13相關(guān)，定義了罕見的免疫亞群 (i.e. MZB1 is a marker for plasmacytoid DCs)就珠。如果不是事先知道的情況下寇壳，很難發(fā)現(xiàn)這些問題。

Serat官當(dāng)文檔因此鼓勵用戶使用不同數(shù)量的PC（10妻怎、15壳炎，甚至50）重復(fù)下游分析。其實也將觀察到的逼侦，結(jié)果通常沒有顯著差異匿辩。因此，在選擇此參數(shù)時榛丢，可以盡量選大一點的維度铲球，維度太小的話對結(jié)果會產(chǎn)生不好的影響。

七晰赞、 細胞聚類

Seurat v3 應(yīng)用基于圖形的聚類方法稼病，例如KNN方法选侨。具有相似基因表達模式的細胞之間繪制邊緣，然后將他們劃分為一個內(nèi)聯(lián)群體然走。

在PhenoGraph中援制，首先基于pca維度中（先前計算的pca數(shù)據(jù)）計算歐式距離（the euclidean distance），然后根據(jù)兩個細胞在局部的重合情況（Jaccard 相似系數(shù)）優(yōu)化兩個細胞之間的邊緣權(quán)值芍瑞。此步驟內(nèi)置于FindNeighbors函數(shù)晨仑，輸入時先前確定的pc數(shù)據(jù)。

為了聚類細胞啄巧，接下來應(yīng)用模塊化優(yōu)化技術(shù)迭代將細胞聚集在一起寻歧。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters 函數(shù)實現(xiàn)這一功能秩仆，其中需要注意resolution參數(shù)码泛，該參數(shù)設(shè)置下游聚類分析的“granularity”，更大的resolution會導(dǎo)致更多的細胞類群澄耍。3000左右的細胞量噪珊，設(shè)置resolution為0.4-1.2是比較合適的。細胞數(shù)據(jù)集越大齐莲，需要更大的resolution參數(shù), 會獲得更多的細胞聚類痢站。
查看細胞屬于那個類群可以使用函數(shù)Idents。

> sce<- FindNeighbors(sce, dims = 1:10) #選取前10個主成分來分類細胞
Computing nearest neighbor graph
Computing SNN


> sce<- FindClusters(sce, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2638
Number of edges: 95965

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

#查看細胞屬于那個類群
head(Idents(sce), 5)
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 
               0                3                2                5                6 
Levels: 0 1 2 3 4 5 6 7 8

八、非線性降維(UMAP/tSNE)

Seurat 提供了一些非線性降維度分析的方法，如 tSNE 和 UMAP鞋喇，以可視化和探索這些數(shù)據(jù)集；這一步使用的數(shù)據(jù)建議與聚類分析使用的pc維度一致呜叫。

# install UMAP： reticulate::py_install(packages ='umap-learn')

sce <- RunUMAP(sce,dims = 1:20)
DimPlot(sce,reduction = 'umap')

# 添加細胞類群 的標(biāo)簽
DimPlot(sce,reduction = 'umap', label = TRUE)
LabelClusters(DimPlot(sce, reduction = "umap"),id = 'ident')

此時可以保存數(shù)據(jù)，方便下次直接導(dǎo)入數(shù)據(jù)修改圖形殿衰。

saveRDS(sce, file = "../output/pbmc_tutorial.rds")

九朱庆、找差異表達基因 ( 聚類標(biāo)志cluster biomarkers )

利用 FindMarkers 命令，可以找到各個細胞類型中與其他類別的差異表達基因闷祥，作為該細胞類型的生物學(xué)標(biāo)記基因娱颊。

• ident.1：設(shè)置待分析的細胞類別
• min.pct：設(shè)定在兩個細胞群中任何一個被檢測到的百分比，通過此設(shè)定不檢測很少表達基因來縮短程序運行時間凯砍。默認0.1箱硕。
• thresh.test：設(shè)定在兩個細胞群中基因差異表達量（平均）∥蝰茫可以設(shè)置為0 颅痊，程序運行時間會更長。
• max.cells.per.ident：每個類群細胞抽樣設(shè)置局待；也可以縮短程序運行時間。通過降低每個類的采樣值，提高計算速度

# find all markers of cluster 1
cluster1.markers <- FindMarkers(sce, 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


# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(sce, 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


# 找出每個cluster的標(biāo)記與所有剩余的細胞相比較钳榨，只報告陽性細胞
# find markers for every cluster compared to all remaining cells, report only the positive ones
sce.markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
sce.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)


       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene    
       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
 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可以通過參數(shù)test.use設(shè)定檢驗差異表達的方法(詳情見 DE vignett)舰罚。

cluster1.markers <- FindMarkers(sce, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
      myAUC  avg_diff power avg_log2FC pct.1 pct.2
RPS12 0.824 0.5094804 0.648  0.7350248 1.000 0.991
RPS6  0.821 0.4712446 0.642  0.6798622 1.000 0.995
RPS27 0.819 0.4996079 0.638  0.7207819 0.999 0.992
RPL32 0.815 0.4238952 0.630  0.6115515 0.999 0.995
RPS14 0.808 0.4296946 0.616  0.6199183 1.000 0.994

FindAllMarkers（）參數(shù)意義：

• only.pos = TRUE：只尋找上調(diào)的基因
• min.pct = 0.1：某基因在細胞中表達的細胞數(shù)占相應(yīng)cluster細胞數(shù)最低10%
• logfc.threshold = 0.25 ：fold change倍數(shù)為0.25
  該函數(shù)是決速步，執(zhí)行比較耗時薛耻。

有多種方法可視化標(biāo)記基因的方法

• VlnPlot: 基于細胞類群的基因表達概率分布
• FeaturePlot：在tSNE 或 PCA圖中畫出基因表達情況
• RidgePlot营罢，CellScatter，DotPlot

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

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

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

DoHeatmap為指定的細胞和基因花表達熱圖饼齿。每個類群默認展示top 10 標(biāo)記基因饲漾。

#每個聚類前10個差異基因表達熱圖(如果小于10，則繪制所有標(biāo)記)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(sce, features = top10$gene) + NoLegend()

十缕溉、鑒定細胞類型

數(shù)據(jù)集的 markers 可以通過查閱相關(guān)文獻及網(wǎng)站 (CellMarker) 人工注釋考传，或者利用singleR自動注釋。singleR 比較雞肋证鸥，最好的辦法目前還是查閱文件和cellmarker(http://biocc.hrbmu.edu.cn/CellMarker/)進行區(qū)分僚楞，比較考驗個人及課題組經(jīng)驗及搜索能力。

比如研究多形性膠質(zhì)母細胞瘤(GBM)枉层，下圖中的文獻給出了利用單細胞轉(zhuǎn)錄組測序確定膠質(zhì)母細胞瘤中腫瘤泉褐、基質(zhì)及免疫細胞亞群：通過scRNA測序確定了GBM中CD45-和CD45+細胞(即非免疫和免疫細胞群)。表達EGFR和iCre的簇定為腫瘤細胞鸟蜡。采取Johnson-Verhaak命名法標(biāo)記EGFR+腫瘤簇膜赃。通過marker確定細胞類型，如Oligodendrocytes(Mog揉忘，Plp1)跳座，Microglia(P2ry12，Tmem119)癌淮，Astrocytes(Aqp4)躺坟，Perivascular macrophages(Cd163，Mrc1)乳蓄，cDC1(Xcr1)等咪橙。

那我們就可以根據(jù)文中的基因試一下：

FeaturePlot(sce,
            features = c('MOG','PIP1','AQP4','CD163','P2RY12','TMEM119','XCR1'),
            reduction = 'umap')

列表

Cluster ID Markers Cell Type
8 P2ry12 Microglia

十一、細胞注釋

Assigning cell type identity to clusters
根據(jù)已知細胞類型與基因標(biāo)記的對應(yīng)關(guān)系虚倒，可以為細胞類群找到對應(yīng)的細胞類型美侦。

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", 
    "NK", "DC", "Platelet")

names(new.cluster.ids) <- levels(sce)
sce<- RenameIdents(sce, new.cluster.ids)
DimPlot(sce, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

# 保存分析的結(jié)果
saveRDS(sce, file = "../output/sce_final.rds")

官網(wǎng)此部分介紹：
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

參考：
http://www.reibang.com/p/728125be7c53
http://www.reibang.com/p/5b26d7bc37b7
https://zhuanlan.zhihu.com/p/359020041
http://www.reibang.com/p/728125be7c53
http://www.reibang.com/nb/52917361
http://www.reibang.com/p/fef17a1babc2
Single cell