写在前面的话
拟时序分析已经是单细胞和空间转录组中的应用文章中的必需分析，干货记录分享，简单易懂！
1读取数据提取亚群
library(Seurat)
scRNA<- readRDS("scRNA.rds")
DimPlot(scRNA, reduction = "umap",label = TRUE)
Mono <- subset(mono, idents = c("1"))
Mono <- SCTransform(Mono,method="glmGamPoi",vars.to.regress=c("percent.MT","G2M.Score","S.Score")) %>% RunPCA()
library(harmony)
Mono <- RunHarmony(Mono, group.by.vars = "orig.ident",
assay.use = "SCT", max.iter.harmony = 10)
Mono <- RunUMAP(Mono,reduction = "harmony", dims = 1:30, verbose = FALSE)
Mono <- FindNeighbors(Mono, reduction = "harmony",
dims = 1:30)
2细胞类型注释
library(SingleR)
library(celldex)
ref <- celldex::HumanPrimaryCellAtlasData()
Mono.count=Mono[["RNA"]]@counts
common_gene <- intersect(rownames(Mono.count), rownames(ref))
ref <- ref[common_gene,]
Mono.count_forref <- Mono.count[common_gene,]
Mono.count_forref.se <- SummarizedExperiment(assays=list(counts=Mono.count_forref))
Mono.count_forref.se <- logNormCounts(Mono.count_forref.se)
pred.main.ref <- SingleR(test = Mono.count_forref.se, ref = ref, labels = ref$label.main)
result_main_ref <- as.data.frame(pred.main.ref$labels)
result_main_ref$CB <- rownames(pred.main.ref)
colnames(result_main_ref) <- c('ref_Main', 'CB')
DimPlot(Mono, reduction = "umap", group.by = "ref_Main",label = TRUE,repel = TRUE,ncol=1)+
DimPlot(Mono, reduction = "umap",label = TRUE,repel = TRUE)&
theme(axis.line = element_blank(), axis.ticks= element_blank(),axis.text = element_blank())
## 注释回第一次降维聚类展示
metadata <- scRNA@meta.data
metadata$cell_type <- NA
metadata$cell_type <- Mono@meta.data[match(rownames(scRNA@meta.data), rownames(Mono@meta.data)), 'ref_Main']
celltype_names <- NULL
for(i in 1:dim(metadata)[1]){
# print(i)
sub_data <- metadata[i,]
if(is.na(sub_data$cell_type)){
# print('Change value')
sub_data$cell_type <- sub_data$RNA_snn_res.0.3
celltype_names <- c(celltype_names, sub_data$cell_type)
}else{
celltype_names <- c(celltype_names, sub_data$cell_type)
}
}
metadata$cell_type <- celltype_names
metadata$cell_type <- factor(metadata$cell_type)
scRNA@meta.data <- metadata
DimPlot(scRNA,reduction = "umap",group.by = "cell_type",label = TRUE,repel = TRUE)&
theme(axis.line = element_blank(), axis.ticks= element_blank(),axis.text = element_blank())
3slingshot推断细胞分化轨迹
BiocManager::install("slingshot")
library(slingshot)
BiocManager::install("tradeSeq")
library(tradeSeq)
library(BiocParallel)
#准备分析文件
sce <- as.SingleCellExperiment(Mono, assay = "RNA") # SCT
#run slingshot
sce_slingshot1 <- slingshot(sce, #这里是SingleCellExperiment单细胞对象
reducedDim = 'UMAP', #降维方式
clusterLabels = sce$ref_Main.y) #celltype
#start.clus = ''#轨迹起点
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sce_slingshot1$slingPseudotime_1, breaks=100)]
plot(reducedDims(sce_slingshot1)$UMAP, col = plotcol, pch=16, asp = 1, cex = 0.8)
lines(SlingshotDataSet(sce_slingshot1), lwd=2, col='black')
4CytoTRACE推断细胞干性
下载源文件https://cytotrace.stanford.edu/CytoTRACE_0.3.3.tar.gz
devtools::install_local("CytoTRACE_0.3.3.tar.gz")
library(CytoTRACE)
#检查是否有缺失值，
table(is.na(Mono$ref_Main.y))
#将其转换为字符型，并将其命名为 phenotype这个变量，用于后续的分析。
phenotype = as.character(Mono$ref_Main.y)
names(phenotype) <- rownames(Mono@meta.data)
#将 Mono 对象中的 RNA 计数矩阵提取出来，并转换为一个普通的矩阵
mat <- as.matrix(Mono@assays$RNA@counts)
# 接着，用CytoTRACE 函数分化潜能评估
results <- CytoTRACE(mat = mat,ncores = 1) #可以多核运行
# 使用umap的坐标映射干性得分
emb <- Mono@reductions[["umap"]]@cell.embeddings
# plotCytoTRACE(results,
phenotype = phenotype,#细胞类型注释
emb = emb,
outputDir = "./")
5基因拟时变化表达
slingsce<-SlingshotDataSet(sce_slingshot1)
pseudotimeED <- slingPseudotime(slingsce, na = FALSE)
cellWeightsED <- slingCurveWeights(slingsce)
counts<-sce_slingshot1@assays@data@listData$counts
# 拟合负二项式模型确定nknots
icMat <- evaluateK(counts = counts,
sds = slingsce,
nGenes = 50,
k = 3:10,
verbose = T,
plot = TRUE)
sce_slinghot <- fitGAM(counts = counts, pseudotime = pseudotimeED, cellWeights = cellWeightsED, nknots = 5, verbose = T
,BPPARAM = MulticoreParam(20),parallel=T)
# 与伪时间存在显著关联的基因
rowData(sce_slinghot)$assocRes <- associationTest(sce_slinghot, lineages = TRUE, l2fc = log2(2))
assocRes <- rowData(sce_slinghot)$assocRes
oAsso <- order(assocRes$waldStat, decreasing = TRUE)
Geneasso <- names(sce_slinghot)[oStart[1]]
plotSmoothers(sce_slinghot, assays(sce_slinghot)$counts,
gene =Geneasso, alpha = 0.6, border = T, lwd = 2)+
ggtitle(Geneasso)
# 谱系内比较
startRes <- startVsEndTest(sce_slinghot,lineages=T)
oStart <- order(startRes$waldStat, decreasing = TRUE)
sigGeneStart <- names(sce_slinghot)[oStart[1]]
plotSmoothers(sce_slinghot, assays(sce_slinghot)$counts, gene = sigGeneStart)+
ggtitle(sigGeneStart)
plotGeneCount(slingsce, assays(sce_slinghot)$counts, gene = sigGeneStart)+
ggtitle(sigGeneStart)
# 谱系间比较
endRes <- diffEndTest(sce_slinghot,pairwise=T)
oEnd <- order(endRes$waldStat, decreasing = TRUE)
sigGene <- names(sce_slinghot)[o[1]]
plotSmoothers(sce_slinghot, assays(sce_slinghot)$counts, gene = sigGene)+
ggtitle(sigGene)
plotGeneCount(slingsce, assays(sce_slinghot)$counts, gene = sigGene)+
ggtitle(sigGene)
6通路富集评分拟时变化表达
# 把表达矩阵换成通路富集评分的矩阵，以AUcell评分为例
library(AUCell)
library(clusterProfiler)
cells_rankings <- AUCell_buildRankings(Mono@assays$RNA@data)
Hallmarker <- read.gmt("h.all.v2023.1.Hs.symbols.gmt")
geneSets <- lapply(unique(Hallmarker$term),
function(x){print(x);Hallmarker$gene[Hallmarker$term == x]})
names(geneSets) <- unique(Hallmarker$term)
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings,
aucMaxRank=nrow(cells_rankings)*0.1)
sc_AUC <- getAUC(cells_AUC)
sce_slinghot <- fitGAM(counts = sc_AUC, pseudotime = pseudotimeED, cellWeights = cellWeightsED, nknots = 5, verbose = T
,BPPARAM = MulticoreParam(20),parallel=T)
plotSmoothers(sce_slinghot, assays(sce_slinghot)$counts,
gene ="HALLMARK_HYPOXIA", alpha = 0.6, border = T, lwd = 2)
