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