【原】单细胞亚群比例变化和表达量差异分析

健明 2022-05-15

展开全文

单细胞的多组的设计（比如正常组与给药组）可以为细胞类型水平比较提供以往Bulk RNA-seq分析所不能达到的精度。对此一般有两种进阶分析思路：

（1）DE(Differential expression)--两组样本的同一细胞类型的基因表达差异分析；
（2）DA(Differential abundance)--两组样本的同一细胞类型的丰度差异分析

我们以Nov 2020的文献：《VEGF-B Promotes Endocardium-Derived Coronary Vessel Development and Cardiac Regeneration》为例，链接是：https://www./doi/10.1161/CIRCULATIONAHA.120.050635

文章里面的单细胞戏份并不多，但是很容易理解，适合初学者入门。首先是单细胞亚群比例，如下所示：

两个分组的单细胞亚群比例

可以看到，AAV9-Vegfb组相对于 AAV9-S2来说，主要是处于增殖状态的内皮细胞这个单细胞亚群比例增多了。

降维聚类分群，参考前面的例子即可：人人都能学会的单细胞聚类分群注释，它数据在：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE128509 ，可以看到是4个样品:

GSM3678221 aMHC-Vegfb p7 TG
GSM3678222 aMHC-Vegfb p7 WT
GSM3678223 AAV9-Vegfb
GSM3678224 AAV9-S2

每个样品都给出来了标准的3个文件，如下所示：

GSM3678221_aMHC-Vegfb_TG_barcodes.tsv.gz 11.1 Kb
GSM3678223_AAV9-Vegfb_barcodes.tsv.gz 11.9 Kb
GSM3678223_AAV9-Vegfb_genes.tsv.gz 212.7 Kb
GSM3678223_AAV9-Vegfb_matrix.mtx.gz 14.5 Mb
GSM3678224_AAV9-S2_barcodes.tsv.gz 14.0 Kb
GSM3678224_AAV9-S2_genes.tsv.gz 212.7 Kb
GSM3678224_AAV9-S2_matrix.mtx.gz 12.9 Mb

这样的文件需要首先改名，然后读入，走标准流程，大概是十几分钟就可以做完全部的分析。

首先是文件批量改名；

这个 https://ftp.ncbi.nlm./geo/series/GSE128nnn/GSE128509/suppl/GSE128509_RAW.tar 自己下载，然后解压，在GSE128509_RAW.tar 解压后的文件夹里面运行如下所示代码进行批量改名；

fs=list.files('./','genes.tsv.gz')
fs
samples1=gsub('_genes.tsv.gz','',fs)
samples1
library(stringr)
samples2=str_split(fs,'_',simplify = T)[,2]
samples2 = samples1

lapply(1:length(samples2), function(i){
  x=samples2[i]
  y=samples1[i]
  dir.create(x,recursive = T)
  file.copy(from=paste0(y,'_genes.tsv.gz'),
            to=file.path(x,  'features.tsv.gz' )) 
  file.copy(from=paste0(y,'_matrix.mtx.gz'),
            to= file.path(x, 'matrix.mtx.gz' ) ) 
  file.copy(from=paste0(y,'_barcodes.tsv.gz'),
            to= file.path(x, 'barcodes.tsv.gz' )) 
  
})

得到如下所示的文件夹和路径结构；

outputs/
├── [ 160]  GSM3678221_aMHC-Vegfb_TG
│   ├── [ 11K]  barcodes.tsv.gz
│   ├── [213K]  features.tsv.gz
│   └── [ 20M]  matrix.mtx.gz
├── [ 160]  GSM3678222_aMHC-Vegfb_WT
│   ├── [ 12K]  barcodes.tsv.gz
│   ├── [213K]  features.tsv.gz
│   └── [ 20M]  matrix.mtx.gz
├── [ 160]  GSM3678223_AAV9-Vegfb
│   ├── [ 12K]  barcodes.tsv.gz
│   ├── [213K]  features.tsv.gz
│   └── [ 15M]  matrix.mtx.gz
└── [ 160]  GSM3678224_AAV9-S2
    ├── [ 14K]  barcodes.tsv.gz
    ├── [213K]  features.tsv.gz
    └── [ 13M]  matrix.mtx.gz

4 directories, 12 files

可以看到总共是4个文件夹，12个文件，每个文件夹都是一个单细胞10x样品哦，每个样品都是3个文件，必须文件名字一模一样，路径一模一样的哦.

然后是批量读取和降维聚类分群

批量读取的代码很简单的，我省略了部分质量控制的代码，目前看来是无伤大雅。

library(Seurat)
###### step1:导入数据 ######    
library(data.table)
  dir='GSE128509_RAW/outputs' 
samples=list.files( dir  )
samples 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1]
  folder=file.path( dir ,pro) 
  print(pro)
  print(folder)
  print(list.files(folder))
  sce=CreateSeuratObject(counts = Read10X(folder),
                         project =  pro ,
                         min.cells = 5,
                         min.features = 300)
  
  return(sce)
})
names(sceList)  
 
