scRNA基礎(chǔ)分析-1:安裝包、導(dǎo)入數(shù)據(jù)、過濾質(zhì)控 - 簡書
scRNA基礎(chǔ)分析-2：降維與聚類 - 簡書
scRNA基礎(chǔ)分析-3：鑒定細(xì)胞類型 - 簡書
scRNA基礎(chǔ)分析-4：細(xì)胞亞類再聚類、注釋 - 簡書
scRNA基礎(chǔ)分析-5：偽時間分析 - 簡書
scRNA基礎(chǔ)分析-6：富集分析 - 簡書

library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
scRNA <- readRDS("scRNA.rds")

1、首先選擇2000個（默認(rèn)）表達(dá)值變化大的基因代表細(xì)胞轉(zhuǎn)錄譜

• Seurat包負(fù)責(zé)篩選高變基因的函數(shù)是FindVariableFeatures()，它并不刪除scRNA對象中的非高變基因。

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

• 找出的結(jié)果可以通過VariableFeatures()函數(shù)獲取

top10 <- head(VariableFeatures(scRNA), 10) 
top10
# [1] "GNLY"   "IGLC2"  "S100A9" "IGLC3"  "FCGR3A" "S100A8" "CDKN1C"
# [8] "GZMB"   "ITM2C"  "LYZ" 

• 可視化

plot1 <- VariableFeaturePlot(scRNA) 
LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
#如下圖，橫坐標(biāo)是某基因在所有細(xì)胞中的平均表達(dá)值，縱坐標(biāo)是此基因的方差;紅點即為高變基因（2000個）

1-1

補(bǔ)充：細(xì)胞周期相關(guān)基因

• 上一步找到的高變基因，可能會包含一些細(xì)胞周期相關(guān)基因；
• 它們會導(dǎo)致細(xì)胞聚類發(fā)生一定的偏移，即相同類型的細(xì)胞在聚類時會因為細(xì)胞周期的不同而分開。
• 因此有必要查看是否有細(xì)胞周期相關(guān)基因的存在；若有，則剔除

#細(xì)胞周期有關(guān)基因
head(c(cc.genes$s.genes,cc.genes$g2m.genes))
# [1] "MCM5" "PCNA" "TYMS" "FEN1" "MCM2" "MCM4"

#查看我們選擇的高變基因中有哪些細(xì)胞周期相關(guān)基因
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))

• 在scRNA@meta.data中添加S.Score、G2M.Score和Phase三列有關(guān)細(xì)胞周期的信息。

g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = rownames(scRNA))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA))
scRNA <- CellCycleScoring(object=scRNA,  g2m.features=g2m_genes,  s.features=s_genes)
head(scRNA@meta.data)

1-2

• 觀察細(xì)胞周期相關(guān)基因是否影響聚類

scRNAa <- RunPCA(scRNA, features = c(s_genes, g2m_genes))
DimPlot(scRNA, reduction = "pca", group.by = "Phase")
#如下圖結(jié)果，細(xì)胞周期基因?qū)?xì)胞聚類的影響不大，不需要去除。

# 如需去除，代碼如下
# scRNAb <- ScaleData(scRNA, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA))

1-3

2、PCA降維（線性降維）

• 使用主成分分析將2000維的信息投射到50個維度，并提取前10-20個維度（人工選擇）的信息代表細(xì)胞的轉(zhuǎn)錄特征

scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 
plot1 <- DimPlot(scRNA, reduction = "pca", group.by="orig.ident")
#（左圖）根據(jù)主成分1和2的值將細(xì)胞在平面上展示出來
plot2 <- ElbowPlot(scRNA, ndims=20, reduction="pca") 
#（右圖）展示前20個主成分的解釋量
plot1+plot2

• 重點關(guān)注下右圖，后續(xù)分析要根據(jù)右圖選擇提取的pc軸數(shù)量，一般選擇斜率平滑的點之前的所有pc軸；此圖，作者的建議是選擇前18個pc軸。

  
  
  2-1

3、非線性降維（tSNE或UMAP）

• 最后使用非線性降維方法（tSNE或UMAP）將這10-20個PC值降維到二維空間。
• 經(jīng)過上述操作，轉(zhuǎn)錄譜的特征信息會損失一些，但是大部分轉(zhuǎn)錄特征會在二維空間呈現(xiàn)出來。

先聚類cluster

pc.num=1:18
scRNA <- FindNeighbors(scRNA, dims = pc.num) 
# dims參數(shù)，需要指定哪些pc軸用于分析；這里利用上面的分析，選擇18
scRNA <- FindClusters(scRNA, resolution = 0.5)
# resolution參數(shù)，需要指定0.1-0.9之間的一個數(shù)值，用于決定clusters的相對數(shù)量；
#數(shù)值越大，cluters越多。
table(scRNA@meta.data$seurat_clusters) #分成了0-9，共10個cluster
#  0   1   2   3   4   5   6   7   8   9 
#296 184 120 119  70  55  53  50  44  22 

tSNE

scRNA = RunTSNE(scRNA, dims = pc.num)
embed_tsne <- Embeddings(scRNA, 'tsne')
plot1 = DimPlot(scRNA, reduction = "tsne") 

3-1

UMAP

scRNA <- RunUMAP(scRNA, dims = pc.num)
embed_umap <- Embeddings(scRNA, 'umap')
plot2 = DimPlot(scRNA, reduction = "umap")

3-2

plotc <- plot1+plot2+ plot_layout(guides = 'collect')
saveRDS(scRNA, file="scRNA.rds")