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

生信技能树

发布于 2022-06-08 13:40:26

2.7K00

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

单细胞的多组的设计（比如正常组与给药组）可以为细胞类型水平比较提供以往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.ahajournals.org/doi/10.1161/CIRCULATIONAHA.120.050635

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

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

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

降维聚类分群，参考前面的例子即可：人人都能学会的单细胞聚类分群注释 ，它数据在：https://www.ncbi.nlm.nih.gov/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.nih.gov/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")
)