ixxmu
commented
2 weeks ago 单细胞细节笔记汇总 by 生信漫漫学
写在开头
咱就是说，这个笔记可不得了，汇总了#单细胞细节十六篇文章的内容
所以呢，这个字数呀是比较多的。但是正因为字数多，所以这个笔记非常之详细！也是从系列推文内容的学习中，明白了如果详细解析一篇文章，可以得到这么多的知识点！
果然是单细胞的细节！每一个细节都值得学习和思考！
单细胞转录组降维聚类分群过滤基因和过滤细胞的区别
过滤细胞
首先是创建seurat对象时候的两个参数：
min.cells：
指定了一个基因需要在至少多少个细胞中表达（即检测到的表达量大于0）的最小细胞数。例如，如果设置 min.cells = 5，那么只有那些在至少5个细胞中表达的基因才会被保留在最终的表达矩阵中。这个过滤步骤有助于去除那些可能由于随机噪声或测序误差而在极少数细胞中被检测到的基因，这些基因不太可能对生物学解释有贡献
min.features：
指定了细胞需要表达的最小基因数。例如，如果设置 min.features = 300，那么只有表达超过300个基因的细胞才会被保留在分析中。这个过滤步骤有助于去除那些测序深度不足或质量较低的细胞，这些细胞可能无法提供足够的生物学信息，或者可能代表了技术或生物学上的异常情况。
使用这些过滤参数的目的是从原始数据中移除不可靠或不相关的信息，从而提高后续分析（如聚类、差异表达分析等）的准确性和可解释性。过滤后的数据集应该更加聚焦于生物学上有意义的变异，减少分析结果中的假阳性。
更多细胞过滤参数：
一般比较宽松，不设置过滤的上限，但是有些数据比较特殊或者过滤完发现结果质量不好，可以重新调整一下质控标准
min.features = 300 是规定了细胞的下限，它就是计算得到的nFeature_RNA值，如果这个值过低说明细胞质量肯定是不行，但是过高也不是好事，因为有可能是双细胞。但是如果是处于细胞增殖状态的细胞或者恶性的肿瘤细胞，它的nFeature_RNA值本来就会很高。所以我们很少会设置nFeature_RNA值的上限
3个常见的过滤参数及其意义：
nFeature_RNA：
代表每个细胞中检测到的基因特征（如基因或转录本）的数量。过滤掉低的nFeature_RNA值可以帮助去除那些测序深度不够或背景噪声较高的细胞。
nCount_RNA： 表示细胞中所有基因的测序读段总数。它是一个衡量每个细胞测序数据总量的指标。去除 nCount_RNA过低的细胞可以排除那些可能由于技术原因导致数据不足的细胞。percent_mito： 代表线粒体基因表达占总基因表达的比例。高比例的线粒体基因表达可能指示细胞应激或损伤。通过设置过滤阈值，可以去除可能不健康或死亡的细胞。
tmp=subset(sce.all.filt,subset = percent_mito < 20 &
#nCount_RNA <10000 &
nFeature_RNA <5000 &
nFeature_RNA>200)
percent_mito < 20：过滤掉线粒体基因表达占比超过20%的细胞，这有助于去除可能由于细胞应激或损伤而异常表达线粒体基因的细胞。
nFeature_RNA < 5000：过滤掉表达特征数超过5000的细胞。这可能是为了去除高表达细胞，它们可能代表了技术异常或生物学上的异常状态，如细胞周期的某些阶段。nFeature_RNA > 20：过滤掉表达特征数少于200的细胞，这有助于去除测序深度不足或背景噪声较高的细胞。
过滤基因
过滤基因的操作就是单细胞转录组标准流程里面的降维聚类分群环节的降维，通过筛选高变基因以及PCA降维对基因进行过滤
过滤基因（Filtering Genes）：
主要是设置：min.cells = 5：有助于过滤掉只在极少数细胞中表达的基因，这些基因可能由于随机噪声或技术变异而出现。 根据特定的标准（如表达水平、检测到的细胞数等）去除那些可能不可靠或生物学意义不大的基因。
挑选高变基因（Highly Variable Genes, HVGs）：
高变基因是指在细胞间表达差异较大的基因。这些基因在单细胞数据集中可能特别有信息量，因为它们可能与细胞状态、功能或身份有关。 通常选择一定数量的高变基因进行后续分析，比如2000个。
降维（Dimensionality Reduction）
降维技术用于减少数据集中的变量数量，同时尽可能保留数据的重要结构和变异。 在单细胞转录组分析中，常用的降维方法包括主成分分析（PCA）和t-分布随机邻域嵌入（t-SNE）。
没有绝对正确的单细胞转录组质量控制指标
默认过滤不会设置nFeature_RNA值的上限
第一层次就把髓系免疫细胞区分成为巨噬细胞和中性粒细胞，还有mast细胞，然后呢把b细胞是可以区分成为浆细胞和普通b细胞 很容易看到不同单细胞亚群在不同样品的数量的分布情况：
设置nFeature_RNA < 5000
会删除掉处于细胞增殖状态的细胞或者恶性的肿瘤细胞，它的nFeature_RNA值本来就会很高
两种质控标准对比
两次的结果相差3000多个细胞，其中1000多个是PT2这个样品的上皮细胞被删除，因为是恶性的肿瘤细胞。另外的500多是其它样品的上皮细胞里面的恶性的肿瘤细胞被删除 看到不同的质量控制参数并不会影响大局观，仅仅是被过滤掉的细胞可能会有不同的生物学意义。
学习单细胞亚群命名的层次结构
第一层次降维聚类分群
单细胞亚群主要是肿瘤微环境里面的3大亚群，就是上皮细胞和免疫细胞，以及间质细胞。
免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群 stromal 基质细胞里面的 fibro 和endo进行细分
B细胞细分
B细胞的几个不同阶段和类型的描述：
Naive B细胞（未经激活的B细胞）： 这些是尚未遇到其特定抗原的B细胞。它们在骨髓中成熟，并在没有遇到病原体的情况下循环于血液中。 Memory B细胞（记忆B细胞）： 记忆B细胞是在先前的免疫反应中遇到特定抗原后形成的。它们具有长期存活的能力，并能够快速响应再次出现的相同抗原，提供更快速和有效的免疫反应。 Germinal Center (GC) B细胞： GC B细胞是在淋巴器官的生发中心经历快速增殖和突变的B细胞。它们参与体细胞高频突变（somatic hypermutation, SHM）和抗体类别转换（class-switch recombination, CSR），以产生更高亲和力的抗体。 Plasma Cells（浆细胞）： 浆细胞是B细胞激活和分化的最终阶段，主要功能是大量产生和分泌抗体。它们可以在骨髓中长期存活，并作为长期的抗体来源。
不同B细胞之间的区别：
成熟阶段：Naive B细胞是成熟的早期阶段，而GC B细胞和记忆B细胞是成熟过程中的更晚阶段。 功能：Naive B细胞尚未发挥功能，GC B细胞参与抗体的亲和力成熟，记忆B细胞提供长期免疫记忆，浆细胞负责抗体的产生。 持久性：记忆B细胞和浆细胞在体内有较长的寿命，而Naive B细胞和GC B细胞的寿命相对较短。 分布：Naive B细胞主要在血液和脾脏中循环，GC B细胞主要在淋巴结的生发中心，记忆B细胞和浆细胞可以在骨髓、脾脏和其他组织中找到。 抗原经验：Naive B细胞没有遇到抗原的经验，而记忆B细胞和浆细胞是由先前接触过抗原的B细胞发展而来。
不是造假胜似造假的单细胞降维聚类分群
GSE176078乳腺癌数据细胞分群，以及基于marker基因进行归一化后降维聚类分群
直接降维聚类分群
会得到一个不太好看的umap图，细胞亚群之间好像也没有分的很开
并不需要在乎这个颜值，因为我们对单细胞表达量矩阵做降维聚类分群其实是为了给每个单细胞一个身份而已，只需要身份本身没有错误，或者没有很大的冲突， 就不需要介意umap图上面的斑点。
如果一定要让每个亚群泾渭分明呢
如果我们本来就是要做有监督的分析，比如降维聚类分群后想把不同单细胞亚群泾渭分明的区分开， 那么就可以在数据前期处理做有监督的挑选，比如我们仅仅是挑选那些不同单细胞亚群的特异性高表达量基因去做降维聚类分群。
load('check-by-celltype/qc-_marker_cosg.Rdata')
cg=unique(as.character(unlist(marker_cosg$names)))
input_sce <- ScaleData(input_sce,features = cg)
input_sce <- RunPCA(input_sce, features = cg)
seuratObj <- RunHarmony(input_sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj, dims = 1:15,
reduction = "harmony")
p = DimPlot(seuratObj,reduction = "umap",label=T,group.by = 'celltype' ) ;p
ggsave(filename='umap-by-cosg-genes-by-celltype.pdf',plot = p)
首先你需要对自己的单细胞表达量矩阵进行降维聚类分群和命名，然后依据不同的单细胞亚群进行cosg找各个亚群的特征基因。
上面的上皮细胞会混入一点点到SMC里面，其实这部分上皮细胞就是乳腺basal上皮细胞，它强烈的具有平滑肌的特征
是否需要抹除细胞周期对单细胞降维聚类分群的影响呢
细胞周期背景知识
细胞周期是细胞生长、复制DNA并分裂为两个新细胞的过程。细胞周期主要分为四个阶段：G1期（生长期）、S期（合成期）、G2期（准备期）和M期（分裂期）：
G1期（生长期）： 在G1期，细胞生长并合成蛋白质，为DNA复制做准备。这个阶段是细胞周期中最长的阶段，细胞会评估其环境和内部状态，以决定是否继续进行细胞周期或进入静止状态（G0期）。 S期（合成期）： S期是细胞复制其DNA的阶段。在这个阶段，每个染色体都被精确复制，形成两个相同的染色单体，为后续的细胞分裂做准备。 G2期（准备期）： G2期是细胞在进行分裂前的最后准备阶段。在这个阶段，细胞继续合成蛋白质，并检查DNA复制过程中可能出现的错误，进行修复。 M期（分裂期）： M期，也称为有丝分裂，是细胞实际分裂成两个新细胞的阶段。M期进一步分为前期（Prophase）、中期（Metaphase）、后期（Anaphase）和末期（Telophase），最后是细胞质分裂的胞质分裂（Cytokinesis）。 处于增殖期的状态： 在生物学概念的细胞周期中，只有处于G1期、S期和G2期的细胞才被认为是处于增殖期，因为它们正在准备或进行DNA复制和细胞分裂。 M期是细胞周期的最后阶段，细胞在这个阶段实际分裂成两个新的细胞。 而G0期是细胞周期的一个可选阶段，细胞可以暂时退出细胞周期，进入静止状态，直到接收到适当的信号重新开始细胞周期。 在某些文献中，研究人员可能会将G2期和M期一起讨论，因为它们在细胞准备分裂的过程中是连续的，所以有时会简称为"G2M期"。这种表述方式强调了细胞从G2期过渡到M期的连续性，以及这两个阶段在细胞周期中的重要性。
第一层次降维聚类结果
第一层次降维聚类分群中单独命名了一个细胞增殖亚群，如果看它的top10的特异性基因，是可以比较清晰的看到我们的定义的cycle细胞亚群主要是t细胞和b细胞的混合体。而其他很多癌症单细胞数据，主要是肿瘤细胞会有恶性增殖的现象
细胞周期去除与否的结果
抽样进行细胞周期计算
抽样后细胞聚类还是和原始的一样
细胞周期计算
细胞周期计算靠CellCycleScoring函数：
s.genes <- str_to_upper(cc.genes$s.genes)
g2m.genes <- str_to_upper(cc.genes$g2m.genes)
s.genes;g2m.genes
sce=sce.all.int
sce <- CellCycleScoring(sce,
s.features = s.genes,
g2m.features = g2m.genes,
set.ident = TRUE)
colnames(sce@meta.data)
p1=VlnPlot(sce,features = "S.Score" ,pt.size = 0,group.by = 'celltype') +NoLegend()
p2=VlnPlot(sce,features = "G2M.Score" ,pt.size = 0,group.by = 'celltype')+NoLegend()
p1+p2
table(sce$Phase,sce$celltype)
cycle的亚群很明显主要是t细胞，但是它细胞周期打分就非常高
ScaleData函数回归"S.Score","G2M.Score"
# 回归"S.Score","G2M.Score"
input_sce <- ScaleData(sce, features = rownames(sce),
vars.to.regress = c("S.Score","G2M.Score") )
input_sce<-RunPCA(input_sce )
p1=DimPlot(input_sce,reduction="pca") + DimPlot(sce,reduction="pca")
p2=DimPlot(input_sce,reduction="pca",group.by='celltype') + DimPlot(sce.all.int,reduction="pca") +NoLegend()
p1/p2
seuratObj<-RunHarmony(input_sce,"orig.ident")
names(seuratObj@reductions)
seuratObj<-RunUMAP(seuratObj,dims=1:15,
reduction="harmony")
p=DimPlot(seuratObj,reduction="umap",label=T,group.by='celltype');p
基于ScaleData函数对"S.Score","G2M.Score"回归后，确实是影响了表达量矩阵以及后面的pca分析
但是这个操作并不会把cycle细胞亚群给打散到其它单细胞亚群，如下所示这个cycle细胞亚群仍然是自己聚集在一起
因为这个数据集中定义的cycle细胞亚群主要是t细胞和b细胞的混合体，所以其实这个cycle细胞亚群确实是被打散到其它单细胞亚群，不过是分散到t细胞和b细胞的亚群里面。
每个单细胞亚群取子集后继续降维聚类分群标准操作（以b细胞为例）
提取B细胞细分
提取命名为Bcells和plasma的细胞亚群，提取需要的counts矩阵以及metadata信息重新创建seurat对象，走基本的降维聚类分群流程
rm(list=ls())
options(stringsAsFactors = F)
source('../../scRNA_scripts/lib.R')
set.seed(123456789)
###### step1:导入数据 ######
sce.all.int = readRDS('../../2-harmony/sce.all_int.rds')
colnames(sce.all.int@meta.data)
load('../../phe.Rdata')
table(phe$celltype)
# 因为前面的b细胞和浆细胞被我拆开命名，所以这个时候细分子集需要把两个都提取到
choose_cellIds=rownames(phe)[phe$celltype %in% c('Bcells','plasma') ]
sce.all=sce.all.int[,choose_cellIds]
sce.all=CreateSeuratObject(
counts = sce.all@assays$RNA$counts,
meta.data = sce.all@meta.data
)
取了子集后的单细胞转录组数据对象，仍然是进行了harmony整合哦。因为harmony操作针对的是pca分析结果，而取子集后相当于是一个从零开始的表达量矩阵了没有做过pca也没有任何降维聚类分群操作。也就是说，harmony操作并不会改变表达量矩阵本身。 结果中混入了T细胞亚群
还需要B细胞细分亚群里面的基因列表来区分naive和memory这两种b细胞哦，在免疫学中，B细胞是适应性免疫系统的关键组成部分，负责产生抗体来识别和中和外来病原体。
细分后进行富集以及高级分析
GO和KEGG富集分析
拟时序分析——需调整根节点
细胞通讯分析
上皮细胞里面混入了淋巴系和髓系免疫细胞呢
在选择某个细胞亚群进行细分的时候，会遇到混入别的细胞亚群的情况，比如在选了B细胞亚群细分的时候，混入了少量T细胞
提取上皮细胞进行细分亚群
提取上皮细胞亚群之后，提取之后基于counts和metadata重新创建seurat对象然后走降维聚类分群步骤
仍然是使用了harmony算法去除个体差异，如果是要考虑肿瘤病人间的异质性就不推荐使用harmony哦！ 分群的结果中也是混入了T细胞
进一步确定混入的亚群身份
查看每个亚群的top10marker基因 如果是从上面的标记基因来看，确实是全部的细分亚群都是有上皮细胞特异性的基因的，但是如果看每个亚群的top基因和功能，就可以确定这些是淋巴系和髓系免疫细胞而不是上皮细胞啦。 查看细胞亚群的功能，通过GO和kegg富集分析得到结果
单细胞亚群取子集后的细分亚群再命名的两个难题
基质细胞亚群细分
如果成纤维细胞亚群中混入了平滑肌细胞（Smooth Muscle Cells, SMCs）和周细胞（Pericytes），这个混合的单细胞亚群可以被称为“混合间充质细胞群”（Mixed Mesenchymal Cell Cluster）或“血管相关间充质细胞群”（Vascular-Associated Mesenchymal Cell Cluster），具体命名取决于细胞的来源和所处的组织微环境。
以下是一些可能的命名方式：
成纤维细胞/平滑肌细胞/周细胞亚群（Fibroblast/Smooth Muscle Cell/Pericyte Subpopulation）：
这种命名直接指出了亚群中包含的细胞类型。
由于周细胞和平滑肌细胞都与血管结构相关，这个名称强调了它们的血管周围定位。
如果这些细胞具有干细胞特性，这个名称可能适用。
如果这些细胞是从血管壁中分离出来的，这个名称可能更准确。
如果这些细胞来自特定的组织或器官，可以在名称中指定组织来源。
如果这些细胞具有多种功能或分化潜能，这个名称可以反映它们的多功能性。
合理的命名的第一个难题：亚群应该是细分到什么程度
我们默认使用的是RNA_snn_res.0.1，很简单的提高一下分辨率到RNA_snn_res.0.5，就可以看到b细胞4大亚群其实是可以裂变的
在张泽民老师的单细胞文章：《Pan-cancer single-cell dissection reveals phenotypically distinct B cell subtypes》就是针对这些更细致的b细胞亚群进行了解释
如果没有生物学背景，仅仅是靠每个人在自己的数据集里面去定义，会出现各种千奇百怪的亚群，而且大多数都很难复现出来。
合理的命名的第二个难题：没有统一的生物学背景
算不上什么大错误的成纤维细胞亚群的细分操作
胃癌文章中对成纤维细胞细分结果
胃癌单细胞数据集GSE163558对应的文献“Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing”，**作者首先是默认第一层次降维聚类分群里面的成纤维细胞是纯粹的成纤维。然后就细分成为了 iCAF和myCAF **
但是如果成纤维里面的混入了平滑肌细胞（Smooth Muscle Cells, SMCs）和周细胞（Pericytes）并不能说它们是mCAF的： 附件有对他们细分成为了 iCAF和myCAF的特异性的go功能注释，都是有ECM，很明显就是iCAF不纯粹导致的：
其他文章中对成纤维细胞细分
“Classifying cancer-associated fibroblasts-The good, the bad, and the target”的文章，通过单细胞表型和CAF在不同癌症类型中的空间位置来简化和标准化CAF的分类。分类包括以下的9个caf类型：
基质CAF（mCAFs）：参与细胞外基质的产生和维护。 炎症CAF（iCAFs）：在肿瘤微环境中促进炎症反应。 血管相关CAF（vCAFs）：与血管相关，参与血管生成。 分裂CAF（dCAFs）：处于增殖状态并积极分裂。 肿瘤样CAF（tCAFs）：表现出类似肿瘤细胞的特征。 热休克蛋白高表达的肿瘤样CAF（hsp_tCAFs）：一种特殊类型的tCAF，表达高水平的热休克蛋白。 干扰素反应CAF（ifnCAFs）：参与干扰素反应途径。 抗原呈递CAF（ap-CAFs）：具有呈递抗原的能力，可能影响免疫反应。 网状CAF（rCAFs）：在肿瘤内形成网络，通常与结构支持相关。
这一分类系统为理解CAF的功能多样性及其作为癌症治疗潜在靶点提供了框架。通过基于表型和位置定义清晰的类别，研究人员可以更好地研究CAF在不同癌症背景下的作用，并开发调节其活性的策略。
成纤维细胞背景知识：
正常组织的成纤维细胞和肿瘤相关成纤维细胞（CAF），是两种在组织结构和功能上具有重要作用的细胞类型，成纤维细胞在正常组织中发挥结构和修复作用，而CAF在肿瘤微环境中具有更复杂的功能，包括促进肿瘤生长和调节免疫反应。
成纤维细胞（Fibroblasts）：
形态：长梭形的间质细胞。 历史：由Rudolf Virchow于1858年首次描述，Ernst Ziegler于1881年命名。 功能：在所有组织中负责产生、维护和重塑细胞外基质（ECM）。 静止状态：在健康组织中通常处于静止状态。 激活：在伤口愈合和炎症过程中被激活，伴随着基因表达的变化，如ACTA2和FAP的上调，可以命名为肌成纤维细胞。 肌成纤维细胞：是激活后的成纤维细胞，分泌促炎细胞因子，促进伤口愈合和组织修复。 凋亡：在成功愈合后，通常会经历凋亡，以恢复组织稳态。
相似性：与肌成纤维细胞有许多相似之处。 肿瘤微环境（TME）：在TME中扮演关键角色，参与肿瘤的发生、发展、转移和治疗抵抗。 功能：分泌生长因子、细胞因子和细胞外基质组分，影响肿瘤细胞行为。 免疫调节：通过与免疫细胞的相互作用，调节肿瘤微环境中的免疫反应。 双重作用：在TME中具有促进和抑制肿瘤生长的双重作用，被称为“双刃剑”。 治疗靶点：CAF的功能多样性使它们成为癌症治疗的潜在靶点。 异质性：CAF的异质性给研究和临床应用带来挑战。 分类：依赖于单细胞RNA测序等技术，但缺乏统一的分类标准，导致不同研究中对CAF的命名和定义存在差异。
单细胞亚群绝对数量和相对比例的探索
批量组合不同变量（分组，样品）看细胞类型比例
cal_table = function(x,y,prefix ){
# x = phe$orig.ident
# y = phe$celltype
library(sur)
library(reshape2)
tbl = table(x,y)
pdf(paste0(prefix,'-table.pdf'),width = 10,height = 10)
gplots::balloonplot( tbl )
dev.off()
df = dcast(as.data.frame(tbl),x~y)
head(df)
write.csv( df ,paste0(prefix,'-table.csv'))
# ptb = round(sur::percent.table(x,y),2)
ptb = round(100*tbl/rowSums(df[,-1]),2)
pdf(paste0(prefix,'-percent-table.pdf'),width = 10,height = 10)
gplots::balloonplot( ptb )
dev.off()
write.csv( dcast(as.data.frame(ptb),x~y) ,paste0(prefix,'-percent-table.csv'))
}
cal_table(phe$orig.ident,phe$celltype,prefix = 'celltype-vs-orig.ident')
cal_table(phe$group,phe$celltype,prefix = 'celltype-vs-group')
cal_table(phe$group,phe$seurat_clusters,prefix = 'seurat_clusters-vs-group')
cal_table(phe$group,phe$RNA_snn_res.0.8,prefix = 'RNA_snn_res.0.8-vs-group')
直接查看绝对数量：可以看到的是t细胞占比非常夸张，而中性粒细胞和b细胞呢，有两个组的特异性高表达这两个细胞亚群
查看对应的比例：应该是在具体的某个分组内部，或者在某个样品内部计算不同单细胞亚群的比例。而不是在某个单细胞亚群内部去计算它来源于不同组或者不同样品的比例
多种可视化方式
折线图
x='celltype';y='group'
plot_data <- data.frame(table(phe[, y ],
phe[, x ]))
head(plot_data)
plot_data$Total <- apply(plot_data,1,function(x)sum(plot_data[plot_data$Var1 == x[1],3]))
plot_data <- plot_data %>% mutate(Percentage = round(Freq/Total,3) * 100)
colnames(plot_data) <- c("group","cluster","Freq","Total","Percentage")
head(plot_data)
ggplot(plot_data,aes(group,Percentage,group= cluster,color =cluster))+geom_point(size=4)+
geom_line(position = position_dodge(0.1),cex=2)+theme_bw()+theme_test(base_size = 30)
ggsave("celltype-vs-group-percent-plot.pdf",width = 9,height = 7)
带比例的柱状图
有监督的挑选了特征之后的无监督的分析还可靠吗
层次聚类和主成分分析
层次聚类（Hierarchical Clustering）和主成分分析（PCA）都是无监督学习方法，它们可以用来探索样品之间的关系，而不需要预先定义的类别或标签 层次聚类更多地用于探索样品的分类和分组，而PCA则常用于数据的可视化和预处理，以便于进一步分析。 在实际应用中，这两种方法可以结合使用，先通过PCA降维以减少计算复杂性，然后使用层次聚类来探索样品之间的关系。 但是这个操作仍然是无偏倚的，比如过滤那些在所有的样品都不表达的基因或者都低表达量基因，或者按照表达量的sd或者mad排序后选择top的基因即可
层次聚类（Hierarchical Clustering）：
层次聚类是一种聚类算法，通过构建一个多层次的嵌套聚类树来组织数据。 它不需要预先指定聚类的数量，而是生成一个聚类树（树状图或dendrogram），可以从中选择不同的切割点来得到不同数量的聚类。 层次聚类可以是凝聚的（从单个样本开始，逐步合并聚类）或分裂的（从所有样本作为一个聚类开始，逐步分割）。
主成分分析（PCA）：
PCA是一种降维技术，通过线性变换将数据投影到较低维度的空间，同时尽可能保留原始数据的变异性。 它通过找到数据中方差最大的方向（主成分），并将数据投影到这些方向上，从而揭示样品之间的结构关系。 PCA通常用于数据可视化和预处理，帮助识别数据中的模式和结构。
同时量化不同样品的不同单细胞亚群的相关性
sce.all = readRDS('./2-harmony/sce.all_int.rds')
sp='human'
load('./phe.Rdata')
identical(rownames(phe) , colnames(sce.all))
sce.all@meta.data = phe
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident)
DimPlot(sce.all)
# v4 是 AverageExpression
av <-AggregateExpression(sce.all ,
group.by = c("orig.ident","celltype"),
assays = "RNA")
av=as.data.frame(av[[1]])
head(av) # 可以看到是整数矩阵
#write.csv(av,file = 'AverageExpression-0.8.csv')
cg=names(tail(sort(apply(log(av+1), 1, sd)),1000))
df =cor(as.matrix(log(av[cg,]+1)))
ac=as.data.frame(str_split(colnames(df),'_',simplify = T))
rownames(ac)=colnames(df)
pheatmap::pheatmap(df ,
show_colnames = F,show_rownames = F,
annotation_col = ac,
filename = 'cor_celltype-vs-orig.ident.pdf' )
dev.off()
save(av,file = 'av_for_pseudobulks.Rdata')
不同样品的同一个单细胞亚群还是会优先聚集到一起的 不同的样品虽然是有生物学来源差异，不足以让一个t细胞去跟上皮细胞更加相似的。除非它们都处于某种生物学状态，比如细胞增殖，或者说它们都有双细胞的特性等等。
到底是量化样品还是单细胞之间的相关性呢
GSE206785胃癌单细胞数据集
单细胞之间的相关性
不同单细胞亚群之间的这个差异仍然是最大的，其次是各个亚群里面的样品差异或者分组差异
样品直接的相关性
如果忽略了不同样品其实是由不同的单细胞亚群组成的， 那么它们这些样品虽然说来源于两个很明显的分组，癌症组织和正常组织，它们其实并不会泾渭分明
为什么胃癌并不使用拷贝数来判断恶性的肿瘤上皮细胞呢
恶性肿瘤上皮细胞的组学特性：
拷贝数变化（CNVs）：许多恶性肿瘤细胞存在基因组的拷贝数变化，这可能导致癌基因的扩增或肿瘤抑制基因的丢失。 其他组学改变：除了拷贝数变化，恶性肿瘤细胞还可能表现出其他类型的基因组变异，例如单核苷酸变异（Single Nucleotide Variants, SNVs）、插入和缺失（Indels）、基因重排和易位等。 表观遗传调控：表观遗传学改变，如DNA甲基化、组蛋白修饰和非编码RNA的表达变化，也在肿瘤发生和发展中起重要作用。这些改变可能影响基因表达而不一定涉及DNA序列的改变。 转录组学改变：肿瘤细胞的转录组学分析可以揭示基因表达模式的改变，这些改变可能与肿瘤的增殖、侵袭、转移和药物抗性有关。 蛋白质组学改变：蛋白质组学分析可以揭示肿瘤细胞中蛋白质表达水平、翻译后修饰和蛋白质-蛋白质相互作用的变化。 代谢组学改变：肿瘤细胞可能表现出代谢途径的改变，以支持其快速增殖和生存。 肿瘤微环境：肿瘤细胞与其微环境的相互作用，包括免疫细胞、成纤维细胞、血管细胞和细胞外基质，也对肿瘤的恶性表型有重要影响。 肿瘤异质性：肿瘤内部的细胞可能表现出不同的基因组和表型特征，这种异质性可能影响肿瘤对治疗的响应。 罕见肿瘤类型：某些罕见的肿瘤类型可能不表现出显著的拷贝数变化，而是通过其他机制获得恶性表型。 肿瘤发展阶段：在肿瘤发展的早期阶段，拷贝数变化可能较少，随着肿瘤进展，基因组不稳定性可能增加。
拷贝数变化是许多恶性肿瘤的共同特征，但它们并不是唯一的分子机制。肿瘤细胞可能通过多种机制获得恶性表型，包括基因组、转录组、蛋白质组和代谢组层面的变化，以及表观遗传调控和肿瘤微环境的影响。这些复杂的生物学过程共同促进了肿瘤的发生、发展和异质性。
单细胞转录组水平看拷贝数是最佳实践
目前的最佳实践是依据单细胞转录组表达量矩阵来推断染色体级别的拷贝数变化情况，一般来说， 典型的恶性肿瘤上皮细胞会出现大面积的染色体扩增和缺失的。
基于单细胞进行恶性细胞的鉴定方法有3种：——恶性细胞的鉴定方法
簇状分布：非恶性的分布常按细胞类型聚集在一起，恶性细胞的分布常按肿瘤的样本类型聚集在一起。
CNV：恶性细胞会有异常的拷贝数变异。应用 inferCNV 以T细胞为参考对照，推断每个细胞中的大规模CNV。下图中的颜色越红表示拷贝数越高，越蓝表示拷贝数越低
标记基因表达：在恶性细胞中，恶性细胞的标记基因会高表达，正常的谱系特异性基因低表达或无表达。
胃癌文章里面的恶性细胞判断方式
胃癌单细胞数据集GSE163558的单细胞转录组在降维聚类分群后，针对里面的上皮细胞进行恶性推断的时候，就放弃了拷贝数方法。使用的是类似于标记基因表达的策略
inferCNV辅助判断
首先是从这个单细胞表达量数据集里面挑选'Tcells','Bcells'这样的免疫细胞，我们默认上皮细胞癌症组织里面的免疫细胞仍然是没有拷贝数变异的，可以作为正常的二倍体参考信息： 很明显的看到上皮细胞的拷贝数确实是有变化，但是居然中性粒细胞等也有拷贝数：
恶性细胞评分
走inferCNV流程的时候只需要针对上皮细胞即可
一些总结
对于胃癌，它跟其它上皮细胞来源的癌症有一点点不一样的地方是，它的上皮细胞相对来说没有那么的恶性，以至于走inferCNV流程的时候它甚至没办法在拷贝数变异程度这个指标去跟免疫细胞或者基质细胞进行区分！ 有一个很简单的改进方法， 就是不要试图去对整个表达量矩阵进行走inferCNV流程。因为大多数癌症里面的恶性的肿瘤细胞确实起源于上皮细胞，这类癌症被称为上皮源性癌症或 carcinomas。 上皮细胞是覆盖在身体表面（包括皮肤）和体内腔道（如消化道、呼吸道、泌尿生殖道）的细胞，它们也构成了器官的内衬。由于上皮组织在身体中的广泛分布，上皮源性癌症是最常见的癌症类型。
上皮源性癌症的例子包括但不限于：
肺癌 乳腺癌（虽然乳腺癌发生在乳腺组织，但乳腺也是由上皮细胞组成的） 胃癌 结直肠癌 前列腺癌 宫颈癌 皮肤癌（如黑色素瘤）
其他类型的癌症：
肉瘤（Sarcomas）：起源于间充质细胞，包括肌肉、骨骼、软骨、血管、脂肪等。 淋巴瘤和白血病：起源于血液和淋巴系统。 中枢神经系统肿瘤：起源于脑和脊髓的细胞。 神经内分泌肿瘤：起源于神经内分泌系统的细胞。
基于胃癌数据中上皮细胞进行inferCNV分析
输入数据的格式要求：
主要是表达量矩阵，其中行名是基因名字，如果是人类，就需要根据这些基因名字去获取它们的染色体和坐标。 列表是细胞，有属性，比如来源于不同样品的细胞，或者细胞本身可以是不同的生物学名字！
step1:制作reference_mat.Rdata
需要指定好两种不同的非上皮细胞的单细胞亚群作为正常的二倍体参考，这个可以任意选择哈。 上面的案例选择了成纤维细胞和浆细胞，其实大部分情况下我比较喜欢选择内皮细胞和肥大细胞，但是有一些情况下这两个细胞亚群的数量是少于800的，为了不修改后面的代码，我就懒得调整里面的变量名字了。
step2:根据三种格式的文件来运行myInferCNV函数
人为的选择好的两种不同的非上皮细胞的单细胞亚群，确实是正常的二倍体，上面没有明显的染色体级别的拷贝数变异。 上皮细胞来源于每个不同样品，它们代表了肿瘤病人的个体异质性，所以拷贝数变化情况都不一样！ 上面的染色体级别的拷贝数扩增和缺失热图，如果仅仅是肉眼看，肯定是过于主观了，比如上面的NT1样品来源的上皮细胞粗看起来确实是没什么拷贝数变异。
step3: 解读infercnv结果
重新载入上皮细胞的亚群，对亚群进行打分和可视化
可以比较清晰的看到我们的定义的上皮细胞里面混入了淋巴系和髓系免疫细胞，以及 正常个体来源的上皮细胞，这3个群都是拷贝数相对来说比较低的 而标记为epi系列的细胞亚群，拷贝数变异较高
虽然对每个上皮细胞都是进行了inferCNV打分，但实际上我们并没有使用这个inferCNV打分来探索最佳阈值去划分细胞的恶性与否！
我们是先分组，然后看不同组之间的inferCNV打分的差异情况来划分正常和对照的组，并不是针对具体的每个单细胞在进行判断！
癌症和癌旁的差异基因能在单细胞层面区分上皮细胞的恶性与否吗
癌症和癌旁的差异基因综合打分
胃癌文章中去tcga数据库里面定位到胃癌的转录组测序数据集，然后根据分组做癌症和癌旁的差异分析后，拿到上下调各自的top50基因列表去单细胞表达量矩阵里面打分！
文献的附件里面有这些基因的名字，所以可以基于这些基因名字去我们的单细胞转录组表达量矩阵里面进行打分
sce.all.int = readRDS('2-harmony/sce.all_int.rds')
sce.all.int
f='phe.Rdata'
load(f)
colnames(phe)
sce.all.int=sce.all.int[,colnames(sce.all.int) %in% rownames(phe)]
identical(colnames(sce.all.int), rownames(phe))
sce.all.int$celltype = phe$celltype
th = theme(axis.text.x=element_text(angle=45,hjust = 1))
table(sce.all.int$celltype)
DimPlot(sce.all.int,group.by = 'celltype',
repel = T,label = T)
# 文章里面的 tcga的转录组差异top50基因
up='IBSP CST4 HOXC9 FEZF1 HOXA11 RDH8 ALPP HOXC10 MAGEA12 TAS2R38 EVX1 CCL7 C5orf46 DMBX1 HOXC13 FGF19 MAGEA6 HOXC12 R3HDML ALPG MMP8 MAGEA3 HOXC11 PRAC2 SP8 DCSTAMP CST1 PIWIL1 KRTAP4-1 CSF2 CSAG3 ACTL8 HOXC8 BAAT'
up=trimws(strsplit(up,' ')[[1]]);up
down='SLC5A7 PCSK2 DCAF12L1 SLC7A3 GCG CHIA RPRM PLP1 PTF1A KIAA0408 NR0B1 PEBP4 CCKBR LGI3 HSPB3 PLD5 CCKAR SYT4 NUPR2'
down=trimws(strsplit(down,' ')[[1]]);down
使用Seurat::AddModuleScore函数进行打分：
gene_vector=list(up=up,down=down)
sc_dataset <- Seurat::AddModuleScore(sce.all.int,
features = gene_vector)
#signature.names <- paste0(names(gene_vector), "_UCell")
options(repr.plot.width=6, repr.plot.height=4)
colnames(sc_dataset@meta.data)
table(sc_dataset$orig.ident)
Idents(sc_dataset)
p1=VlnPlot(sc_dataset, features = 'Cluster1',
group.by = "orig.ident",pt.size = 0 ) + NoLegend() ;p1
FeaturePlot(sc_dataset,'Cluster1')
p2=VlnPlot(sc_dataset, features = 'Cluster2',
group.by = "orig.ident",pt.size = 0 ) + NoLegend() ;p2
FeaturePlot(sc_dataset,'Cluster2')
p1+p2
可以比较清晰的看到了文献里面的上调基因在所有的癌症样品里面的打分都比较高，但是在NT1这个正常样品里面就打分很低 （下图的上半部分小提琴图）,文献里面的下调基因在NT1这个正常样品里面就打分很高（上图的下半部分小提琴图）。
使用了上调基因打分减去下调基因的打分，取一个 0.02的阈值，来划分正常细胞和恶性细胞。但是这个打分也是很难完美的区分正常细胞和恶性细胞，因为那个NT1 正常样品里面的上皮细胞理论上都是正常的，被这个算法误判了大部分。
根据文献内容进行调整
文献里面的做法是去tcga数据库里面定位到胃癌的转录组测序数据集，然后根据分组做癌症和癌旁的差异分析后，拿到上下调各自的top50基因列表去单细胞表达量矩阵里面打分！ 每次差异分析上下调基因非常多，不能说是仅仅是靠变化倍数极端值来挑选基因，应该各自独立的差异分析后使用rra算法去汇总：
上面的这些在多个癌症和癌旁差异分析都稳定的表现出来上下调的基因，就可以比较好的去单细胞层面区分上皮细胞的恶性与否
bulk层面的癌症和癌旁的表达量差异主要是因为？
癌症和癌旁组织的表达量芯片的差异分析，上下调基因可能的来源
通过差异分析，研究者可以识别在癌症组织中特异性表达的基因，这些基因可能在肿瘤的发生、发展、转移和预后中发挥重要作用。进一步的研究可以探索这些差异表达基因的功能，为癌症的诊断和治疗提供潜在的靶点。
https://mp.weixin.qq.com/s/48aVJv41p32WnGnNEInazQ