跟洛克菲勒的老师再跑一次Seurat流程

入门是跟着B站人气最高的Biomamba兄视频入门的，保姆级教程，而且提前把要踩的坑都给你想到。

手把手带教感觉一切都很顺利，到自己的数据一下就懵了……有个研三快毕业的群友发了很长的感想，总结一下就是看了三年的视频教程，感觉啥都会了，但是啥都写不出来，最后还是靠纯临床分析毕业。这是学习的常态，总要有个机械记忆和慢慢理解的过程，不能总靠教练扶着的。要保持输出，记在小本本上也好，私密的公开的都可以，关键还得记下来，不要过于相信记忆。我就是半年前看了视频觉得自己行了，确实跑了好几次流程也都没问题，但隔了一段时间又啥都忘了，特别是拿到那些千奇百怪的文件，愁啊。

跟着洛克菲勒生信中心重新学习Seurat流程，他们的网页在：

Session2.knit

下载数据

这次拿到的是10x cellranger处理后的标准三个文件，GSE147415 数据集做了两组，但是示例中用的仅仅是其中的对照组。
没有他们Dropbox网盘的权限，所以用的数据集比他们的要大一个数量级，出图效果也有差别。

构建Seurat对象

以下内容引自Seurat对象结构：https://www.jianshu.com/p/238976158dcc

简书外链似乎不够友善：https://www.bilibili.com/video/BV1Qf4y1s7h1

用的比较多的

Seurat Object

@meta.data

对所有细胞做注释的数据框。每一行代表一个细胞，每一列代表细胞的属性。当需要根据细胞的属性和类型对细胞进行筛选的时候，经常会用到mata.data。当然也可以把分析得到的结果，添加到mata.data中。

@assay

一个Seurat对象可以包括多个assay对象，但是在某个时刻，只有一个assay对象是默认激活的。

@ident

可以理解为细胞的类型，在Seurat对象中，细胞可能有好几种不同方法注释的类型，但是在某一时刻，只有一种细胞类型是默认激活的。

@reduction 和assay一样，reduction返回的也是一个列表。里面包含的是一个或多个 DimReduc object 对象。

Assay Object

三个matrix：

@counts 保存未经处理的原始数据，适合存放稀疏矩阵

@data 存放经过标准化后的数据

@scale.data 存放scale后的数据

@var.features 向量，存放高表达变异的基因名

@meta.features 对每个 features 做的注释

DimReduc Object

每个DimReduc object 对象对应的是 assay 对象进行某种降维分析后得到的结果。降维也就是PCA 、tsne 、umap 三种。

标准流程

library(Seurat)
####创建Seurat对象及生成QC图####
## 读入CellRanger的三个文件，创建对象
tar_dir <- "./"
samID <- "CTRL"
X10_file <- Seurat::Read10X(tar_dir)
obj <- CreateSeuratObject(counts = X10_file, project = samID, min.cells = 5, min.features = 200)
obj$dset <- samID
obj <- Seurat::RenameCells(obj, add.cell.id = samID)
## 估计线粒体含量
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^mt-")
## 对象元数据包含细胞barcode、线粒体含量、数据集和聚类等
head(obj@meta.data)
## 原始计数矩阵
head(obj@assays$RNA@counts)
## 缩放数据或归一化数据(没做出来)
head(obj@assays$RNA@data)
## 评估数据
## - 计数：nCount_RNA - 基因数：nFeature_RNA - 线粒体含量： percent.mt
VlnPlot(obj, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"), pt.size = 0.2)
## 比较计数和基因数
FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
## 同样的 比较计数和线粒体含量
FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "percent.mt")
## 去除细胞碎片：过滤掉高线粒体含量和低读取计数的
summary(obj@meta.data$percent.mt)
mt_cutH <- 10
obj_unfiltered <- obj
obj <- subset(obj, subset = percent.mt < mt_cutH)
## 过滤前后比较
VlnPlot(obj_unfiltered, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),
        pt.size = 0.2)
