參考：
1.第六章 scRNA-seq數(shù)據(jù)分析
2.初生牛犢--模仿cell research (IF 17+) 文章的一張圖
3.Seurat包學(xué)習(xí)筆記（一）：Guided Clustering Tutorial
4.2019-10-22 R語言Seurat包下游分析-1
目錄
1.數(shù)據(jù)下載
2.數(shù)據(jù)讀取
3.質(zhì)量控制
?3.1 載入注釋
?3.2 質(zhì)控及其細(xì)胞選取
4.標(biāo)準(zhǔn)化
5.高差異基因 (HVGs)
6. 二次標(biāo)準(zhǔn)化
?6.1 縮放的標(biāo)準(zhǔn)
?6.2 刪除不需要的數(shù)據(jù)
7. PCA 分析
8. 決定需要考慮的PC
9.細(xì)胞分集
10.降維后近迁，用聚類的類別可視化【UMAP/tSNE】
11.尋找分集標(biāo)記基因
12.對分集細(xì)胞進(jìn)行標(biāo)記
Seurat是一個常用的clustering的軟件包拴疤。它的最大貢獻(xiàn)在于將tSNE圖引入到scRNA-seq cluter的圖形化展示中來碍讨。tSNE是一種降維算法，它將多維的PCA分析降維到2維圖形中碱鳞，使用人們很方便的在二維圖形中區(qū)分不同類型的細(xì)胞。
安裝并加載所需的R包
- 目前最新版的Seurat包都是基于3.0以上的踱蛀，可以直接通過install.packages命令進(jìn)行安裝窿给。
# 設(shè)置R包安裝鏡像
options(repos="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
install.packages('Seurat')
library(Seurat)
library(dplyr)
library(patchwork)
# 查看Seurat包的版本信息
packageVersion("Seurat")
- 如果想安裝早期2.0版本的Seurat包可以通過devtools包進(jìn)行安裝
# 先安裝devtools包
install.packages('devtools')
# 使用devtools安裝指定版本的Seurat包
devtools::install_version(package = 'Seurat', version = package_version('2.3.4'))
library(Seurat)
Seurat的原教程
1.數(shù)據(jù)下載
10X genomics 數(shù)據(jù)下載地址 : 可以在這里下載到
注釋文件： ensembl數(shù)據(jù)庫中下載 human gtf文件 。為下文計算每個細(xì)胞中不同類型的表達(dá)基因比例做裝備率拒。(如：編碼基因崩泡，線粒體基因，長鏈非編碼基因等等比例)
2.數(shù)據(jù)讀取
Seurat教程中就是10X Genomics數(shù)據(jù)猬膨。通過Read10X來讀取角撞。
library(Seurat)
library(dplyr)
library(ggplot2)
# 讀取Peripheral Blood Mononuclear Cells (PBMC)數(shù)據(jù)
# 解壓文件應(yīng)該是hg19文件夾,切換到這個目錄中：
setwd("C:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
# 直接用Read10X 讀取文件。
pbmc.data <- Read10X(data.dir = "./")
pbmc是一個稀疏矩陣.用. 代替0,優(yōu)化內(nèi)存.
建立Seurat 對象
pbmc <- CreateSeuratObject(counts = pbmc.data, min.cells = 3,
min.features = 200, project = "10X_PBMC")
- counts = pbmc.data : 未標(biāo)準(zhǔn)化的數(shù)據(jù)寥掐，如原始計數(shù)或TPMs .
- min.cells =3 : 基因過濾指標(biāo), 基因至少在3個細(xì)胞中表達(dá).
- min.feature =200 ： 細(xì)胞過濾指標(biāo)靴寂，一個細(xì)胞中表達(dá)基因至少有200.
- project = "10X_PBMC" : 設(shè)置項(xiàng)目名稱，給每一個細(xì)胞打一個標(biāo)簽召耘，都是來自10X_PBMC 項(xiàng)目
獲得的pbmc 對象百炬，有13714個基因，2700個細(xì)胞用過了初始的質(zhì)控標(biāo)準(zhǔn)（min.cells =3;min.feature =200）
Tips: 對于拿到UMI計數(shù)表格的情況污它，我們只要能夠夠建好一個dgTMatrix就可以了剖踊。問題的關(guān)鍵是如何在有限的內(nèi)存下構(gòu)建dgTMatrix∩辣幔可以采取分段讀入的辦法德澈。
比如下面的方法：
library(Seurat)
library(Matrix)
nUMI.summary.file <- "data/umi_expression_matrix.tsv"
counts <- list()
con <- file(nUMI.summary.file, open="r")
## 讀取文件頭
header <- readLines(con, n=1, warn=FALSE)
header <- strsplit(header, "\\t")[[1]]
## 讀取計數(shù)
i <- 1
while(length(buf <- readLines(con, n=1e6, warn=FALSE))>0){
buf <- as.matrix(read.delim(text=buf, header=FALSE, row.names=1))
counts[[i]] <- Matrix(buf, sparse=TRUE)
i <- i + 1
}
counts <- do.call(rbind, counts)
rownames(counts)[1:5]
colnames(counts) <- header
## 初始化一個Seurat對象
d <- CreateSeuratObject(counts = counts, min.cells = 10,
min.features = 200, project = "test_project")
3.質(zhì)量控制
質(zhì)量控制從一開始構(gòu)建Seurat對象的時候就已經(jīng)開始了。那個時候固惯， 我們依據(jù)的是基因表達(dá)的細(xì)胞數(shù)以及細(xì)胞內(nèi)被檢測到的基因數(shù)來進(jìn)行初步篩選的(min.cells =3;min.feature =200)梆造。 Seurat還提供了更多的篩選手段，比如線粒體DNA含量等葬毫。 所以需要gtf注釋文件镇辉，得到哪些基因是線粒體基因，哪些基因是編碼基因.
3.1 載入注釋
對于線粒體基因贴捡，通常為 MT-啟始的忽肛，所以人們常常用這一特點(diǎn)來標(biāo)定線粒體DNA。 但不是全部情形烂斋，下面是示例如何通過 ensembl GFF文件 進(jìn)行讀取文件屹逛。
library(rtracklayer)
# BiocManager::install("rtracklayer")
gtf <- import("Homo_sapiens.GRCh37.87.chr.gtf.gz")
gtf <- gtf[gtf$gene_name!=""]
gtf <- gtf[!is.na(gtf$gene_name)]
## protein coding
protein <-
gtf$gene_name[gtf$transcript_biotype %in%
c("IG_C_gene", "IG_D_gene", "IG_J_gene", "IG_LV_gene",
"IG_M_gene", "IG_V_gene", "IG_Z_gene",
"nonsense_mediated_decay", "nontranslating_CDS",
"non_stop_decay", "polymorphic_pseudogene",
"protein_coding", "TR_C_gene", "TR_D_gene", "TR_gene",
"TR_J_gene", "TR_V_gene")]
## mitochondrial genes
mito <- gtf$gene_name[as.character(seqnames(gtf)) %in% "MT"]
## long noncoding
lincRNA <-
gtf$gene_name[gtf$transcript_biotype %in%
c("3prime_overlapping_ncrna", "ambiguous_orf",
"antisense_RNA", "lincRNA", "ncrna_host", "non_coding",
"processed_transcript", "retained_intron",
"sense_intronic", "sense_overlapping")]
## short noncoding
sncRNA <-
gtf$gene_name[gtf$transcript_biotype %in%
c("miRNA", "miRNA_pseudogene", "misc_RNA",
"misc_RNA_pseudogene", "Mt_rRNA", "Mt_tRNA",
"Mt_tRNA_pseudogene", "ncRNA", "pre_miRNA",
"RNase_MRP_RNA", "RNase_P_RNA", "rRNA", "rRNA_pseudogene",
"scRNA_pseudogene", "snlRNA", "snoRNA",
"snRNA_pseudogene", "SRP_RNA", "tmRNA", "tRNA",
"tRNA_pseudogene", "ribozyme", "scaRNA", "sRNA")]
## pseudogene
pseudogene <-
gtf$gene_name[gtf$transcript_biotype %in%
c("disrupted_domain", "IG_C_pseudogene", "IG_J_pseudogene",
"IG_pseudogene", "IG_V_pseudogene", "processed_pseudogene",
"pseudogene", "transcribed_processed_pseudogene",
"transcribed_unprocessed_pseudogene",
"translated_processed_pseudogene",
"translated_unprocessed_pseudogene", "TR_J_pseudogene",
"TR_V_pseudogene", "unitary_pseudogene",
"unprocessed_pseudogene")]
annotations <- list(protein=unique(protein),
mito=unique(mito),
lincRNA=unique(lincRNA),
sncRNA=unique(sncRNA),
pseudogene=unique(pseudogene))
names(annotations )
annotations <- annotations[lengths(annotations)>0]
3.2 質(zhì)控及其細(xì)胞選取
# Seurat會計算基因數(shù)以及UMI數(shù) (nFeature and nCount).
# Seurat會將原始數(shù)據(jù)保存在RNA slot中础废，
# 每一行對應(yīng)一個基因，每一列對應(yīng)一個細(xì)胞.
slotNames(pbmc)
## [1] "assays" "meta.data" "active.assay" "active.ident" "graphs"
## [6] "neighbors" "reductions" "project.name" "misc" "version"
## [11] "commands" "tools"
names(pbmc@assays) ## Assays(pbmc)
## [1] "RNA"
head(Idents(pbmc)) ## 查看一下當(dāng)前的細(xì)胞分組
## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG
## 10X_PBMC 10X_PBMC 10X_PBMC 10X_PBMC 10X_PBMC
## AAACGCACTGGTAC
## 10X_PBMC
## Levels: 10X_PBMC
下面將運(yùn)用到注釋文件罕模，計算每一個細(xì)胞中各種類型基因评腺，占所有基因比例。比如說編碼基因比例0.85淑掌，線粒體基因比例0.10等等歇僧。
- PercentageFeatureSet ： 函數(shù)就是計算各種基因所占比例
# 在計算比例時，使用目標(biāo)基因中的數(shù)值除以總的數(shù)值锋拖。
# `[[`運(yùn)算符可以給metadata加column.
percent <- lapply(annotations, function(.ele){
PercentageFeatureSet(pbmc, features = rownames(pbmc)[rownames(pbmc) %in% .ele])
})
# AddMetaData 會在object@meta.data中加一列诈悍。這些信息都會在QC中使用到。
for(i in seq_along(percent)){
pbmc[[paste0("percent.", names(percent)[i])]] <- percent[[i]]
}
- 也可以手動對meta.data進(jìn)行編輯. 添加一列對細(xì)胞進(jìn)行標(biāo)示（s_A,s_B）
id <- sample(c("s_A", "s_B"), nrow(pbmc@meta.data), replace = TRUE)
names(id) <- rownames(pbmc@meta.data) #必須有names
pbmc@meta.data$f1 <- id
注釋后結(jié)果如下：
head(pbmc@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.protein percent.mito
## AAACATACAACCAC 10X_PBMC 2419 779 97.02356 3.0177759
## AAACATTGAGCTAC 10X_PBMC 4903 1352 96.45115 3.7935958
## AAACATTGATCAGC 10X_PBMC 3147 1129 93.83540 0.8897363
## AAACCGTGCTTCCG 10X_PBMC 2639 960 98.78742 1.7430845
## AAACCGTGTATGCG 10X_PBMC 980 521 96.53061 1.2244898
## AAACGCACTGGTAC 10X_PBMC 2163 781 91.77069 1.6643551
## percent.lincRNA percent.sncRNA percent.pseudogene f1
## AAACATACAACCAC 84.82844 0.20669698 0 s_B
## AAACATTGAGCTAC 86.33490 0.12237406 0 s_A
## AAACATTGATCAGC 84.30251 0.09532888 0 s_A
## AAACCGTGCTTCCG 75.74839 0.07578628 0 s_B
## AAACCGTGTATGCG 76.73469 0.30612245 0 s_B
## AAACGCACTGGTAC 86.54646 0.04623209 0 s_A
orig.ident 是細(xì)胞標(biāo)示.初始化seurat 時候兽埃，orig.ident 為10x_pbmc 侥钳， 當(dāng)我們需要改變orig.ident時，我們需要改變的是ident slot.
# 試著改變分組柄错。
olID <- Idents(pbmc)
Idents(pbmc) <- as.factor(id)
head(Idents(pbmc))
## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG
## s_B s_A s_A s_B s_B
## AAACGCACTGGTAC
## s_A
## Levels: s_A s_B
Idents(pbmc) <- olID ## 變回來
Tips : 上面是通過自己gtf中提取信息舷夺，進(jìn)行注釋，同時Seurat 里面也可以大體上按照特征進(jìn)行標(biāo)記.
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
可視化開始啦~
- 【一維可視化】對上面pbmc 所有注釋的信息售貌，進(jìn)行可視化
VlnPlot(object = pbmc,
features = c("nFeature_RNA", "nCount_RNA",
paste0("percent.", names(percent))),
ncol = 4)
可以看到各種特征分布情況给猾，進(jìn)行設(shè)置特征的閾值，過濾細(xì)胞颂跨。
- nFeature_RNA ：每個細(xì)胞中檢測到的基因數(shù)目分布情況敢伸。
- nCount_RNA ： 每個細(xì)胞中檢測到的read 數(shù)目分布情況。
- percent.protein : 每個細(xì)胞中編碼蛋白基因所占比例分布情況恒削。
- percent.mito : 每個細(xì)胞中線粒體相關(guān)基因所占比例分布情況池颈。
- ......
-
【二維可視化】選取兩個特征進(jìn)行畫圖 ： FeatureScatter
# FeatureScatter是一個用來畫點(diǎn)圖的工具，可以比較不同的metadata钓丰。 plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mito") plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot1 + plot2image.png
開始真正的細(xì)胞過濾啦~
從VlnPlot中躯砰，我們可以得到大多數(shù)細(xì)胞集中的范圍，通過這些信息携丁，可以設(shè)置過濾條件琢歇。
比如可以使用unique gene counts, percent.protein, percentage.mito來進(jìn)行過濾。
這里我們設(shè)置了200 < nGene < 2500, 0.90 < percent.protein < Inf, -Inf < percent.mito < 0.05.
## 查看維度信息
dim(pbmc)
## [1] 13714 2700
## 開始過濾
pbmc <- subset(x = pbmc,
subset = nFeature_RNA > 200 & nFeature_RNA < 2500 &
percent.protein > 95 &
percent.mito < 5) ## something wrong
## 查看維度信息
dim(pbmc)
## [1] 13714 2142
4.標(biāo)準(zhǔn)化
標(biāo)準(zhǔn)化是一個比較簡單的過程梦鉴，使用的是"logNormalize", 就是將每個基因的表達(dá)量對該細(xì)胞總表達(dá)量進(jìn)行平衡李茫，然后乘以一個因子(scale factor, 默認(rèn)值為10,000), 然后中進(jìn)行對數(shù)轉(zhuǎn)換。
pbmc <- NormalizeData(pbmc)
# 等同于
pbmc <- NormalizeData(seurat.obj,
normalization.method = "LogNormalize", scale.factor = 10000)
# 方法是LogNormalize
# scale.factor是它默認(rèn)的
5.高差異基因 (HVGs)
在單細(xì)胞測序中尚揣，有些基因會在不同的細(xì)胞表現(xiàn)出高度的差異表達(dá)涌矢。Seurate可以對這些基因進(jìn)行標(biāo)記掖举。
? 如果手動計算快骗，相當(dāng)于計算變異系數(shù)，而不是標(biāo)準(zhǔn)差。因?yàn)闃?biāo)準(zhǔn)差對于數(shù)據(jù)呈現(xiàn)對稱分布方篮，及單峰情況體現(xiàn)差異比較好名秀，但是對于雙峰情況，變異系數(shù)cv 可能更好藕溅。
-
Tips:
使用
FindVariableFeatures完成差異分析匕得，選擇數(shù)據(jù)集中差異較高的特征基因（默認(rèn)2000）并用于下游分析。 FindVariableFeatures 在2.0是FindVariablegene巾表。包更新的很快.
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst",
nfeatures = 2000)
# 找到變化最大的十個基因
top10 <- head(VariableFeatures(pbmc), 10)
top10
## [1] "PPBP" "GNLY" "S100A9" "LYZ" "IGLL5" "PF4" "FTL" "FTH1"
## [9] "GNG11" "FCER1A"
- selection.method ： 默認(rèn)是“vst” ,也可以選擇"roc"
- nfeatures : 需要得到多少個高變異的基因汁掠。此處為2000個基因，下圖可視化中標(biāo)紅的基因集币。
對HVGs 進(jìn)行可視化
plot1 <- VariableFeaturePlot(pbmc)
LabelPoints(plot = plot1, points = top10, repel = TRUE)
- 紅色的HVGs 有2000個考阱，從上面 FindVariableFeatures ，命令參數(shù)就可以知道.
6. 二次標(biāo)準(zhǔn)化
從上圖中鞠苟，我們可以看到很多基因都表現(xiàn)出很大的差異性乞榨，但是這些差異性并不一定都是事實(shí)，有些可能是因?yàn)閷?shí)驗(yàn)当娱，或者數(shù)據(jù)采集的方式導(dǎo)入的吃既。因此人們想盡量的消除這些因素帶來的差異。這些差異有可能是因?yàn)閎atch-effect跨细， 或者cell alignment rate鹦倚，或者其它。這些因素都可以寫入metadata中冀惭，然后通過ScaleData函數(shù)來消除申鱼。
? ScaleData 可以通過metadata信息，消除很多其他因素影響云头。
6.1 縮放的標(biāo)準(zhǔn)
? 應(yīng)用線性變換縮放數(shù)據(jù)捐友，這是一個標(biāo)準(zhǔn)的預(yù)處理步驟，比PCA等降維技術(shù)更重要溃槐。使用ScaleData函數(shù):
（1）使每個基因在所有細(xì)胞間的表達(dá)量均值為0
（2）使每個基因在所有細(xì)胞間的表達(dá)量方差為1
縮放操作給予下游分析同等的權(quán)重匣砖，這樣高表達(dá)基因就不會占主導(dǎo)地位，其結(jié)果存儲在pbmc[["RNA"]]@scale.data中昏滴，這里默認(rèn)ScaleData 進(jìn)行了中心化猴鲫，標(biāo)準(zhǔn)化.
## 全局縮放
all.genes <- rownames(pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
head(pbmc[["RNA"]]@scale.data[, 1:3])
## AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG
## AL627309.1 -0.06452497 -0.06452497 -0.06452497
## AP006222.2 -0.03726409 -0.03726409 -0.03726409
## RP11-206L10.2 -0.04153370 -0.04153370 -0.04153370
## RP11-206L10.9 -0.03734170 -0.03734170 -0.03734170
## LINC00115 -0.08343360 -0.08343360 -0.08343360
## NOC2L -0.32126699 -0.32126699 -0.32126699
6.2 刪除不需要的數(shù)據(jù)
使用ScaleData函數(shù)從單細(xì)胞數(shù)據(jù)集中刪除不需要的變量。例如谣殊，我們可以“regress out”與細(xì)胞周期階段或線粒體污染相關(guān)的異質(zhì)性（去除不需要的變量拂共，即校正協(xié)變量，校正線粒體基因比例的影響）姻几。如下：
# 刪除不需要的數(shù)據(jù)宜狐，校正線粒體基因引起的影響
# 這個步驟非常耗時
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
Tips
我們對矩陣進(jìn)行了細(xì)胞水平質(zhì)控势告，但是如何對基因進(jìn)行質(zhì)控，當(dāng)然我們創(chuàng)建seurat 對象時候抚恒，已經(jīng)設(shè)定一個cutoff. 基因必須早多少個細(xì)胞中表達(dá). 但是如果直接用這些基因進(jìn)行聚類和降維（T-sne/UMAP 可視化 可能引起偏差咱台，所以需要先進(jìn)行一下PCA 初步選取主要的主成分，對應(yīng)著權(quán)重不同的基因俭驮，再用這些基因進(jìn)行聚類回溺，降維可視化效果更好）
7. PCA 分析
在PCA分析中，前文標(biāo)記的高差異基因object@var.genes會拿來做分析的數(shù)據(jù)混萝。當(dāng)然人們也可以指定自己希望的使用的基因遗遵。當(dāng)細(xì)胞數(shù)越多，使用標(biāo)記的高差異基因越接近考慮全基因組時得到的結(jié)論逸嘀。所以瓮恭，在這一步的時候我們可以看出，之前標(biāo)記的高差異基因如果數(shù)量過低的話厘熟，會對這一步的影響很大屯蹦。
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(pbmc))
## PC_ 1
## Positive: CST3, TYROBP, AIF1, FTL, LST1, FTH1, LYZ, FCER1G, TYMP, FCN1
## S100A9, CFD, LGALS1, CTSS, SERPINA1, S100A8, LGALS2, COTL1, PSAP, IFITM3
## SPI1, SAT1, CFP, S100A11, NPC2, IFI30, LGALS3, GSTP1, S100A6, OAZ1
## Negative: LTB, IL32, IL7R, CD2, B2M, ACAP1, STK17A, CTSW, AQP3, GIMAP5
## RARRES3, CCL5, HOPX, GZMA, SELL, CST7, PPP2R5C, BEX2, CD8B, MYC
## TNFAIP8, ETS1, GZMK, TCF7, ZAP70, RORA, PRF1, KLRG1, GIMAP7, LDLRAP1
......
......
......
## PC_ 5
## Positive: LTB, IL7R, CKB, VIM, AQP3, MS4A7, CYTIP, RGS10, RP11-290F20.3, CD2
## SIGLEC10, HMOX1, HN1, NDUFA4, ATP1A1, TRADD, LILRB2, PTGES3, FYB, ANXA5
## TCF7, CORO1B, FAM110A, IL32, ABRACL, PPA1, TIMP1, GDI2, GIMAP7, VMO1
## Negative: GZMB, FGFBP2, NKG7, GNLY, CST7, PRF1, CCL4, GZMA, SPON2, GZMH
## S100A8, XCL2, CTSW, CLIC3, AKR1C3, CCL5, CCL3, IGFBP7, XCL1, S100A9
## LGALS2, GPR56, HOPX, TTC38, GSTP1, RBP7, TYROBP, CD14, S100A12, MS4A6A
Seurat提供了多種手段來查看PCA分析的結(jié)果，其中包括PrintPCA, VizPCA, PCAPlot, 以及PCHeatmap
【一維可視化】 對每一個維度绳姨，權(quán)重較大的基因進(jìn)行可視化
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
【二維可視化】 將細(xì)胞降維到二維空間
DimPlot(object = pbmc, reduction = "pca")
Tips: 下圖以PC_1 來看登澜，是解釋細(xì)胞分布最大的基因組合方向，所有的基因是維度飘庄，有些基因?qū)@個方向共享很大袱衷，也就是說更容易分開細(xì)胞丹禀，這些基因可能更加重要斋扰。下面是熱圖1展示這些基因?qū)τ趨^(qū)分細(xì)胞的效果陵刹，可以看到將很多細(xì)胞區(qū)分成紫色和黃色.效果不錯.
可以看到5,6維度貌似區(qū)分不開了.
DimHeatmap(pbmc, dims = 1:12, cells = 500, balanced = TRUE)
8. 決定需要考慮的PC
Suerat使用了jackStraw方法[4]來估計應(yīng)該使用多少個principal components的基因來進(jìn)行下游分析。
# 注意: 這可能是一個非常費(fèi)時的過程碾盐。
pbmc <- JackStraw(pbmc, num.replicate = 100) # 重復(fù)一百次
pbmc <- ScoreJackStraw(pbmc, dims = 1:20) # 選擇20個PC
JackStrawPlot以及PCElbowPlot函數(shù)提供給我們兩種可視化手段來查看jackStraw的結(jié)果晃跺。使用JackStrawPlot可以查看每個PC的p value，這樣方便我們通過p值進(jìn)行篩選毫玖。PCElbowPlot可以方便的查看從哪個PC開始掀虎，差分度就進(jìn)入的平臺期。
- 1.
JackStrawPlot
JackStrawPlot(pbmc, dims = 1:15)
-
2.
PCElbowPlot: 拐點(diǎn)圖 # 最大可設(shè)置到50ElbowPlot(pbmc,ndims = 50, reduction = "pca") # 最大可設(shè)置到50可以直觀看到7 后面共享比較平均.
先擇PC付枫，對于得到正確的結(jié)果具體關(guān)鍵性的作用烹玉。除了上面的方法以外，還可以使用其它的一些辦法阐滩，比如說使用GSEA等富集分析的辦法來選擇PC二打，或者說使用一個隨機(jī)模型來查看具體哪些PC會比較有效果。但是后者這一方法在實(shí)現(xiàn)起來需要消耗過多的資源掂榔。
在本例中继效，因?yàn)槭荢eurat挑選的例子症杏，所以通過上面的JackStraw方法，只要把cut.off值設(shè)置在7-10之間莲趣，就可以得到差不多的結(jié)果。
Tips: 確定進(jìn)行聚類和降維的基因和細(xì)胞后饱溢，下一步就是聚類了喧伞，在這里需要搞情況聚類和降維作用是不一樣的。特別是T-sne是降維可視化作用绩郎，而不是用它來聚類.為什么看到T-sne 圖是不同顏色的色塊呢潘鲫？常見的聚類包括kmeans/hclust .
9.細(xì)胞分集
當(dāng)決定了使用哪些PC中的基因?qū)?xì)胞進(jìn)行分類之后，就可以使用FindClusters來對細(xì)胞分集了肋杖。分集后的結(jié)果會保存在object@ident slot中溉仑。
（Step-1）先在PCA空間中構(gòu)造一個基于歐氏距離的KNN圖，然后根據(jù)任意兩個細(xì)胞在局部鄰近區(qū)域的共享重疊(Jaccard相似性)來細(xì)化它們之間的邊權(quán)值状植。此步驟使用FindNeighbors函數(shù)執(zhí)行浊竟，并將之前定義的數(shù)據(jù)集維度作為輸入。
（Step-2）使用FindClusters函數(shù)對數(shù)據(jù)進(jìn)行優(yōu)化津畸，并包含一個分辨率參數(shù)振定，該參數(shù)設(shè)置下游集群的“granularity”，增加的值將導(dǎo)致更多的集群肉拓。將該參數(shù)設(shè)置在0.4-1.2之間后频，對于3K左右的單細(xì)胞數(shù)據(jù)集通常會得到良好的結(jié)果。對于較大的數(shù)據(jù)集暖途，最佳分辨率通常會增加卑惜。可以使用Idents函數(shù)找到clusters驻售。
? 下面選取了1到10 維度的基因進(jìn)行聚類分析
pbmc <- FindNeighbors(pbmc, dims = 1:10) ## dims 作為輸入的維度
pbmc <- FindClusters(pbmc, resolution = 0.6) ## 分辨率參數(shù)露久，和cluster 個數(shù)有關(guān)系，如何你想獲得粗糙的分類欺栗，需要調(diào)成大一些.
# 查看一下前幾個的cluster,結(jié)果聚成3類.可能發(fā)現(xiàn)新的細(xì)胞系.
head(Idents(pbmc), 5)
## AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG AAACCGTGTATGCG AAACGCTGGTTCTT
## 4 3 0 6 4
## Levels: 0 1 2 3 4 5 6 7 8
10.降維后抱环，用聚類的類別可視化【UMAP/tSNE】
使用UMAP或者tSNE的好處在于將多維的PC圖降維至2維圖像中來，這樣方便分集細(xì)胞的可視化纸巷。 UMAP皺縮成團(tuán)分不開镇草，而tSNE圖則分散得多。
pbmc <- RunUMAP(pbmc, dims = 1:10)
P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
pbmc <- RunTSNE(pbmc, dims = 1:10)
P2 <- DimPlot(pbmc, reduction = "tsne",label = TRUE)
P1 + P2
得到聚類分組結(jié)果瘤旨，通常需要繼續(xù)分析梯啤，在每組里面那些標(biāo)記基因。由于10x genomic 對每一個細(xì)胞都進(jìn)行標(biāo)簽化存哲，不知道每一個標(biāo)簽對應(yīng)那種細(xì)胞因宇，可以通過每組的標(biāo)記基因?qū)γ拷M進(jìn)行細(xì)胞類型判別（根據(jù)文獻(xiàn)說明的已知marker）,進(jìn)而發(fā)現(xiàn)新的標(biāo)記基因.
11.尋找分集標(biāo)記基因
生成了tSNE圖像后七婴，我們還需要確認(rèn)這一堆堆的細(xì)胞都分別是些什么細(xì)胞，這時察滑，可以使用特定細(xì)胞的marker來對不同的細(xì)胞群進(jìn)行標(biāo)記打厘。但是想要對它們進(jìn)行標(biāo)記，首先需要得到每個集群差異表達(dá)的基因贺辰。
在比較時户盯，Seurat使用一比多的辦法。如果需要對cluster 1的細(xì)胞(ident.1)進(jìn)行差異分析時饲化，它比較的對象將是除了cluster 1以外的其它所有的細(xì)胞莽鸭。也可以使用ident.2參數(shù)來傳遞需要比較的對象。 在FindMarkers函數(shù)中吃靠，min.pct是指的marker基因所占細(xì)胞數(shù)的最低比例硫眨，這里使用1/4。 使用max.cells.per.ident可以讓函數(shù)對細(xì)胞進(jìn)行抽樣巢块，以加速計算礁阁。 FindAllMarkers可以一次性找到所有的高表達(dá)的marker。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## RPS12 4.371171e-118 0.5264969 1 0.990 5.994624e-114
## RPS6 1.681008e-114 0.5234111 1 0.994 2.305334e-110
## RPS27 1.719993e-107 0.5242477 1 0.991 2.358799e-103
## RPS14 1.546910e-105 0.4500125 1 0.993 2.121432e-101
## RPS25 5.848778e-98 0.5427189 1 0.973 8.021015e-94
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = 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 3.849393e-139 2.723391 0.981 0.098 5.279058e-135
## RHOC 1.149445e-91 1.665090 0.865 0.130 1.576350e-87
## CDKN1C 6.799344e-75 1.061123 0.519 0.022 9.324621e-71
## IFITM2 1.347125e-73 1.560680 1.000 0.644 1.847448e-69
## HES4 3.616785e-69 1.107357 0.596 0.052 4.960059e-65
計算所有類和其他類的標(biāo)記基因
# find markers for every cluster compared to all remaining cells,
# report only the positive ones
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE,
min.pct = 0.25, thresh.use = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
## # A tibble: 18 x 7
## # Groups: cluster [9]
## p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 0. 3.72 0.969 0.136 0. 0 S100A8
## 2 6.30e-298 3.77 0.996 0.236 8.63e-294 0 S100A9
## 3 7.44e- 94 0.807 0.938 0.587 1.02e- 89 1 LDHB
## 4 3.14e- 61 0.848 0.426 0.105 4.31e- 57 1 CCR7
## 5 1.22e- 80 0.975 0.981 0.635 1.67e- 76 2 LTB
## 6 5.02e- 69 0.956 0.783 0.315 6.88e- 65 2 IL7R
## 7 0. 2.95 0.931 0.045 0. 3 CD79A
## 8 5.62e-214 2.46 0.617 0.024 7.71e-210 3 TCL1A
## 9 3.86e-164 2.39 0.973 0.22 5.30e-160 4 CCL5
## 10 2.97e-142 2.12 0.595 0.052 4.07e-138 4 GZMK
## 11 1.66e-171 2.34 0.981 0.135 2.27e-167 5 FCGR3A
## 12 1.68e-106 1.90 1 0.368 2.30e-102 5 LST1
## 13 2.24e-211 3.48 0.972 0.063 3.07e-207 6 GZMB
## 14 1.84e-127 3.52 0.943 0.138 2.52e-123 6 GNLY
## 15 1.11e-180 2.68 0.806 0.013 1.52e-176 7 FCER1A
## 16 5.01e- 20 1.93 1 0.548 6.87e- 16 7 HLA-DPB1
## 17 4.08e-171 5.02 1 0.012 5.60e-167 8 PF4
## 18 6.09e- 99 5.96 1 0.026 8.35e- 95 8 PPBP
還可以使用不同的算法來尋找差異表達(dá)的基因族奢。
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 0,
thresh.use = 0.25, test.use = "roc",
only.pos = TRUE)
想了解更多關(guān)于查找差異表達(dá)基因的信息氮兵，可以查看Seurat的vignette。
【一維可視化】標(biāo)記基因在不同cluster 表達(dá)量
拿到差異表達(dá)基因后歹鱼，可以使用VlnPlot來查看結(jié)果泣栈。
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw UMI counts as well
VlnPlot(object = pbmc, features= c("NKG7", "PF4"),
slot = "counts", log = TRUE)
## RidgePlot
RidgePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
## DotPlot
DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
【二維可視化】標(biāo)記基因在不同細(xì)胞表達(dá)量
FeaturePlot(object = pbmc,
features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A",
"FCGR3A", "LYZ", "PPBP", "CD8A"),
cols = c("blue", "red"))
提取所有cluster 的top10 標(biāo)記基因可視化
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
# 設(shè)置slim.col.label=TRUE將避免打印每個細(xì)胞的名稱，而只打印cluster的ID.
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
目前我們通過聚類知道了弥姻，將細(xì)胞分成8個類別南片，但是并不知道具體是那種細(xì)胞。所有需要進(jìn)一步通過已知的特定細(xì)胞標(biāo)記基因?qū)luster 打上細(xì)胞類型標(biāo)簽.
12.對分集細(xì)胞進(jìn)行標(biāo)記
對于示例數(shù)據(jù)庭敦，我們已經(jīng)有一些比較常見的細(xì)胞marker了疼进，很方便的就可以對分集的細(xì)胞進(jìn)行標(biāo)記了。
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | IL7R, S100A4 | Memory CD4+ |
| 2 | CD14, LYZ | CD14+ Mono |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
將原來的cluster 1-8 秧廉，修改成對應(yīng)的細(xì)胞名稱
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(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
P2 <- DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 0.5) + NoLegend()
P1 + P2
當(dāng)我們修改了聚類的參數(shù)伞广，結(jié)果會不同。比如分辨率參數(shù).
# 生成一個快照, 其實(shí)就是把pbmc@ident中的數(shù)據(jù)拷貝到pbmc@meta.data中疼电。
# 當(dāng)我們想把pbmc@meta.data中的值賦值給pbmc@ident的話嚼锄，可以使用SetAllIdent函數(shù)。
# 比如pbmc <- SetAllIdent(object=pbmc, id = 'orig.ident')
pbmc <- StashIdent(object = pbmc, save.name = "ClusterNames_0.6")
# 現(xiàn)在我們使用不同的resolution來重新分群蔽豺。當(dāng)我們提高resolution的時候区丑，
# 我們可以得到更多的群。
# 之前我們在FindClusters中設(shè)置了save.snn=TRUE，所以可以不用再計算SNN了沧侥。
# 下面就直接提高一下區(qū)分度可霎，設(shè)置resolution = 0.8
pbmc <- FindClusters(pbmc, resolution = 0.8)
分辨率左圖是0.8，右圖0.6
# 并排畫兩個tSNE,右邊的是用ClusterNames_0.6 (resolution=0.6) 來分組顏色的宴杀。
# 我們可以看到當(dāng)resolution提高之后癣朗，CD4 T細(xì)胞被分成了兩個亞群。
plot1 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE) + NoLegend()
plot2 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE,
group.by = "ClusterNames_0.6") +
NoLegend()
plot1 + plot2
類別不同旺罢，我們又可以按照這一次聚類結(jié)果旷余，尋找marker基因
# 我們再一次尋找不同分集中的marders。
tcell.markers <- FindMarkers(object = pbmc, ident.1 = 0, ident.2 = 1)
# Most of the markers tend to be expressed in C1 (i.e. S100A4).
# However, we can see that CCR7 is upregulated in
# C0, strongly indicating that we can differentiate memory from naive CD4 cells.
# cols.use demarcates the color palette from low to high expression
FeaturePlot(object = pbmc, features = c("S100A4", "CCR7"),
cols = c("blue", "red"))