samples =  gsub('^GSM[1-9]*_','',samples)
samples
names(sceList)  = samples
names(sceList) 

sce.all <- merge(sceList[[1]], y= sceList[ -1 ] ,
                 add.cell.ids =  samples) 

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident)

降维聚类分群，参考前面的例子即可：人人都能学会的单细胞聚类分群注释，这里需要注意的是因为来源于4个样品，所以harmony处理一下：


sce <- NormalizeData(sce, 
                         normalization.method = "LogNormalize",
                         scale.factor = 1e4) 
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))

library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T ) 

sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
                     dims = 1:15) 
sce.all=sce
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
  sce.all=FindClusters(sce.all, #graph.name = "CCA_snn", 
                       resolution = res, algorithm = 1)
}
colnames(sce.all@meta.data)
apply(sce.all@meta.data[,grep("RNA_snn",colnames(sce.all@meta.data))],2,table)
 
p2_tree=clustree(sce.all@meta.data, prefix = "RNA_snn_res.")
ggsave(plot=p2_tree, filename="Tree_diff_resolution.pdf") 
#接下来分析，按照分辨率为0.8进行 
sel.clust = "RNA_snn_res.0.8"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
saveRDS(sce.all, "sce.all_int.rds")

差不多就是在我写这个公众号的同时，代码就跑完了，得到了全部的文件夹和图表。另外，就是降维聚类分群仅仅是算法层面的事情，后续仍然是需要自己根据自己的背景知识，对这些算法上面的单细胞亚群进行生物学命名哦。

因为文章里面提到了，就是纯粹的内皮细胞，所以得到如下所示降维聚类分群：

降维聚类分群

可以看到，仅仅是11群这个并不是内皮，其它都是，而 6,7,12是处于细胞增殖状态的内皮，而 8,10是高表达量 Vwf的内皮，其它都是普通内皮。

如果你具体看AAV9-Vegfb组相对于 AAV9-S2来说，确实是处于增殖状态的内皮细胞这个单细胞亚群比例增多了。

具体可以去读一个解读：Circulation | 冠心病治疗新策略：VEGF-B可促进心内膜源性冠脉血管发育和心肌再生，那里面说的很清楚：

①细胞增殖相关转录本如Ccnb2, Birc5, Top2a等在冠脉ECs和心内膜ECs的表达均增加；
②通常表达于冠脉ECs的Pecam1与发育过程中肌小梁形成所必需的Nrg1在心内膜ECs中的表达均上调；
③心内膜特异性标记基因（Cdh11和Npr3）也表达上调。

AAV9-Vegfb组相对于 AAV9-S2来说

如果有做到原文作者那样的单细胞亚群命名，就需要自己去看文章里面的基因了，反正我使用另外一个整理好的内皮细胞亚群基因，发现非常的不好用：

genes_to_check =list(
  capillaryA=c("Pitp","Plavp","Cd93","Cd74","Tspan7"),
  capillaryB=c("Igfbp7","Car4","Emp2","Hspb1","Kdr"),
  cyc=c("Mki67","Top2a"),
  Vein=c("Vwf","Bgn","Cytl1","Slc6a2","Cpe"),
  Lymphatic=c("Mmrn1","Fxyd6","Fabp4","Ccl21a","Gng11"),
  Artery=c("Plac8","Mgp","Apce","Cxcl12","Tm4sf1"),
  Proliferating=c("Hmgb2","Stmn1","Hmgn2","Tubb5","Tuba1b"),
  SFTP=c("Sftpb","Sftpc","Sftpa1","Scgb1a1","Lyz2")
)

这样的单细胞亚群，太不稳定了。

表达量差异分析

如果仅仅是看单细胞亚群比例变化，其实流式细胞等技术就绰绰有余，而且成本更低可以做多个分组大量样品，单细胞比例也有足够的统计学说服力。相比起来，单细胞转录组因为成本太高，每个分组一般来说也就是两三个样品而已，这样的降维聚类分群后简单的比较不同分组的单细胞亚群比例往往是无法达到统计学显著性。

到目前（2022-01）为止，一个10x单细胞样品费用已经从三四万降价到了2.5万。不过特殊情况下可以更低，比如我们在2021的尾巴在《生信技能树》和《单细胞天地》等公众号推出来的10X单细胞转录组钜惠套餐，详见：2个分组的单细胞项目标准分析，原价15~20万的6个10x单细胞转录组套餐，现价10万。就相当于一个10x单细胞样品仅需1.6万，不过活动价并不是常规价格。

目前，绝大部分小伙伴手上的单细胞转录组数据仍然是一个样品3万人民币左右，而单个项目通常是好几个甚至十几个或者几十个10x样品，取决于财力。

所以我们会看同一个单细胞亚群在不同分组的表达量差异，因为它单细胞转录组虽然每个单细胞本身就成百上千个基因，但是每个单细胞亚群都是有成百上千个细胞，合起来就是两万多个基因基本都是会涉及到，差异分析起来也可以走常规转录组流程啦。

文章里面也是有两个差异分析的火山图，如下所示：