前言

我之前写过这篇文章的精读笔记：
单细胞好文2--Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in PDAC
接着跟着Seurat官网学了两个教程：
盲人摸象--single cell sequencing的Seurat学习--详细版
趁热打铁--scRNA cell-type specific responses

本事没多少，胆子倒挺肥，就想试一试水，所以今天就挑文章：Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in PDAC 的数据来试着完成基础部分。先放一张图来馋一下

image.png

1. 准备工作

1.1 数据下载

老规矩，在文章中找到数据存放的地方：

data available

意外发现一个类似于GEO的数据库NGDC；谷歌搜索GSA: CRA001160

image.png

1.2 载入数据并构建Seurat对象

由于数据非常大，文件count-matrix.txt就有2.6 Gb，如此大的文件我的小小电脑运行起来非常吃力，这里提醒需要使用服务器的R来运行数据。

（1）载入数据，生成矩阵

library(Seurat)
library(data.table)

df <- fread("count-matrix.txt",  sep = " ")
?fread
# 快读读取数据

df[1:10,1:4]  # 简单查看数据

mt <- as.matrix(df[,-1]) # 去除第一列
row.names(mt) <- df$V1 # 改列名
rm(df); gc()   # 及时清除不需要的数据和垃圾
?gc    # Garbage Collection，不懂就要打问号好好学习

（2）创建Seurat对象
Seurat有自己的数据存放形式，由于一个小型数据库，构建的方式是使用函数CreateSeuratObject

seurat.obj <- CreateSeuratObject(mt, min.features = 3, names.delim = "_")
# 一个细胞至少有3个基因？文章里面这样写的呀

image.png

2. QC

2.1 细胞过滤即QC条件设置

在前面的推文中已经提出QC需要注意的3个点：
（1）一个基因能在多少个细胞中被检测到（min.cells）;
（2）一个细胞能检测到多少基因（min.features）;
（3）线粒体基因的比例要足够低（percent.mt）。
在这里我们根据文章描述内容进行QC的条件设置

# 添加线粒体基因的比例,加到meta.data里面的percent.mt
seurat.obj[["percent.mt"]] <- PercentageFeatureSet(seurat.obj, pattern = "^MT-")


# 筛选cells and features，以下是文章中的标准
# Low quality cells (<200 genes/cell, <3 cells/ gene and >10% mitochondrial genes) were excluded
seurat.obj <- subset(seurat.obj, subset = nFeature_RNA > 200  & percent.mt < 10)

## 小结：我们的标准可以更高一下，线粒体基因的百分比下调到5，min.cells = 3, min.features = 200（这部分在导入数据的时候就要界定）

2.2 可视化QC结果

喜闻乐见的检查方法（文章中没有这部分）

（1）小提琴图

# 可视化QC
VlnPlot(seurat.obj, 
        features = c("nFeature_RNA", 
                     "nCount_RNA", 
                     "percent.mt"), 
        ncol = 1)  

image.png

通过查看QC plot我们知道，这些数据的MT比例真的很高，如果以5为标准，将会丢失非常多的数据。

（2）散点图查看特征之间的关系
我尝试使用之前的代码来实现这部分内容

plot1 <- FeatureScatter(seurat.obj,  feature1 = "nCount_RNA",  feature2 = "percent.mt")

plot2 <- FeatureScatter(seurat.obj,  feature1 = "nCount_RNA",  feature2 = "nFeature_RNA")
plot1 + plot2

image.png

出来的结果让我想发笑，R语言太“南”了，个人不知是什么原因导致这样的报错，但不要为难自己，分开画plot1和plot2也是可以的。

image.png

这个图其实就是看一眼而已，好的测序数据不管怎么看都完美的。

3. Normalizing the data 标准化数据

节选我的学习笔记：默认情况下，使用全局缩放规范化方法LogNormalize，该方法通过总表达式对每个单元格的特征表达式度量进行标准化，并将其乘以一个缩放因子(默认为10,000)，然后对结果进行log转换。标准化值存储在pbmc[["RNA"]]@data中。

seurat.obj <- NormalizeData(seurat.obj, 
                            normalization.method = "LogNormalize", scale.factor = 10000)

# 方法是LogNormalize
# scale.factor是它默认的

4. Identification of highly variable features (feature selection) 寻找高差异基因

使用FindVariableFeatures完成差异分析，选择数据集中差异较高的特征基因（默认2000）并用于下游分析。

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

小小的标记一下

