ixxmu
commented
2 years ago bulk转录组当做单细胞转录组去做SCENIC分析 by 东林的扯淡小屋
YDL<- readRDS("mm_ternimalE.rds")YDLset.seed(123)library(Seurat)library(tidyverse)library(patchwork)library(SCENIC)##==分析准备==##dir.create("SCENIC")dir.create("SCENIC/int")scRNA <- YDLsetwd("./SCENIC")##准备细胞meta信息##准备细胞meta信息cellInfo <- data.frame(scRNA@meta.data)colnames(cellInfo)[which(colnames(cellInfo)=="orig.ident")] <- "sample"#colnames(cellInfo)[which(colnames(cellInfo)=="seurat_clusters")] <- "cluster"colnames(cellInfo)[which(colnames(cellInfo)=="celltype")] <- "celltype"cellInfo <- cellInfo[,c("sample","celltype")]saveRDS(cellInfo, file="int/cellInfo.Rds")##准备表达矩阵#为了节省计算资源，随机抽取1000个细胞的数据子集#subcell <- sample(colnames(scRNA),1000)#scRNAsub <- scRNA[,subcell]#saveRDS(scRNAsub, "scRNAsub.rds")exprMat <- as.matrix(scRNA@assays$RNA@counts)dim(exprMat)# [1] 32285 1000#head(cellInfo)##设置分析环境mydbDIR <- "/data/yudonglin/scenic/cisTarget_databases/cisTarget_databases/"mydbs <- c("mm9-500bp-upstream-7species.mc9nr.feather","mm9-tss-centered-10kb-7species.mc9nr.feather")names(mydbs) <- c("500bp", "10kb")scenicOptions <- initializeScenic(org="mgi",nCores=20,dbDir=mydbDIR,dbs = mydbs,datasetTitle = "SCENIC")saveRDS(scenicOptions, "int/scenicOptions.rds")genesKept <- geneFiltering(exprMat, scenicOptions,minCountsPerGene = 3 * 0.01 * ncol(exprMat),minSamples = ncol(exprMat) * 0.01)exprMat_filtered <- exprMat[genesKept, ]##计算相关性矩阵runCorrelation(exprMat_filtered, scenicOptions)##TF-Targets相关性回归分析exprMat_filtered_log <- log2(exprMat_filtered+1)runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)#这一步消耗的计算资源非常大，个人电脑需要几个小时的运行时间##推断共表达模块runSCENIC_1_coexNetwork2modules(scenicOptions)scenicOptions <- initializeScenic(org="mgi",nCores=9,dbDir=mydbDIR,dbs = mydbs,datasetTitle = "erythrogenesis")##推断转录调控网络（regulon）runSCENIC_2_createRegulons(scenicOptions)#以上代码可增加参数coexMethod=c("w001", "w005", "top50", "top5perTarget", "top10perTarget", "top50perTarget"))#默认6种方法的共表达网络都计算，可以少选几种方法以减少计算量scenicOptions <- initializeScenic(org="mgi",nCores=1,dbDir=mydbDIR,dbs = mydbs,datasetTitle = "erythrogenesis")##==regulon活性评分与可视化==####regulons计算AUC值并进行下游分析exprMat_all <- as.matrix(scRNA@assays$RNA@counts)exprMat_all <- log2(exprMat_all+1)saveRDS(exprMat_all,"exprMat_all.rds")# rm(list=ls())# scenicOptions<-readRDS("scenicOptions.rds")# exprMat_all<-readRDS("exprMat_all.rds")runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)# #使用shiny互动调整阈值# aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_all)# savedSelections <- shiny::runApp(aucellApp)# #保存调整后的阈值# newThresholds <- savedSelections$thresholds# scenicOptions@fileNames$int["aucell_thresholds",1] <- "int/newThresholds.Rds"# saveRDS(newThresholds, file=getIntName(scenicOptions, "aucell_thresholds"))# saveRDS(scenicOptions, file="int/scenicOptions.Rds")runSCENIC_4_aucell_binarize(scenicOptions, exprMat=exprMat_all)savehistory("scenic_code.txt")
只有20个样本，视为20个单细胞：
也算是一个骚操作了hhh，不过要求每个样本至少三个重复才行。
https://mp.weixin.qq.com/s/WvbszojyMMIoMus0smMy7g