Seurat数据处理

之前的学习中，我们现在已经知道如何搞一个scRNA-seq的表达矩阵了，那我们现在就从读取数据开始，进行一系列的分析，最后可视化分析结果。

读取数据

之前我们讲过 filtered_feature_bc_matrix 目录，在这里不多做赘述，直接一个加载滴干活

data.list <- list()
for (i in 1:5) {
  data.dir <- paste0("dataset", i, "/filtered_feature_bc_matrix/")
  data.list[[i]] <- Read10X(dataDir = data.dir)
}

初始化一个空列表data.list，用于存储不同数据集的Seurat对象。
使用for循环遍历1到5（假设有5个数据集），构造每个数据集的路径，并使用Read10X函数读取数据，然后将结果（Seurat对象）存储到data.list中。

创建Seurat对象并进行标准化

seurat.objects <- list()
for (i in 1:length(data.list)) {
  seurat.objects[[i]] <- CreateSeuratObject(counts = data.list[[i]], project = paste0("Sample", i))
}

初始化一个空列表seurat.objects，用于存储处理后的Seurat对象。
遍历data.list中的每个元素（每个数据集），使用CreateSeuratObject创建Seurat对象，并指定项目名称。每个项目名称通过将”Sample”和索引i连接而成，确保每个Seurat对象有一个唯一的项目名称。

整合数据

FindIntegrationAnchors和IntegrateData函数是Seurat包提供的数据整合工具，通过寻找跨数据集的相似细胞群体（锚点）并利用这些信息来校正数据，从而减少批次效应，使得来自不同实验条件或平台的数据能够在细胞水平上进行有效比较。

integration.anchors <- FindIntegrationAnchors(object.list = seurat.objects, dims = 1:30)
integrated.data <- IntegrateData(anchorset = integration.anchors, dims = 1:30)

FindIntegrationAnchors函数在多个数据集之间寻找锚点，这些锚点用于后续的数据整合。dims = 1:30表示使用前30个主成分进行锚点的寻找。
IntegrateData函数使用上一步找到的锚点来整合数据，生成一个新的Seurat对象integrated.data。

质量控制

subset是Seurat包中用于筛选Seurat对象中的细胞（观测值）和/或特征（基因）的函数。这个函数允许基于各种条件对数据进行子集化，非常有用于去除低质量的细胞、选择特定表达模式的细胞，或是根据其他标准（如实验分组）进行筛选。

## 1. 筛选细胞
# 只保留基因数量在200到2500之间且线粒体基因表达比例小于5%的细胞
integrated.data <- subset(integrated.data, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mito < 0.05)

## 2. 筛选特征（没明白，有案例可以单开一篇）
# 计算特定基因的表达水平
gene.expression <- FetchData(seurat.data, vars = "GeneX")
# 筛选基于特定基因表达阈值的细胞
selected.cells <- WhichCells(seurat.data, expression = GeneX > threshold)
# 使用selected.cells来筛选Seurat对象
seurat.obj.filtered <- subset(seurat.data, cells = selected.cells)

## 3. 筛选特定身份的细胞
# 如果只对特定的聚类结果或身份感兴趣，仅保留聚类标识为1和2的细胞
seurat.obj <- subset(seurat.data, idents = c(1, 2))

使用subset函数根据基因数量和线粒体基因表达比例筛选细胞。主要参数:

x: 待筛选的Seurat对象。
subset: 用于筛选的逻辑条件。可以基于细胞的元数据（如基因表达量、线粒体基因比例等）或特征表达设置条件。
cells: 一个可选的细胞标识符向量，用于进一步限制筛选的范围。
features: 一个可选的特征（基因）标识符向量，用于筛选特定的基因。
idents: 用于筛选的一个或多个身份类别（例如聚类结果），只有这些类别中的细胞会被保留。

标准化、查找高变异基因、缩放数据、进行PCA

在单细胞RNA测序数据分析中，标准化、查找高变异基因、数据缩放和主成分分析（PCA）是关键的预处理步骤:

## 1. 标准化
# 标准化是调整不同细胞之间测序深度差异的过程。它将原始的计数数据转换为可比较的表达水平，常见的方法包括对每个细胞的基因表达计数进行对数转换。
# normalization.method = "LogNormalize"指定了标准化方法，这里使用对数归一化。
# scale.factor = 10000是一个比例因子，用于在对数转换之前将每个细胞的基因表达计数缩放至平均每细胞10,000个转录本。
integrated.data <- NormalizeData(integrated.data, normalization.method = "LogNormalize", scale.factor = 10000)

## 2.  查找高变异基因
# 在单细胞数据中，部分基因在细胞间表达具有显著变异性，这些基因往往与细胞状态或类型有关。查找这些高变异基因对于后续分析至关重要
# selection.method = "vst"选择使用方差稳定转换（VST）方法查找高变异基因。
# nfeatures = 2000指定选取变异性最高的2000个基因。
integrated.data <- FindVariableFeatures(integrated.data, selection.method = "vst", nfeatures = 2000)

