https://www.jianshu.com/p/5ae987ba68e0

    1. setwd("~/project")
    2. load('singleBB.Rdata')
    3. experiment.aggregate@meta.data$celltype4 = "NA"
    4. # 更改 celltype 信息，设置细胞群新名称
    5. experiment.aggregate@meta.data[which(experiment.aggregate@meta.data$celltype3 %in% c('GC B cells','Memory B cells','Naive B cells')), "celltype4"] = "CD20+ B cells"
    6. experiment.aggregate@meta.data[which(experiment.aggregate@meta.data$celltype3 %in% c('IGA+ B cells','IGG+ B cells')), "celltype4"] = "Plasma B cells"
    9. cellInfo <- data.frame(experiment.aggregate@meta.data)
    10. view(experiment.aggregate@meta.data)
    11. colnames(cellInfo)[which(colnames(cellInfo)=="seurat_clusters")] <- "cluster"
    12. colnames(cellInfo)[which(colnames(cellInfo)=="celltype3")] <- "celltype3"
    13. colnames(cellInfo)[which(colnames(cellInfo)=="celltype4")] <- "celltype4"
    14. cellInfo <- cellInfo[,c("cluster","celltype3","celltype4")]
    15. saveRDS(cellInfo, file="int/cellInfo.Rds")
    16. ##准备表达矩阵
    17. #为了节省计算资源，随机抽取1000个细胞的数据子集
    18. subcell <- sample(colnames(experiment.aggregate),1000)
    19. exprMat1 <- experiment.aggregate[,subcell]
    20. DimPlot(exprMat1, reduction = "tsne", label = T,group.by = 'celltype3')
    21. table(exprMat@meta.data$celltype3)
    24. library(SCENIC)
    25. rm(list=ls())
    26. dir.create("SCENIC")
    27. dir.create("SCENIC/int")
    28. setwd("~/project/SCENIC")
    30. ##准备scenic输入文件
    32. exprMat <- GetAssayData(exprMat1, assay = 'RNA', slot = 'data') %>% as.matrix()
    35. mydbDIR <- "../Resource/hg38_scenic/"
    36. dir(mydbDIR)
    37. mydbs <- c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather",
    38.            "hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather")
    39. names(mydbs) <- c("500bp", "10kb")
    40. #小鼠org="mgi"
    41. scenicOptions <- initializeScenic(org = "hgnc", 
    42.                                   nCores = 16,
    43.                                   dbDir = mydbDIR, 
    44.                                   dbs = mydbs,
    45.                                   datasetTitle = "HNC")
    46. saveRDS(scenicOptions, "int/scenicOptions.rds")
    47. #scenicOptions = readRDS("int/scenicOptions.rds")
    48. #scenicOptions@settings$nCores <- 10
    50. ##如果需要高性能计算服务，请保存以下文件联系Kinesin
    53. ##基因过滤
    54. genesKept <- geneFiltering(exprMat, scenicOptions, 
    55.                            minCountsPerGene = 0.015 * ncol(exprMat), 
    56.                            minSamples = ncol(exprMat) * 0.01)
    57. exprMat_filtered <- exprMat[genesKept, ]
    59. ##计算相关性矩阵
    60. runCorrelation(exprMat_filtered, scenicOptions)
    62. ##计算TF-targets相关性
    63. runGenie3(exprMat_filtered, scenicOptions, nParts = 20)
    64. runSCENIC_1_coexNetwork2modules(scenicOptions)
    66. ##推断转录调控网络（regulon）
    67. #此步运行时间长且极耗内存，慎用多线程！！！
    68. scenicOptions@settings$nCores <- 2
    69. runSCENIC_2_createRegulons(scenicOptions)
    71. exprMat_all <- as.matrix(experiment.aggregate@assays$RNA@data)
    73. scenicOptions@settings$nCores <- 6
    74. runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)
    75. runSCENIC_4_aucell_binarize(scenicOptions, exprMat=exprMat_all)
    1. cellInfo <- readRDS("int/cellInfo.Rds")
    2. celltype = subset(cellInfo,select = 'celltype4')
    3. celltype <- celltype %>%
    4.   arrange(celltype$celltype4)
    5. dd <- rownames(celltype )
    7. Regulon <-readRDS("int/3.4_regulonAUC.Rds")
    8. Regulon <- Regulon@assays@data@listData$AUC
    9. Regulon_all <- Regulon[,dd]
    11. library(pheatmap)
    13. pheatmap(Regulon_all, show_colnames=F, annotation_col=celltype,
    14.          filename = 'scenic_seurat/myAUCmatrix_heatmap.png',
    15.          width = 6, height = 20)
    16. #整体看一下，然后挑选你觉得有意义且差异大的转录因子继续做热图
    17. my.regulons <- c('SPIB (795g)','JUNB_extended (1298g)','NR3C1 (1424g)','BCLAF1 (1908g)','ELF1 (2242g)',
    18.                  'CREM (1318g)','TAF7_extended (2245g)','BCL11A_extended (874g)','CHD2_extended (1353g)',
    19.                  'JUND (413g)','HMGB1 (18g)','REL (317g)',
    20.                  'IRF8 (103g)','CREB3L2 (1268g)','XBP1 (1737g)','PRDM1_extended (28g)',
    21.                  'KLF13_extended (353g)')
    22. myAUCmatrix <- Regulon_all[rownames(Regulon_all)%in%my.regulons,dd]
    23. result3<- t(scale(t(myAUCmatrix)))
    24. result3[result3>1.5]=1.5
    25. result3[result3< -1.5]= -1.5
    26. E1 <- pheatmap(result3, show_colnames=F,annotation_col=celltype,
    27.          cluster_rows = T,#行聚类
    28.          cluster_cols = F,
    29.          #scale = "row",
    30.          color =colorRampPalette(c('blue','white', "red"))(100), cellwidth = 0.018, cellheight = 15,# 格子比例
    31.          fontsize = 10)
    33. E1
    34. pdf("E1.pdf",width = 7,height = 4)
    35. print(E1)
    36. dev.off()
    1. AUCmatrix <- readRDS("int/3.4_regulonAUC.Rds")
    2. AUCmatrix <- AUCmatrix@assays@data@listData$AUC
    3. AUCmatrix <- data.frame(t(AUCmatrix), check.names=F)
    4. identical(rownames(AUCmatrix),rownames(experiment.aggregate@meta.data))
    5. RegulonName_AUC <- colnames(AUCmatrix)
    6. RegulonName_AUC <- gsub(' \\(','_',RegulonName_AUC)
    7. RegulonName_AUC <- gsub('\\)','',RegulonName_AUC)
    8. colnames(AUCmatrix) <- RegulonName_AUC
    9. scRNAauc <- AddMetaData(experiment.aggregate, AUCmatrix)
    10. E3 = FeaturePlot(scRNAauc, features='REL_317g', label=T, reduction = 'tsne',cols =   c("lightgrey", "#CC0033"),)
    11. E3
    12. pdf("E3.pdf",width = 5,height = 4)
    13. print(E3)
    14. dev.off()
    BINmatrix <- readRDS("int/4.1_binaryRegulonActivity.Rds")
BINmatrix <- data.frame(t(BINmatrix), check.names=F)
RegulonName_BIN <- colnames(BINmatrix)
RegulonName_BIN <- gsub(' \\(','_',RegulonName_BIN)
RegulonName_BIN <- gsub('\\)','',RegulonName_BIN)
colnames(BINmatrix) <- RegulonName_BIN
scRNAbin <- AddMetaData(experiment.aggregate, BINmatrix)
E10 = FeaturePlot(scRNAbin, features='CREM_1318g', label=T, reduction = 'tsne',cols =   c("lightgrey", "#CC0033"))
E10
pdf("E10.pdf",width = 5,height = 4)
print(E10)
dev.off()