癌症和癌旁的差异基因能在单细胞层面区分上皮细胞的恶性与否吗

癌症和癌旁的差异基因能在单细胞层面区分上皮细胞的恶性与否吗 by 生信技能树

前面在这个笔记：为什么胃癌并不使用拷贝数来判断恶性的肿瘤上皮细胞呢，我提到了对于胃癌来说，它跟其它上皮细胞来源的癌症有一点点不一样的地方是，它的上皮细胞相对于免疫细胞或者间质细胞来说没有那么的恶性。所以我们推荐走inferCNV流程的时候只需要针对上皮细胞即可，这样的话可以比较清晰的看到我们的定义的上皮细胞里面混入了淋巴系和髓系免疫细胞，以及正常个体来源的上皮细胞，这3个群都是拷贝数相对来说比较低的。

也就是说，来源于normal样品的上皮细胞以及少量来源于其它癌症样品的上皮细胞，组合成为了一个可能的normal上皮细胞亚群，而它又恰好是被我们的inferCNV流程推断成为了拷贝数比较低的亚群：

虽然我们对每个上皮细胞都是进行了inferCNV打分，但实际上我们并没有使用这个inferCNV打分来探索最佳阈值去划分细胞的恶性与否！我们是先分组，然后看不同组之间的inferCNV打分的差异情况来划分正常和对照的组，并不是针对具体的每个单细胞在进行判断！

判断恶性与否的另外一个方法：癌症和癌旁的差异基因综合打分

大家可以详细的阅读这个胃癌单细胞数据集GSE163558对应的文献“Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing”，研究者写的蛮清楚这个方法学，就是去tcga数据库里面定位到胃癌的转录组测序数据集，然后根据分组做癌症和癌旁的差异分析后，拿到上下调各自的top50基因列表去单细胞表达量矩阵里面打分！

文献的附件里面有这些基因的名字，所以很容易去我们的单细胞转录组表达量矩阵里面进行打分，前面的走inferCNV流程的时候只需要针对上皮细胞即可和上皮细胞里面混入了淋巴系和髓系免疫细胞呢做好了上皮细胞的细分：

rm(list=ls())
options(stringsAsFactors = F) 
getwd()
source('../../scRNA_scripts/lib.R') 
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的阈值，来划分正常细胞和恶性细胞。

ss=sc_dataset$Cluster1 -sc_dataset$Cluster2;head(ss)
plot(sort(ss))
x=1:length(ss)
y=sort(ss);  
plot(x, y, type = "l", col = "black")
fivenum(ss)
table(ss < -0.02)  
table(ss < -0.02,sce.all.int$orig.ident)

可以看到绝大部分细胞的这个打分都是在0 附近，所以为什么取0.02也是未解之谜：

> table(ss < 0.02) 

FALSE  TRUE 
  283  2044 
> table(ss < -0.02) 

FALSE  TRUE 
 2258    69 
> table(ss < -0.02,sce.all.int$orig.ident) 
       
         Li1  Li2  LN1  LN2  NT1   O1   P1  PT1  PT2  PT3
  FALSE   10   21    2   56   86  291  133  181 1442   36
  TRUE     0    0    0    0   44    1    0    1   14    9

而且这个打分也是很难完美的区分正常细胞和恶性细胞，因为那个NT1 正常样品里面的上皮细胞理论上都是正常的，被这个算法误判了大部分。

但是文章里面秀出来的数据就很完美，可以看到那个NT1 正常样品里面的上皮细胞绝大部分都被判定为正常：

其实是因为上下调基因不应该是直接简单粗暴挑选top50或者top100

文献里面的做法是去tcga数据库里面定位到胃癌的转录组测序数据集，然后根据分组做癌症和癌旁的差异分析后，拿到上下调各自的top50基因列表去单细胞表达量矩阵里面打分！

但实际上，每次差异分析上下调基因非常多，不能说是仅仅是靠变化倍数极端值来挑选基因，我个人比较推崇的一个方法学是RRA算法，如下所示做多个不同国家地区不同芯片平台的数据集，各自独立的差异分析后使用rra算法去汇总：

上面的这些在多个癌症和癌旁差异分析都稳定的表现出来上下调的基因，就可以比较好的去单细胞层面区分上皮细胞的恶性与否啦。

 #   https://doi.org/10.3389/fgene.2018.00265
# Identification of Potential Key Genes Associated With the Pathogenesis and Prognosis of Gastric Cancer Based on Integrated Bioinformatics Analysis
up= c("CST1", "SFRP4", "THY1", "SULF1", "FAP", "SPP1", "CEACAM6", "TIMP1", "CLDN1", 
      "MMP7", "THBS2", "BGN", "ASPN", "CLDN7", "CDH17", "PLA2G7", "OLFM4", "CLDN4", 
      "MFAP2", "SPARC")
down <- c( "GIF", "GKN1", "ATP4A", "PSCA", "KCNE2", "LIPF", "ATP4B", 
           "SST", "CHIA", "GKN2", "CA9", "CPA2", "TFF2", "FBP2", "PGC", "SOSTDC1", "CA2", 
           "ANXA10", "ALDH3A1");

gene_vector=c(up,down)
# remotes::install_github("carmonalab/UCell")
library(UCell)
scRNA=sce.all.int
scRNA
gene_vector=list(up=up,down=down)
gene_vector
names(gene_vector)
colnames(scRNA@meta.data)
sc_dataset <- Seurat::AddModuleScore(scRNA, 
                                     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

ss=sc_dataset$Cluster1 -sc_dataset$Cluster2;head(ss)
plot(sort(ss))
x=1:length(ss)
y=sort(ss);  
plot(x, y, type = "l", col = "black")
fivenum(ss)
table(ss < -0.25)  
table(ss < -0.25,sce.all.int$orig.ident)

可以看到，这个时候的阈值就不是文章里面的0.02这么诡异的情况了，我们的阈值来判断上皮细胞的恶性与否就比较符合常识，因为绝大部分正常样品里面的上皮细胞都会被我们的打分划分为正常上皮细胞。

> table(ss < -0.25,sce.all.int$orig.ident) 
       
         Li1  Li2  LN1  LN2  NT1   O1   P1  PT1  PT2  PT3
  FALSE    5   20    2   56   20  289  133  177 1451   38
  TRUE     5    1    0    0  110    3    0    5    5    7

当然了，这里面仍然是有一个问题，就是为什么我们选择的是0.25这个值呢，实际上是一个数学问题，就是判断下面的曲线的拐点：

写在文末

如果你也想做单细胞转录组数据分析，最好是有自己的计算机资源哦，比如我们的2024的共享服务器交个朋友福利价仍然是800，而且还需要有基本的生物信息学基础，也可以看看我们的生物信息学马拉松授课（买一得五），你的生物信息学入门课。

ixxmu / mp_duty