## 3. 数据缩放: 
# 数据缩放（或标准化）是指减去特征（如基因表达）的均值并除以其标准差，使得数据在不同特征间具有可比性，这对于PCA等降维方法尤为重要
# features = rownames(seurat.obj)指所有基因
integrated.data <- ScaleData(integrated.data , features = rownames(integrated.data))

## 4. 主成分分析（PCA）
# PCA是一种降维技术，通过保留数据中的主要变异来源来简化数据结构，它能够揭示样本间的关键差异和模式。
# 使用前面找到的高变异基因运行PCA。VariableFeatures(object = seurat.obj)返回Seurat对象中已标记为高变异的基因列表。
integrated.data <- RunPCA(integrated.data, features = VariableFeatures(object = integrated.data))

NormalizeData:对数据进行标准化。
FindVariableFeatures:找到最具变异性的2000个基因，用于后续的分析。
ScaleData:对数据进行缩放，为降维分析做准备。
RunPCA:进行主成分分析（PCA），为数据降维提供基础。

查找邻居和聚类

在Seurat中，查找邻居（FindNeighbors）和聚类（FindClusters）是单细胞RNA-seq数据分析的关键步骤，用于探索和定义细胞的自然群体（即聚类），这些群体通常对应于不同的细胞类型或细胞状态。

查找邻居（FindNeighbors）

FindNeighbors函数基于高维数据（如PCA或其他降维结果）计算每个细胞的局部邻居。这是通过构建一个k近邻图（k-NN图）完成的，该图中的节点代表细胞，如果两个细胞是彼此的最近邻，则它们之间存在一条边。

# 使用Seurat对象中PCA分析的前20个维度来寻找每个细胞的邻居
integrated.data <- FindNeighbors(integrated.data, dims = 1:20)

object: Seurat对象。
dims: 用于邻居搜索的维度（通常是PCA的主成分）。例如，如果使用PCA结果，dims可以设置为前几个主要成分，如1:10或1:20。
k.param: k近邻的数量，默认值通常为20，这意味着每个细胞将寻找其20个最近的邻居。

聚类（FindClusters）


integrated.data <- FindClusters(integrated.data, resolution = 0.5)

object: Seurat对象。
resolution: 控制聚类粒度的参数。较高的值将导致更多的聚类，较低的值导致较少的聚类。这个参数需要根据数据集进行调整。
algorithm: 选择使用的聚类算法，默认为1（Louvain算法）。

进行UMAP降维和可视化

UMAP（Uniform Manifold Approximation and Projection）是一种流行的降维技术，特别适用于大规模数据集，如单细胞RNA测序（scRNA-seq）数据。

UMAP能够有效地揭示高维数据中的结构，例如细胞类型的聚类，同时保持全局数据结构的同时揭示局部群体的细微差异。Seurat包通过RunUMAP函数提供了UMAP降维的实现，并可以通过DimPlot等函数进行可视化

integrated.data <- RunUMAP(integrated.data, dims = 1:20)

RunUMAP 函数接受一个Seurat对象，并基于之前的PCA结果进行UMAP降维。
dims = 1:20 参数指定使用PCA结果的前20个主成分进行UMAP计算。这个参数的选择取决于PCA分析结果和数据集的特性。

完成UMAP降维后，你可以使用DimPlot函数来可视化降维的结果。

DimPlot(integrated.data, reduction = "umap", label = TRUE)

reduction = "umap"指定使用UMAP的结果进行绘图。
group.by = "ident"参数表示按照Seurat对象中存储的细胞身份（ident，通常是聚类结果）进行分组，为每个组分配不同的颜色。
label = TRUE添加每个群组的标签。
pt.size = 0.5调整点的大小。

找到差异表达基因并保存结果

在单细胞RNA测序（scRNA-seq）分析中，寻找差异表达基因（DEGs）是理解细胞间异质性和识别细胞类型的关键步骤。Seurat提供了FindAllMarkers函数来寻找差异表达基因，这个函数能够在不同的细胞群体（通常是聚类结果）之间比较基因表达，以识别每个群体的标志性基因。

FindAllMarkers函数

object: Seurat对象。
only.pos: 是否只返回正向标志性基因（即在特定群体中上调的基因）。
min.pct: 要求基因在至少一组细胞中的表达百分比达到该阈值。
logfc.threshold: 对数折叠变化（log fold change, logFC）阈值，用于定义显著差异表达。
test.use: 使用哪种统计测试进行DEGs分析，常用的有wilcox（Wilcoxon rank sum test）。

# 寻找差异表达基因
integrated.data <- FindAllMarkers(object = integrated.data, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, test.use = 'wilcox')
# 保存差异表达基因结果
write.csv(integrated.data@markers, file = "DEG_clusters.csv")

only.pos = TRUE表示只关心在特定聚类中上调的基因。
min.pct = 0.25要求基因至少在25%的某一聚类细胞中被检测到。
logfc.threshold = 0.25设置对数折叠变化阈值为0.25，意味着差异表达基因的表达量至少需要在一组中比另一组高出25%。
test.use = 'wilcox'指定使用Wilcoxon rank sum test进行差异表达分析。

这个流程涵盖了从数据读取、整合、质量控制，到降维聚类和差异表达基因分析的完整过程。每一步都是为了确保数据的质量，并从多个数据集中提取有意义的生物信息。

bajiu