search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
118 stars 30 forks source link

harmony单细胞整合分析细节真让人纠结：数据集合并后取变异最大的3000个基因（即不取交集的方法）与取交集方法之间的优劣 #5691

Closed ixxmu closed 2 weeks ago

ixxmu commented 2 weeks ago

https://mp.weixin.qq.com/s/M8q1WZ8dgcWYwVWcT-iYAA

ixxmu commented 2 weeks ago

harmony单细胞整合分析细节真让人纠结：数据集合并后取变异最大的3000个基因（即不取交集的方法）与取交集方法之间的优劣 by 生信小博士


01单细胞整合分析的目的

不同的单细胞整合方式都是为了最后的FindNeighbors、FindClusters()两步，得到聚类信息，进行细胞命名，然后单细胞整合的任务就完成了。（其他的所有分析都是基于分好的亚群基础上进行的。和整合之后获取的data矩阵和scale.data矩阵关系不大。这句话不太对，容易引起误会）

    注意：SCTransform不能去除批次效应    visium_slide = visium_slide %>%      PercentageFeatureSet(pattern = "^Mt-", col.name = "percent.mt") %>%      PercentageFeatureSet(pattern = "^Hb-", col.name = "percent.hb") %>%      PercentageFeatureSet(pattern = "^RP[LS]", col.name = "percent.ribo") %>%
      SCTransform(assay = "RNA",vars.to.regress = c("percent.mt","percent.ribo","percent.hb")) %>%      RunPCA() %>%     RunUMAP(dims = 1:30) %>%      FindNeighbors(dims = 1:30) %>%  FindClusters()


开头一段话种需要澄清和扩展的地方：

1. 整合分析的目标

单细胞RNA测序的整合分析主要目的是对来自不同实验条件、时间点或技术平台的细胞进行整合，使得这些细胞可以在一个统一的空间中进行比较。整合分析的两个主要目标是：

  • 消除批次效应: 通过整合分析减少不同实验条件或平台之间的技术性差异（批次效应），使得来自不同来源的相似细胞可以更好地对齐。

  • 揭示细胞亚群: 通过聚类分析识别不同的细胞亚群，并为这些亚群进行命名，赋予生物学意义。


2. 整合后的数据矩阵（data矩阵和scale.data矩阵）

整合之后得到的data矩阵和scale.data矩阵，虽然不直接用于后续的某些分析（如差异表达分析等），但它们在整合过程和后续的特定类型的分析中仍然具有重要作用：

  • data矩阵: 通常表示的是归一化后的表达值，这些值可以用于后续的基因表达模式分析、差异表达分析等。

  • scale.data矩阵: 在整合分析中，scale.data矩阵用于PCA和聚类分析。这些分析帮助你识别和定义细胞亚群。因此，整合后的scale.data矩阵和聚类信息之间存在直接关系。

3. 聚类和细胞命名

确实，聚类分析和细胞亚群命名是整合分析的关键步骤之一，这为后续的所有分析奠定了基础。所有其他的生物学分析，如差异表达分析、功能富集分析等，都是基于这些聚类信息进行的。因此，准确的聚类和命名对于后续分析至关重要。

4. 整合分析后的用途

虽然整合后的数据矩阵和scale.data矩阵不直接用于某些后续分析，但它们仍然在以下几个方面具有重要作用：

  • 维度缩减和可视化: UMAP、t-SNE等方法依赖于整合后的数据矩阵或其衍生物，用于在低维空间中可视化细胞分布。

  • 整合数据的一致性验证: 整合后的矩阵可以用于评估整合效果，验证不同数据集之间的细胞是否有效对齐，是否消除了批次效应。

  • 整合模型的生成: 这些矩阵可以作为生成整合模型的基础，这些模型可以用来整合新的数据集，或在不同的分析中重复使用。

5. 总结

  • 聚类信息和细胞命名: 是整合分析的核心目标，后续分析确实依赖于这一信息。

  • 整合后的数据矩阵: 虽然在某些分析中不直接使用，但它们在整合过程、可视化、模型生成等方面起到重要作用。