VlnPlot(obj, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"), pt.size = 0.2)
FeatureScatter(obj_unfiltered, feature1 = "nCount_RNA", feature2 = "percent.mt")

####细胞周期测定####
##在Seurat，S期和G2/M期的默认标记基因被储存在cc.gene list
##基于人的list
##小鼠可以通过bioMart转换
##报了502错误 搜了一下是biomaRt的问题，那只能换一个HOST
## 转换基因列表
library(biomaRt)
convertHumanGeneList <- function(x){
  require("biomaRt")
  human = useMart("ensembl", dataset = "hsapiens_gene_ensembl",host = "<https://dec2021.archive.ensembl.org>")
  mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl",host = "<https://dec2021.archive.ensembl.org>")
  genesV2 = getLDS(attributes = c("hgnc_symbol"), 
                   filters = "hgnc_symbol", 
                   values = x , 
                   mart = human,
                   attributesL = c("mgi_symbol"), 
                   martL = mouse, uniqueRows=T)
  humanx <- unique(genesV2[, 2])
  return(humanx)}
#
s_gene <- convertHumanGeneList(cc.genes$s.genes)
g2m_gene <- convertHumanGeneList(cc.genes$g2m.genes)
## 评估周期
obj <- CellCycleScoring(obj, s.features = s_gene, g2m.features = g2m_gene, set.ident = TRUE)
obj@meta.data[1:2, ]
yd_dat <- as.data.frame(table(obj@meta.data$dset, obj@meta.data$Phase))
head(yd_dat)
library(ggplot2)
ggplot(yd_dat, aes(x = Var1, y = Freq, fill = Var2)) + geom_bar(stat = "identity",
                                                                position = "stack") + labs(x = "", y = "Counts", fill = "Phase") + theme_classic()
####回归（降维）与聚类####
## 数据归一化
obj <- ScaleData(obj, vars.to.regress = c("percent.mt", "S.score", "G2M.score", "Phase"))
## 主成分分析PCA
set.seed(1000)
obj <- FindVariableFeatures(obj)
obj <- RunPCA(obj, npcs = 30, verbose = FALSE)
## 如果未指定基因列表，默认会用高变基因，前提是已经找了高变基因，不然会报错
## 所以加了一行找高变基因
## 画肘图，估计要采取的PC数量很有用。
## 找肘部，即曲线变平的起始点。这个近似值可以使偏差和PC数量最小化。
ElbowPlot(obj, ndims = 30)
numPC <- 15
## 聚类
set.seed(1000)
maxPC <- numPC
obj <- FindNeighbors(obj, reduction = "pca", dims = 1:maxPC)
obj <- FindClusters(obj, resolution = 0.5)
obj <- RunUMAP(obj, reduction = "pca", dims = 1:maxPC)
obj <- RunTSNE(obj, reduction = "pca", dims = 1:maxPC)
DimPlot(obj, reduction = "umap")  # default it 'umap'
DimPlot(obj, reduction = "tsne")
####注释分群####
obj <- SetIdent(obj, value = "seurat_clusters")
clust.markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(clust.markers)
##select the top 5 marker genes for each cluster
topG <- clust.markers %>%
  group_by(cluster) %>%
  top_n(n = 5, wt = avg_log2FC)
head(topG)
####进阶绘图####
## 调节点的大小
DimPlot(obj, pt.size = 0.2)
## 加上分群标签
DimPlot(obj, pt.size = 0.2, label = TRUE) + NoLegend()
## 做热图看marker基因表达
DoHeatmap(obj, features = topG$gene) + NoLegend()
## Feature plot - 降维可视化marker基因表达
gene_marker <- c("Krt1", "Pthlh", "Krt14", "Cenpa", "Shh")
FeaturePlot(obj, features = gene_marker, pt.size = 0.2)
## 山脊图看marker基因表达
RidgePlot(obj, features = gene_marker)
## 按分群划分的细胞周期阶段
tbl <- table(obj@meta.data$seurat_clusters, obj@meta.data$Phase)
tbl_dat <- as.data.frame(tbl)
to <- rowSums(tbl)
names(to) <- rownames(tbl)
tbl_dat$to <- to[match(names(to), tbl_dat$Var1)]
tbl_dat$prop <- tbl_dat$Freq/tbl_dat$to
tbl_dat[1:2, ]
## 过滤掉第2群，然后保存数据对象
cellID <- rownames(obj@meta.data)[!obj@meta.data$seurat_clusters == 1]
obj_sub <- obj[, cellID]
saveRDS(obj_sub, "scSeq_Seurat_clean.rds")

练习题

就是用KO组再跑一次上面的流程。