|
由于Immujent最近比較忙，單細(xì)胞流程的教程就擱置了一段時(shí)間，今天小編繼續(xù)拾起之前的剩下的內(nèi)容來(lái)介紹一下Seurat包。很多單細(xì)胞的分析都是建立在確定了細(xì)胞分群（類型）以后，而常用來(lái)做細(xì)胞分群的軟件就是Seurat，在2022年尹始，我們就一起來(lái)學(xué)習(xí)一下Seurat的分析流程吧~
1. 軟件包的安裝和加載
install.packages('Seurat') # 如果已安裝則不需要運(yùn)行
library(Seurat)
library(dplyr)
2. 導(dǎo)入數(shù)據(jù)
這里用到的測(cè)試數(shù)據(jù)都是從Seurat官網(wǎng)下載的，大家有需要的也可以去官網(wǎng)下載（鏈接：https:///seurat/）或者通過(guò)本公眾號(hào)獲取，當(dāng)然如果有自己的數(shù)據(jù)也可以用自己的數(shù)據(jù)來(lái)測(cè)試，具體得到矩陣的方法可參考單細(xì)胞分析流程之Cell Ranger和單細(xì)胞分析流程之Cell Ranger結(jié)果解讀
setwd("J:/scRNA-seq/")
#設(shè)置工作目錄，這是事先創(chuàng)建好的目錄，這樣做的好處是方便進(jìn)行數(shù)據(jù)整理
pbmc.data <- Read10X(data.dir = "hg19")
#導(dǎo)入存放數(shù)據(jù)的目錄，hg19是目錄名
#如果是從GEO等網(wǎng)站中下載的單細(xì)胞表達(dá)矩陣，可以直接讀入表達(dá)矩陣：pbmc.data <- read.table("FPKM.txt",header = TRUE)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
# 創(chuàng)建seurat對(duì)象（10X數(shù)據(jù)和表達(dá)矩陣都可以通過(guò)這條命令創(chuàng)建對(duì)象），
# min.cells：一個(gè)基因至少在多少個(gè)細(xì)胞中表達(dá)（默認(rèn)3）
# min.features：一個(gè)細(xì)胞中最少表達(dá)多少個(gè)基因（默認(rèn)200）
# 一般在最開始創(chuàng)建對(duì)象的時(shí)候不會(huì)過(guò)濾掉太多細(xì)胞，建議使用默認(rèn)參數(shù)即可
pbmc#查看pbmc中存的細(xì)胞數(shù)和基因數(shù)；另外可以用str(pbmc)查看pbmc中的包含的數(shù)據(jù)
如上圖可以看到，一共有13,714個(gè)基因，2,700個(gè)細(xì)胞
3. 質(zhì)控并且挑選用于后續(xù)分析的細(xì)胞
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# 計(jì)算線粒體基因的比例，人中線粒體的基因都是以MT開頭，如果是小鼠的細(xì)胞則將MT替換成Mt即可
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0.1)
# pt.size設(shè)置點(diǎn)的大?。ㄕ故镜氖腔驍?shù)，UMI和線粒體比例的結(jié)果）` </pre>
通過(guò)上面的命令可以先初步看一下數(shù)據(jù)質(zhì)量，然后再挑選過(guò)濾的指標(biāo)。
常用的質(zhì)控信息：
a. 每個(gè)細(xì)胞檢測(cè)到的基因數(shù)目的分布:基因數(shù)過(guò)低的細(xì)胞可能是捕獲率低或細(xì)胞破碎或空beads,基因數(shù)過(guò)高可能是doublets或multiplets。
b. 每個(gè)細(xì)胞中總的UMI(reads)數(shù)的分布：與基因數(shù)一般是很強(qiáng)的正相關(guān)，擁有很少的reads或UMI分子數(shù)的樣品可能是細(xì)胞破損或捕獲失敗，應(yīng)該移除。
c. 線粒體基因表達(dá)所占的比例:一般低質(zhì)量或死亡細(xì)胞線粒體污染比例高一些。
比較簡(jiǎn)單的質(zhì)控方法如下圖（去除異常極端值）：
那么我們接下來(lái)就使用這種過(guò)濾異常極端值的方法來(lái)篩選細(xì)胞：
Q1 <- quantile(pbmc$nFeature_RNA)[2]
Q3 <- quantile(pbmc$nFeature_RNA)[4]
IQR <- Q3 - Q1
upper <- Q3+1.5*IQR
lower <- Q1-1.5*IQR
pbmc <- subset(pbmc, subset = nFeature_RNA > lower & nFeature_RNA < upper & percent.mt < 5)
pbmc
# 查看過(guò)濾完以后查看還剩多少細(xì)胞和基因` </pre>
可以看到，通過(guò)這種過(guò)濾方式一共還剩下13,714個(gè)基因，2,491個(gè)細(xì)胞。
過(guò)濾細(xì)胞這一步是需要反復(fù)調(diào)整的，主要看需要解決的科學(xué)問(wèn)題是什么。如果一個(gè)課題如果想研究細(xì)胞是怎么凋亡的，那么在過(guò)濾的時(shí)候就不能把線粒體的比例作為篩選標(biāo)準(zhǔn)，因?yàn)橐话愕蛲龅募?xì)胞都有較高的線粒體比例，一旦去掉以后則會(huì)把想重點(diǎn)研究的細(xì)胞都去掉，所以過(guò)濾這一步大家一定要反復(fù)測(cè)試。
4. 數(shù)據(jù)標(biāo)準(zhǔn)化
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)`
可以把這一步理解為bulk RNA-seq中計(jì)算的TPM/100，因?yàn)槊總€(gè)細(xì)胞中標(biāo)準(zhǔn)化后的數(shù)據(jù)相加總和為10,000（scale.factor的值），標(biāo)準(zhǔn)化后的數(shù)據(jù)存放在pbmc[["RNA"]]@data，可以通過(guò)輸入pbmc[["RNA"]]@data調(diào)用標(biāo)準(zhǔn)化后的值
如果有心的人會(huì)發(fā)現(xiàn)，如果直接計(jì)算pbmc[["RNA"]]@data（后面就簡(jiǎn)稱為data）每一列的和，值并不是10,000，如下圖：
這里就是其中的一個(gè)重點(diǎn)啦，因?yàn)樵谧鰳?biāo)準(zhǔn)化的時(shí)候選擇的標(biāo)準(zhǔn)化方法是normalization.method = "LogNormalize"，這樣會(huì)在計(jì)算完'TPM'后對(duì)數(shù)據(jù)做log處理，所以直接去計(jì)算單個(gè)細(xì)胞中data的和是不等于10,000的。需要在計(jì)算之前對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換，如下命令：
TPM <- expm1(pbmc[["RNA"]]@data)
# 這樣就能得到單細(xì)胞中的TPM矩陣（每一列的和為10,000）
head(colSums(TPM))
5. 尋找高變基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 默認(rèn)情況會(huì)選擇前2000個(gè)基因，nfeatures可以設(shè)置選擇的基因數(shù)
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) #標(biāo)記出top10的var genes是哪些
plot1
plot2
該圖的X軸是基因在所有細(xì)胞中的均值，Y軸是變異程度，Y軸越大變成程度越大，默認(rèn)會(huì)根據(jù)Y軸選擇前2,000個(gè)基因作為高變基因（紅色的點(diǎn)）。
值得注意的是，這里的高變基因也可以根據(jù)課題的需要自己搜集一些關(guān)鍵基因用來(lái)替換這里的VariableFeatures，在后面涉及到VariableFeatures的地方替換成搜集的基因集即可（這也是根據(jù)課題的需要）
6. 數(shù)據(jù)歸一化
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
# 對(duì)所有基因進(jìn)行歸一化，在做下一步PCA之前需要對(duì)數(shù)據(jù)進(jìn)行歸一化處理
# 數(shù)據(jù)縮放后的結(jié)果保存在pbmc[["RNA"]]@scale.data` </pre>
7. 線性降維（PCA）
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc),verbose = FALSE)
# 通過(guò)設(shè)置 features來(lái)選擇用于降維的基因，也可以輸入自己搜集的基因集（見(jiàn) 5. 尋找高變基因）
DimPlot(object = pbmc, reduction = "pca")
# 繪制細(xì)胞PCA的結(jié)果` </pre>
ElbowPlot(pbmc)
# 選擇用于做后續(xù)分析中的數(shù)據(jù)（PC），一般選取趨緊平緩的點(diǎn)，即可選擇10-15 PC.
# 注意：如果選擇的PC數(shù)太少，會(huì)影響數(shù)據(jù)。如果在不知道怎么選擇時(shí)，可以多測(cè)試一下
8. 細(xì)胞聚類
pbmc <- FindNeighbors(pbmc, dims = 1:10)
# dims：通過(guò)上一步確定的PC，這里選擇的是1~10 PC，大家在做的過(guò)程中可以多調(diào)整一下，可以嘗試1:12,1:15等多個(gè)組合，看看PC的增加對(duì)結(jié)果的影響
pbmc <- FindClusters(pbmc, resolution = 0.7)
# resolution設(shè)置分辨率，即分群的數(shù)目，需不斷調(diào)試；同樣的數(shù)據(jù)resolution越大，分的群越多
Number of communities：分群的個(gè)數(shù)，從上圖可以看到一共分了10群，
9. 可視化（TSNE/UMAP）
pbmc <- RunTSNE(pbmc, dims = 1:10)
# dims中的參數(shù)需要和FindNeighbors中選擇的PC保持一致
DimPlot(pbmc, reduction = "tsne")
# TSNE的結(jié)果
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
# UMAP的結(jié)果
從上面兩個(gè)圖大家可以看出，UMAP的結(jié)果會(huì)相對(duì)來(lái)說(shuō)更集中一些，會(huì)把同一個(gè)群中的細(xì)胞集中在一起。當(dāng)然這個(gè)也要看具體的數(shù)據(jù)，現(xiàn)在主要是推薦用UMAP進(jìn)行可視化，不過(guò)小編覺(jué)得TSNE和UMAP都是可視化的方法，是為數(shù)據(jù)服務(wù)的，所以當(dāng)然是哪個(gè)圖好看就用哪個(gè)圖啦 ~不用一味的追求一定要用某種方法。小編在做課題的時(shí)候就會(huì)兩種都跑一下，然后看哪個(gè)圖好看就會(huì)在文章中放哪個(gè)圖（希望這點(diǎn)能幫到你
10. 尋找每個(gè)cluster的Marker gene并鑒定細(xì)胞類型
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# min.pct：基因在最少多大比例的細(xì)胞中檢測(cè)到才用于后續(xù)的檢測(cè)分析；
# logfc.threshold：高于給定的最小的變化倍數(shù)(log2)的基因用于后續(xù)統(tǒng)計(jì)分析。數(shù)值越大，計(jì)算會(huì)越快，差異基因可能會(huì)少；注意的是seurat v4中計(jì)算的差異倍數(shù)是log2，而seurat v3中計(jì)算的差異倍數(shù)是loge,也就是ln
# 計(jì)算每個(gè)cluster相對(duì)于其他cluster 的marker gene。
# 使用FindAllMarkers就能一步把所有cluster的marker都計(jì)算完，比較簡(jiǎn)單的方式
計(jì)算完每個(gè)群的差異基因就可以用這些差異基因去鑒定細(xì)胞類型啦，鑒定細(xì)胞類型的方法這次就不介紹了。
11. 基因表達(dá)可視化
下面就介紹一些展示基因表達(dá)的方法，可以放到文章中哦~
- 小提琴圖
library(ggplot2)``VlnPlot(object = pbmc, features = c("MS4A1")) + theme(plot.title = element_text(face = "bold.italic"))
# theme(plot.title = element_text(face = "bold.italic")) 使得基因名斜體，不過(guò)這個(gè)命令只能對(duì)單個(gè)基因繪圖時(shí)有用，文章中的基因名要斜體，作圖時(shí)候斜體了就不怕后面忘記修改啦~
- 染色圖
FeaturePlot(object = pbmc,features = c("MS4A1"),cols = c("#C0C0C0", "#E41C12"),reduction = "umap") + theme(plot.title = element_text(face = "bold.italic"))
# reduction：選擇在pca,tsne,umap中對(duì)基因進(jìn)行染色
# cols：顏色的設(shè)置
# pt.size：設(shè)置點(diǎn)的大小，當(dāng)細(xì)胞比較稀疏時(shí)，可以將點(diǎn)設(shè)置較大一些
- 氣泡圖
markers <- c("CD3D", "CD8A", "CD8B", "SELL", "CD4", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1", "GPR183", "PPBP", "GNG11", "TSPAN13", "IL3RA", "IGJ")
DotPlot(object = pbmc, features = markers,cols = c("#C0C0C0", "#E41C12")) + RotatedAxis
# 可用于展示每個(gè)cluster中的marker表達(dá)情況，每一列的點(diǎn)是可以直接比較的（對(duì)列進(jìn)行scale）
- 表達(dá)分布圖
RidgePlot(object = pbmc, features = c("LDHB", "FCGR3A"), ncol = 2)` </pre>
|