02如何让单细胞整合分析的umap图更好看

如果仅仅为了自己的umap图好看，甚至可以自己选择或者指定高边基因：不是造假胜似造假的单细胞降维聚类分群

在单细胞转录组数据整合分析中，使用特定的基因集（如gene_selection）作为PCA的输入，而不是直接使用所有基因进行PCA，有几个原因：

1. 数据降维和噪声处理:

  • 单细胞转录组数据通常包含大量的噪声和冗余信息。如果直接使用所有基因进行PCA，可能会引入不相关或噪声信号，影响PCA的效果。通过选择高变异基因（HVGs）或特定的基因集（如gene_selection），可以集中在那些在细胞类型之间差异最大的基因，从而更好地揭示细胞间的异质性。

2. 提高计算效率:

  • 使用所有基因进行PCA计算量巨大，尤其是在大规模数据集中。通过预先选择少量具有生物学意义的基因，减少数据的维度，可以显著提高PCA的计算效率，并减少内存使用。

3. 保留生物学意义:

  • 选取的基因集往往包含那些与特定生物学过程或细胞类型相关的基因，这样做能够更好地保留这些生物学信息，从而在后续分析中（如聚类、UMAP、Harmony整合）中，获得更可靠的结果。

4. 增强整合效果:

  • 在多样本或多批次数据整合中，使用特定的基因集进行PCA，并结合Harmony等批次效应去除工具，有助于增强整合的效果，使得来自不同样本的细胞类型能够更好地对齐，而不受批次效应的干扰。


在代码片段中，通过FindVariableFeatures选取3000个变异最大的基因，形成gene_selection，并在PCA时使用这些基因作为输入。这样做有助于将焦点集中在数据中最具信息量的部分，从而提高分析的效果和效率。

相比之下，直接使用所有基因进行PCA，可能会导致计算冗余、噪声干扰，并降低整合的效果，这也是为什么要选择特定基因集作为PCA输入的主要原因。


03 整合分析前，是否需要对所有数据集取交集


1. 不取交集的方法（合并后选择整体变异最大的3000个基因）

优点:

  • 综合考虑所有数据集的变异性: 通过将所有数据集合并后再选择变异最大的3000个基因，这种方法能够综合考虑所有样本的变异性，从而确保所选择的基因在整个数据集中具有高变异性。这可以捕捉到所有数据集之间共享的主要生物学信号。

  • 最大化信息利用: 这种方法利用了整个数据集的完整信息，确保在选择高变异基因时没有遗漏任何在个别样本中表现显著的基因。

  • 简化数据处理流程: 只需对合并后的数据集进行一次高变异基因选择，减少了在不同数据集上重复选择的工作量，同时确保所有数据集的一致性。

缺点:

  • 可能偏向大样本量数据集: 在合并数据集中，大样本量的数据集可能会对基因变异性的计算产生更大影响，导致选择的基因更偏向于这些数据集，从而忽略了小样本量数据集中可能同样重要的基因。

  • 潜在的批次效应影响: 如果不同数据集之间存在显著的批次效应，这种方法可能会导致批次效应在选择高变异基因时的作用被放大，从而在整合分析中引入偏差。


2. 取交集的方法（从每个数据集各自的高变异基因中取交集）

优点:

  • 减少批次效应的影响: 通过从每个数据集的高变异基因中取交集，确保只使用在所有数据集中都表现出高变异的基因，从而减少了批次效应对整合分析的影响。

  • 提高数据一致性: 使用的基因集合在所有数据集中一致，这有助于在整合分析中获得更加稳健和一致的结果，特别是在对齐细胞类型和去除批次效应时。

缺点:

  • 信息丢失: 取交集方法可能会排除那些在某些数据集中重要但在其他数据集中没有表现出显著变异的基因，导致潜在的生物学信息丢失，尤其是在数据集之间存在生物学差异的情况下。

  • 基因集大小受限: 交集往往会导致基因集合的缩小，特别是在数据集之间差异较大时。这可能会降低PCA和后续分析的维度，进而影响分析的全面性和深度。

