高歌课题组绘制完成 63 种植物功能性转录调控图谱
PlantTFDB – Plant Transcription Factor Database
植物转录因子数据库【planttfdb】的使用
植物比较基因组学和数据库
SCENIC 分析的主要目的是：把单细胞转录组数据结合motif数据库，去构建每个cluster的细胞的regulons，得到每个细胞的regulon activity scores，从而构建转录调控网络，鉴定细胞状态。
- 每种类型的细胞的单细胞转录组数据，进行基因共表达分析，找到可能的转录因子-靶基因的共表达关系从而推动可能发生调控。
- 通过motif分析，把在数据库富集到的结果用于推断转录因子-靶基因的调控关系，即把之前共表达感觉的转录因子-靶基因的关系进一步修剪，找到比较强的关系，对于特定转录因子，总有一系列的下游调控基因，这些基因就统称regulon。
- 挖掘出一系列的因子之后寻找文献支持，并通过实验来验证寻找到的新的因子。
基于共表达和 DNA 模基序 （motif）分析推断基因调控网络。
- step 1：采用 grnboost2 识别并筛选出与 TFs 共表达（co-expression）的基因。注意：其中部分基因仅与 TFs 表达相关，而非靶基因； GENIE3：推断共表达网络
- step 2：使用 RcisTarget 对每个共表达模块进行显著性 Motif 富集，筛选得到显著表达的靶基因，我们将每一对 TFs 与靶基因的组合称为调节子（regulons）。 RcisTarget：TF 结合的 motif 分析
- step 3：采用 AUCell 算法对每一组 regulon 的转录活性进行打分，从而确定 TFs 对靶基因的调控强度（图 2）。 AUCell：鉴定具有激活基因集或者基因网络的细胞
如何获取目标基因的转录因子–Biomart 下载基因和 motif 位置信息
- 获取参考基因组和 GTF 文件 http://plants.ensembl.org/index.html
- Ensembl 中 BioMart 获取 gene.bed 文件 https://mp.weixin.qq.com/s?__biz=MzI5MTcwNjA4NQ==&mid=2247486622&idx=1&sn=55781f9929f478e46c74bc5469d48dbd&scene=21#wechat_redirect
- 植物转录因子数据库获取
下载安装
# R版本SCENIC安装
BiocManager::install(c("GENIE3", "AUCell", "RcisTarget"))
install.packages('zoo')
BiocManager::install(c("mixtools", "rbokeh"))
BiocManager::install(c("NMF", "pheatmap", "Rtsne", "R2HTML"))
BiocManager::install(c("doMC", "doRNG"))
devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)
BiocManager::install(c("SingleCellExperiment"))
devtools::install_github("aertslab/SCENIC")
packageVersion("SCENIC")
# Python版本SCENIC安装
pip install pyscenic
输入
SCENIC 需要输入的是单细胞 RNA-seq 表达矩阵—— 每列对应于样品（细胞），每行对应一个基因。基因 ID 应该是 gene-symbol 并存储为 rownames （尤其是基因名字部分是为了与 RcisTarget 数据库兼容）；表达数据是 Gene 的 reads count。提供原始的或 Normalized UMI count，无论是否 log 转换，或使用 TPM 值，结果相差不大。
# 输入.loom文件
# .loom文件可以通过SCopeLoomR包直接导入SCENIC
library(SCopeLoomR)
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
cellInfo <- get_cellAnnotation(loom)
close_loom(loom)
dim(exprMat)
# 10x输出文件夹输入
# 通过Seurat进行输入
singleCellMatrix <- Seurat::Read10X(data.dir="data/pbmc3k/filtered_gene_bc_matrices/hg19/")
cellInfo <- data.frame(seuratCluster=Idents(seuratObject))
# 输入R对象
sce <- load_as_sce(dataPath) # any SingleCellExperiment object
exprMat <- counts(sce)
cellInfo <- colData(sce)
# 从GEO数据库
# 获取数据及数据标注
# dir.create("SCENIC_MouseBrain"); setwd("SCENIC_MouseBrain") # if needed
# 安装GEOquery包
if (!requireNamespace("GEOquery", quietly = TRUE)) BiocManager::install("GEOquery")
library(GEOquery)
# 数据下载
geoFile <- getGEOSuppFiles("GSE60361", makeDirectory=FALSE)
gzFile <- grep("Expression", basename(rownames(geoFile)), value=TRUE)
txtFile <- gsub(".gz", "", gzFile)
gunzip(gzFile, destname=txtFile, remove=TRUE)
# 数据读取
library(data.table)
geoData <- fread(txtFile, sep="\t")
geneNames <- unname(unlist(geoData[,1, with=FALSE]))
exprMatrix <- as.matrix(geoData[,-1, with=FALSE])
rm(geoData)
dim(exprMatrix)
rownames(exprMatrix) <- geneNames
exprMatrix <- exprMatrix[unique(rownames(exprMatrix)),]
exprMatrix[1:5,1:4]
# Remove file downloaded:
file.remove(txtFile)
SCENIC 分析流程
建立基因调控网络（Gene Regulation Network，GRN）：
- 基于共表达识别每个转录因子 TF 的潜在靶标。
过滤表达矩阵并运行 GENIE3 或者 GRNBoost，它们是利用表达矩阵推断基因调控网络的一种算法，能得到转录因子和潜在靶标的相关性网络；
将目标从 GENIE3 或者 GRNBoost 格式转为共表达模块。 - 根据 DNA 模序分析（motif）选择潜在的直接结合靶标（调节因子）（利用 RcisTarget 包：TF 基序分析）
确定细胞状态及其调节因子：
- 分析每个细胞中的网络活性（AUCell） 在细胞中评分调节子（计算 AUC）
详细分析流程
# 创建新目录
## 中间文件和图文件将会保存到int文件夹中
## 重要的输出结果会保存到output文件夹中
dir.create("SCENIC_ara")
setwd("SCENIC_ara")
# 数据输入
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
cellInfo <- get_cellAnnotation(loom)
close_loom(loom)
dim(exprMat)
## 细胞信息或者表行数据
head(cellInfo)
cellInfo <- data.frame(cellInfo)
cellTypeColumn <- "Class"
colnames(cellInfo)[which(colnames(cellInfo)==cellTypeColumn)] <- "CellType"
cbind(table(cellInfo$CellType))
dir.create("int")
saveRDS(cellInfo, file="int/cellInfo.Rds")
colVars <- list(CellType=c("microglia"="forestgreen", # 设置颜色
"endothelial-mural"="darkorange",
"astrocytes_ependymal"="magenta4",
"oligodendrocytes"="hotpink",
"interneurons"="red3",
"pyramidal CA1"="skyblue",
"pyramidal SS"="darkblue"))
colVars$CellType <- colVars$CellType[intersect(names(colVars$CellType), cellInfo$CellType)]
saveRDS(colVars, file="int/colVars.Rds")
plot.new(); legend(0,1, fill=colVars$CellType, legend=names(colVars$CellType))
# 初始化SCEMIC设置
library(SCENIC)
org="mgi" # or hgnc, or dmel
dbDir="cisTarget_databases" # RcisTarget databases location
myDatasetTitle="SCENIC example on Mouse brain" # choose a name for your analysis
data(defaultDbNames)
dbs <- defaultDbNames[[org]]
scenicOptions <- initializeScenic(org=org, dbDir=dbDir, dbs=dbs, datasetTitle=myDatasetTitle, nCores=10)
## Motif databases selected:
## mm9-500bp-upstream-7species.mc9nr.feather
## mm9-tss-centered-10kb-7species.mc9nr.feather
# Modify if needed
scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
scenicOptions@inputDatasetInfo$colVars <- "int/colVars.Rds"
# Databases:
# scenicOptions@settings$dbs <- c("mm9-5kb-mc8nr"="mm9-tss-centered-5kb-10species.mc8nr.feather")
# scenicOptions@settings$db_mcVersion <- "v8"
# Save to use at a later time...
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
# 共表达网络建立
# 根据表达数据推测潜在的转录因子靶点
# 输入文件：经过筛选的表达矩阵 + 转录因子列表
# 输出文件：相关性矩阵
# 存在两种方法：GENIE3 或 GRNBoost(处理大数据时更优)
# 当数据集存在大量低质量细胞时，会消耗大量时间，可以挑选其中一部分的细胞来推断调控网络，在之后的第三步中可以使用全部的细胞在AUCell中重新评估。
## 基因筛选和选择
### 基于基因表达总数和其表达的细胞数目的双筛选标准进行筛选基因
# (Adjust minimum values according to your dataset)
genesKept <- geneFiltering(exprMat, scenicOptions=scenicOptions,
minCountsPerGene=3*.01*ncol(exprMat),
minSamples=ncol(exprMat)*.01)
exprMat_filtered <- exprMat[genesKept, ] # 表达矩阵筛选
dim(exprMat_filtered)
rm(exprMat)
## 相关性分析（spearman方法)
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered, scenicOptions)
# 建立和计分GRN
library(SCENIC)
scenicOptions <- readRDS("int/scenicOptions.Rds")
scenicOptions@settings$verbose <- TRUE
scenicOptions@settings$nCores <- 10
scenicOptions@settings$seed <- 123
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"]
runSCENIC_1_coexNetwork2modules(scenicOptions)
runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) #** Only for toy run!!
runSCENIC_3_scoreCells(scenicOptions, logMat)
# 对网络活性的结果进行二维化（on/off） 可选
aucellApp <- plotTsne_AUCellApp(scenicOptions, logMat)
savedSelections <- shiny::runApp(aucellApp)
# Save the modified thresholds:
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")
# scenicOptions@settings$devType="png"
runSCENIC_4_aucell_binarize(scenicOptions)
# 在调节子活性的矩阵上对细胞降维并聚类
nPcs <- c(5) # For toy dataset
# nPcs <- c(5,15,50)
scenicOptions@settings$seed <- 123 # same seed for all of them
# Run t-SNE with different settings:
fileNames <- tsneAUC(scenicOptions, aucType="AUC", nPcs=nPcs, perpl=c(5,15,50))
fileNames <- tsneAUC(scenicOptions, aucType="AUC", nPcs=nPcs, perpl=c(5,15,50), onlyHighConf=TRUE, filePrefix="int/tSNE_oHC")
# Plot as pdf (individual files in int/):
fileNames <- paste0("int/",grep(".Rds", grep("tSNE_", list.files("int"), value=T), value=T))
scenicOptions@settings$defaultTsne$aucType <- "AUC"
scenicOptions@settings$defaultTsne$dims <- 5
scenicOptions@settings$defaultTsne$perpl <- 15
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
主要运行流程
library(foreach)#注意一定要加载这个包，不然会报错！
library(SCENIC)
### Load data
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
cellInfo <- get_cellAnnotation(loom)
close_loom(loom)
### Initialize settings
library(SCENIC)
scenicOptions <- initializeScenic(org="mgi", dbDir="cisTarget_databases", nCores=10)
# scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
### Co-expression network
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions)
### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
runSCENIC_1_coexNetwork2modules(scenicOptions)
runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
# Export:
export2scope(scenicOptions, exprMat)
# Binarize activity?
# aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_log)
# 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)
### Exploring output
# Check files in folder 'output'
# .loom file @ http://scope.aertslab.org
# output/Step2_MotifEnrichment_preview.html in detail/subset:
motifEnrichment_selfMotifs_wGenes <- loadInt(scenicOptions, "motifEnrichment_selfMotifs_wGenes")
tableSubset <- motifEnrichment_selfMotifs_wGenes[highlightedTFs=="Sox8"]
viewMotifs(tableSubset)
# output/Step2_regulonTargetsInfo.tsv in detail:
regulonTargetsInfo <- loadInt(scenicOptions, "regulonTargetsInfo")
tableSubset <- regulonTargetsInfo[TF=="Stat6" & highConfAnnot==TRUE]
viewMotifs(tableSubset)
SCENIC详细流程分解学习–便于后续植物研究
第一节：GENIE3 共表达网络分析
GENIE3是一种基于提供的表达矩阵来预测基因调控网络的R包，GENIE3通过确定最能解释每个靶基因表达的转录因子表达模式，整合转录因子信息来构建调控网络。重点是需要TF信息
- GENIE3的输入：表达矩阵，可以是count或TPM等，但不应该标准化或者归一化处理。表达矩阵必须为matrix
- 需要考虑基因的个数
相关性分析简单来说就是，如果geneA表达模式变化，geneB同样跟着变化(相同或者相反)那么，geneA和geneB就存在较强相关性。
第二节：cigTarget motif富集分析
常见的motif富集分析软件：HOMER、MEME
RcisTarget是一个用于识别基因列表中转录因子结合基序的R包；
- Gene-motif rankings：提供每个motif的所有基因排名
- 转录因子的motif注释
cisTarget()运行步骤：1. motif富集分析；2. motif-TF注释；3. 选择重要的基因
# HOMER中motif分析原理理解
## 预处理
## 1. 序列提取
若给出的基因组位置信息，则提取出来的是对应的基因组序列；
若给出的是基因accession number，则需要选择适当的promoter区域
## 2. 背景选择
如果使用的是基因组位置，将从整个基因组序列抓取GC含量一致的序列作为背景序列。若进行的是promoter分析，则将所有的promoter作为背景。
## 3. GC含量矫正
将目标序列和背景序列对GC含量进行分bin（5%区间），背景序列通过调节权重得到和目标序列同样的GC含量分布。这可以使得HOMER在分析来自CpG岛时的序列时，不会简单的找到那些GC富集的motif。
## 4. 自动标准化
考虑其他可能的影响因素
## 5. 解析输入序列
根据设置的motif长度，将输入序列解析为寡核苷酸，并生成一个寡核苷酸频度表。该表只记录出现的寡核苷酸和其出现的次数，可以增加motif搜索时的效率，但是也会破坏寡核苷酸与其原始序列的关系。
## 6. 寡核苷酸自动均一化
将那些较短的寡核苷酸与较长的寡核苷酸相平衡。
## 7. 全局搜索
若某个motif富集，则其存在的寡核苷酸也同样富集。
## 8. 矩阵优化
## 9. Mask and Repeat
如何获取目标基因的转录因子（上）——Biomart下载基因和motif位置信息
如何获取目标基因的转录因子（下）——Linux命令获取目标基因TF