# 标记出top 10
top10 <- head(VariableFeatures(seurat.obj), 10)
# 绘图
plot1 <- VariableFeaturePlot(seurat.obj)
plot2 <- LabelPoints(plot = plot1,  points = top10, repel = TRUE)
plot1 + plot2

image.png

5. Scaling the data 缩放数据

5.1 缩放的标准

应用线性变换缩放数据，这是一个标准的预处理步骤，比PCA等降维技术更重要。使用ScaleData函数:
（1）使每个基因在所有细胞间的表达量均值为0
（2）使每个基因在所有细胞间的表达量方差为1

缩放操作给予下游分析同等的权重，这样高表达基因就不会占主导地位，其结果存储在pbmc[["RNA"]]@scale.data中

# 全局缩放
all.genes <- rownames(seurat.obj)
seurat.obj <- ScaleData(seurat.obj, features = all.genes)

5.2 删除不需要的数据

使用```ScaleData``函数从单细胞数据集中删除不需要的变量。例如，我们可以“regress out”与细胞周期阶段或线粒体污染相关的异质性（去除不需要的变量，即校正协变量，校正线粒体基因比例的影响）。如下：

# 删除不需要的数据，校正线粒体基因引起的影响
# 这个步骤非常耗时，连512G的服务器都要跑好久（超过5分钟）
seurat.obj <- ScaleData(seurat.obj, vars.to.regress = "percent.mt")

6. PCA 主成分分析

这个分析非常重要，选择合适的PC数，对后面的分析是关键性的。

seurat.obj <- RunPCA(seurat.obj, features = VariableFeatures(object = seurat.obj))

image.png

PCA可视化有以下方法：
（1）DimPlot

# DimPlot
DimPlot(seurat.obj, reduction = "pca")

image.png

（2）DimHeatmap

DimHeatmap(seurat.obj, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(seurat.obj, dims = 1:20, cells = 500, balanced = TRUE, ncol = 5)

image.png

7. Determine the ‘dimensionality’ of the dataset 计算数据维度

为了克服scRNA-seq数据的technical noise，Seurat基于它们的PCA分数对单元进行分组，每个PC实质上代表一个“metafeature”，它将跨相关特征集的信息组合在一起。但是到底我们该选多少个PC呢？即，我们需要确定数据集的维度，这部分是真正的重点。

（1）ElbowPlot

ElbowPlot(seurat.obj)

image.png

可见，前20个PC数都是可以选择的，看个人的标准。

（2）JackStraw
慎用

# 对于大的数据集，JackStraw会耗费更多时间
seurat.obj <- JackStraw(seurat.obj, num.replicate = 100) # 重复一百次
seurat.obj <- ScoreJackStraw(seurat.obj, dims = 1:20) # 选择20个PC

# 可视化每个PC的P value分布
JackStrawPlot(seurat.obj, dims = 1:20)

image.png

能有多耗费时间？原先跟着教程做的时候，2700个cells耗时3 min 49 s（使用的是我的8G电脑），而换了这个512G内存的服务器，这边竟然也跑了33 m 07 s。

image.png

上图的结果让我傻眼了，竟然有如此多的PC的value是0，只有PC10有确切的pvalue，而且也是相当低，目前还没有解释。但从拐点图的结果看，前20个PC有可能还没到达平缓的拐点，因此我们建议将拐点图的PC数设置多一些：

ElbowPlot(seurat.obj,ndims = 50, reduction = "pca") 
# 最大可设置到50

image.png

将拐点图的PC数设置到50（最大值）后可视化，可见我的猜想是正确的，前20个PC确实还没到拐点，甚至严格来说前30个OC都不算，所以在这里，PC数的选择将考验个人对数据对课题的理解。从数据分析的角度，如果是我，我愿意选择严格的做法，即PC40；而真正情况下是要都做一下，看看哪一个数据维度数的选择对后续分析的效果是最好的。

8. Cluster the cells 细胞聚类

步骤：
（1）先在PCA空间中构造一个基于欧氏距离的KNN图，然后根据任意两个细胞在局部邻近区域的共享重叠(Jaccard相似性)来细化它们之间的边权值。此步骤使用FindNeighbors函数执行，并将之前定义的数据集维度作为输入。
（2）使用FindClusters函数对数据进行优化，并包含一个分辨率参数，该参数设置下游集群的“granularity”，增加的值将导致更多的集群。将该参数设置在0.4-1.2之间，对于3K左右的单细胞数据集通常会得到良好的结果。对于较大的数据集，最佳分辨率通常会增加。可以使用Idents函数找到clusters。

seurat.obj <- FindNeighbors(seurat.obj, dims = 1:20) # 先选择20个
seurat.obj <- FindClusters(seurat.obj, resolution = 0.8) # 越高越精细
head(Idents(seurat.obj), 3)

image.png

如果将resolution的值调低，找到的cluster数也会相应降低。

image.png

可以看到UMAP和tSNE图略有不同，UMAP皱缩成团分不开，而tSNE图则分散得多。

9. Run non-linear dimensional reduction (UMAP/tSNE) 非线性降维

（1）非线性降维可达到数据可视化，常用的UMAP/tSNE的方法。

seurat.obj <- RunUMAP(seurat.obj, dims = 1:20)
seurat.obj <- RunTSNE(seurat.obj, dims = 1:20)
DimPlot(seurat.obj, reduction = "umap",label = T)
DimPlot(seurat.obj, reduction = "tsne",label = T)

image.png

（2）标记cluster
标记上cluster的名字，以下只是示例，并不是本次分析的真正结果。实际上我还需要找到每个cluster的marker，然后定义每个cluster的细胞类型，之后生产相应label，加到Seurat对象中，最后可视化出来。这里假定已经找好marker并且已经定义好细胞，下面即是可视化。为什么我没有自行找marker后定义细胞再label呢？因为这个工作量比较大，并且需要有非常好的生物学背景知识，中间涉及的已经不再是代码问题了。

# 准备label相关数据（假如已经定义好每一个cluster的细胞）
# 注意以下只是为了示范，用了该文章提供的数据
# 我们在实际分析过程中，并没有找到文章的参数情况，因此无法完全炮制
df <-fread("all_celltype.txt")
# 将cluster的名字放进去
seurat.obj[['labels']] <- df$cluster

DimPlot(seurat.obj,  label = FALSE)
DimPlot(seurat.obj, group.by = "labels", label = TRUE, reduction = "tsne")

image.png

（3）根据细胞来源给线性降维图上色

p1 <- DimPlot(seurat.obj, reduction = "tsne", group.by = "orig.ident")
p2 <- DimPlot(seurat.obj, reduction = "umap", group.by = "orig.ident")
p1
p2

image.png

10. Finding differentially expressed features (cluster biomarkers) 差异基因

虽然上面已经假定细胞类型已经定义完了，但这个过程还是需要继续学习的，目前首先要找到每个cluster的marker，之后按照marker或者GO分析定义cluster是什么细胞（依据GO分析定义细胞是在汤富酬教授组的文章看到的）。使用函数FindMarker完成找marker的过程：
(1) min.pct参数要求在两组单元中至少检测到一个特征；
(2) thresh.test参数要求一个特性在两组之间有一定的差异(平均)；
若两个值都设置为0，会非常耗时——因为这将测试大量不太可能具有高度差异性的features。
(3) 为了加速，可以设置max.cells.per.ident参数，这将向下采样每个标识类，使其单元格数不超过设置的单元格数。

举例分析如下
（1）找cluster 5的marker

cluster5.markers <- FindMarkers(seurat.obj, 
                                ident.1 = 1,  # 至少检测到1个feature
                                min.pct = 0.25) # 设置后会加速，默认是0.1

head(cluster5.markers, n = 5) # 显示！

image.png

（2）找某个cluster与其他cluster的差异基因（marker）

# 找cluster 5的markers
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(seurat.obj, 
                                ident.1 = 5, # 找cluster5的marker
                                ident.2 = c(0, 3), # 和cluster 0-3对比
                                min.pct = 0.25)
head(cluster5.markers, n = 5)

image.png

（3）找所有cluster的marker
首先，简单粗暴的找所有cluster的marker，先定义了细胞类型，而后再去查看自己感兴趣的细胞类型。挑出感兴趣的cluster，互相之间找差异基因（marker），这个操作话了我相当久的时间，35个cluster大约花了3-4小时去做Findmarker。

# report only the positive ones
seurat.obj.markers <- FindAllMarkers(seurat.obj, only.pos = TRUE, 
                               min.pct = 0.25, logfc.threshold = 0.25)

seurat.obj.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

11. 可视化cluster marker

（1）小提琴图

# 小提琴图将marker可视化，随便挑上部分结果中的3个基因
# 以下3个基因在cluster0中高表达
VlnPlot(seurat.obj, features = c("DARC", "CLDN5","SPRY1"),ncol = 1)

图片太大就不放了

（2）降维图显示

FeaturePlot(seurat.obj, features = c("DARC", "CLDN5","SPRY1","RGS5"))

image.png

暂时先做到这里啦~

下次见。