总结

  • 不取交集: 适用于希望最大限度地利用所有数据集的整体变异信息时，但需要注意可能引入批次效应和偏向大样本量数据集的风险。


  • 取交集: 适用于希望在多数据集整合中减少批次效应并提高分析结果一致性时，但可能会丢失一些特定数据集中的重要基因信息。

最终，选择哪种方法取决于你对批次效应的敏感度、数据集的规模和多样性，以及对全面生物学信息保留的需求。

#取交集进行整合分析的完整代码All.merge=SCTransform(All.merge, verbose = FALSE)  %>% RunPCA(npcs = 30, verbose = FALSE)library('harmony')All.merge <- All.merge  %>% RunHarmony("stim", plot_convergence = TRUE)harmony_embeddings <- Embeddings(All.merge , 'harmony') #######################clusterdims = 1:30All.merge  <- All.merge  %>%   RunUMAP(reduction = "harmony", dims = dims) %>%   RunTSNE(reduction = "harmony", dims = dims) %>%   FindNeighbors(reduction = "harmony", dims = dims)                       使用gene_selection作为pca的输入的代码 ：#####安装archr包##别处复制.libPaths(c("/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2",            "/home/data/t040413/R/yll/usr/local/lib/R/site-library",             "/usr/local/lib/R/library",            "/home/data/refdir/Rlib/"))
.libPaths()library(tidyverse)library(Seurat)library(optparse)library(scater)library(harmony)library(dplyr)library(optparse)library(tidyverse)library(Seurat)library(harmony)library(cluster)library(clustree)#https://github.com/saezlab/visium_heart/blob/5b30c7e497e06688a8448afd8d069d2fa70ebcd2/st_snRNAseq/02_snuc_integration_harmony/integrate_objects.R#L194



slide_files <- list.files("~/silicosis/spatial/prcessed_visum_for_progeny/data/",recursive = TRUE,                                        all.files = TRUE,full.names = TRUE,pattern = "_")



integrated_data <- map(slide_files, readRDS)print("You managed to load everything")print("Object size")print(object.size(integrated_data))




# Calculate HVG per sample - Here we assume that batch and patient effects aren't as strong# since cell-types and niches should be greater than the number of batches
Assays(integrated_data [[1]])str(integrated_data [[1]])


def_assay="Spatial"

hvg_list <- map(integrated_data, function(x) {    DefaultAssay(x) <- def_assay    x <- FindVariableFeatures(x, selection.method = "vst",                             nfeatures = 3000)    x@assays[[def_assay]]@var.features  }) %>% unlist()
hvg_list <- table(hvg_list) %>%  sort(decreasing = TRUE)



gene_selection_plt <- hvg_list %>% enframe() %>%   group_by(value) %>%   mutate(value = as.numeric(value)) %>%  summarize(ngenes = length(name)) %>%   ggplot(aes(x = value, y = ngenes)) +   geom_bar(stat = "identity")
gene_selection <- hvg_list[1:3000] %>% names()

########## Create merged object ---------------------------------integrated_data <- purrr::reduce(integrated_data,                          merge,                          merge.data = TRUE)
print("You managed to merge everything")print("Object size")print(object.size(integrated_data))
# Default assay ---------------------------------------#DefaultAssay(integrated_data) <- def_assay


# Process it before integration -----------------------integrated_data <- integrated_data %>%  ScaleData(verbose = FALSE) %>%   RunPCA(features = gene_selection,          npcs = 30,          verbose = FALSE) 
original_pca_plt <- DimPlot(object = integrated_data,                             reduction = "pca",                             pt.size = .1,                             group.by = "orig.ident")
head(integrated_data@meta.data)
batch_covars="orig.ident"# Integrate the data -----------------------integrated_data <- RunHarmony(integrated_data,                               group.by.vars =   batch_covars,                               plot_convergence = TRUE,                              assay.use = def_assay,                              max.iter.harmony = 20)