今天开始测试SCENIC的安装，先测试一下R版本的。
先安装一些必要的包。
## Required
BiocManager::install(c("AUCell", "RcisTarget"))
BiocManager::install(c("GENIE3")) # Optional. Can be replaced by GRNBoost
## Optional (but highly recommended):
BiocManager::install(c("zoo", "mixtools", "rbokeh"))
# For various visualizations and perform t-SNEs:
BiocManager::install(c("DT", "NMF", "ComplexHeatmap", "R2HTML", "Rtsne"))
# To support paralell execution (not available in Windows):
BiocManager::install(c("doMC", "doRNG"))
# To export/visualize in http://scope.aertslab.org
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)
#最后安装SCENIC
devtools::install_github("aertslab/SCENIC")
=========下载库文件======
官网给出的下载地址为：
#人类
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather")
#小鼠
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-7species.mc9nr.feather")
#果蝇
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/drosophila_melanogaster/dm6/flybase_r6.02/mc8nr/gene_based/dm6-5kb-upstream-full-tx-11species.mc8nr.feather")
但是，我去根据链接找的时候，在最新版本里面是根本就找不到的，需要去old文件夹中去找，才能下载到需要的文件。
#人类
https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/
#小鼠
https://resources.aertslab.org/cistarget/databases/old/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/
#果蝇
https://resources.aertslab.org/cistarget/databases/old/drosophila_melanogaster/dm6/flybase_r6.02/mc8nr/gene_based/
#官网建议放在一个特定的目录中。
dir.create("cisTarget_databases");
setwd("cisTarget_databases")
for(featherURL in dbFiles)
{
download.file(featherURL, destfile=basename(featherURL))
}
===========例子测试===========
我们先用包里自带的例子做测试，自带的是一个loom文件。
#先加载数据
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
library(SCENIC)
loom <- open_loom(loomPath)
exprMat <- get_dgem(loom)
cellInfo <- get_cell_annotation(loom)
close_loom(loom)
### 初始化设置，相当于导入评分数据库
scenicOptions <- initializeScenic(org="mgi", dbDir="cisTarget_databases", nCores=10)
#但是运行这一步，我这里报错了，说是找不到motifAnnotations_mgi这个对象。
网上查了一下，也有别人碰到了同样的错误。网上建议的解决办法是：
运行上面的那个命令之前，先运行下面的这2个命令。试了试，确实不报错了。
data(list="motifAnnotations_mgi_v9", package="RcisTarget")
motifAnnotations_mgi <- motifAnnotations_mgi_v9
# 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)
#有个警告，不确定是不是库里面用的ID，和这个数据集的ID不是很匹配。
共表达网络产生了一个结果很多都在int文件夹，比如1.2_corrMat.Rds，这里面存储的就是基因之间相关性矩阵。
1.3_GENIE3_weightMatrix_part_1.Rds GENIE3中间结果，矩阵格式，我这里例子产生了10个；
1.4_GENIE3_linkList.Rds GENIE3的最终结果，是把1.3那一堆文件合并在一起了。文件有三列：TF为转录因子名称，Target为潜在靶基因的名称，weight时TF与Target之间的相关性权重。
1.5_weightPlot.pdf 是关于weight的，类似于密度分布图吧。
### Build and score the GRN，主要推断共表达模块
Step 1:
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"]
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
GENIE3计算了转录因子与每个基因的相关性，接下来需要过滤掉低相关性的组合，形成共表达基因集。作者一共推荐6种都用，这六种分别为：
· w001：以每个TF为核心保留weight>0.001的基因形成共表达模块；
· w005：以每个TF为核心保留weight>0.005的基因形成共表达模块；
· top50：以每个TF为核心保留weight值top50的基因形成共表达模块；
· top5perTarget：每个基因保留weight值top5的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
· top10perTarget：每个基因保留weight值top10的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
· top50perTarget：每个基因保留weight值top50的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
结果文件为1.6_tfModules_asDF.Rds，一共四列：method为6种过滤方法；corr为runCorrelation(exprMat_filtered, scenicOptions)得到的结果，1代表激活，-1代表抑制，0代表中性。SCENIC只会采用值为1，即激活的数据用于后续分析。
1.6_tfModules_asDF.Rds 具体结果如下图所示：
Step 2:
#上一步找到了每个转录因子强相关的靶基因，这一步要修建共表达模块形成有生物学意义的调控单元（regulons，我自己倾向于叫module）。对每个共表达模块进行motif富集分析，保留显著富集的motif；这一步依赖gene-motif评分数据库，行为基因，列为motif，值为排名，也就是我们下载的cisTarget数据库。使用数据库对motif进行TF注释，得到高可信度、低可信度。数据库直接注释和同源基因推断的TF为高可信度，使用motif序列相似性注释的TF为低可信度结果。
用保留的motif对共表达模块里的基因进行打分（依据cisTarget数据库），识别显著高分的基因（理解为motif里这些基因的TSS很近）。删除共表达模块中与motif评分不高的基因，剩下的基因集被称为调控单元（regulon）。
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget"))
输出结果在：output/Step2_MotifEnrichment.tsv，共有11列：motifDb geneSet motif NES AUC highlightedTFs TFinDB TF_highConf TF_lowConf nEnrGenes rankAtMax enrichedGenes。
output/Step2_regulonTargetsInfo.tsv 这里存储最终的regulon文件，将第一个文件的数据整合在了一起。
output/Step2_MotifEnrichment_preview.html 还有个html的页面可以浏览。但是不知道为啥，我的流浪器不显示motif的logo图片。
Step 3: Regulon活动评分与可视化
一个regulon=一个TF与其Target gene的基因集，regulon的名称有两种写法：
TF名称+extended+靶基因数目：转录因子与所有靶基因组成的基因调控网络
TF名称+靶基因数目：转录因子与高可信靶基因（即highConfAnnot=TRUE的基因）组成的基因调控网络。AUCell对每个regulon在各个细胞中的活性进行评分。评分的基础是基因表达值，分数越高代表基因集的激活程度越高。这一步要用所有细胞做计算。
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
Step 4：生成二进制的regulon活性矩阵
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
===============几种结果可视化形式=============
1. Regulon Activity heatmap
step3有所有regulon在细胞的AUCscore热图Step3_RegulonActivity_heatmap.pdf
2. Binary Regulon Activity Heatmap
step4则生成多个热图Step4_BinaryRegulonActivity_Heatmap*
3. Cell-type specific regulators (RSS)
当细胞种类比较多时可以用RSS来识别细胞类型特异性regulons。
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo[colnames(regulonAUC), "CellType"], )
rssPlot <- plotRSS(rss)
rssPlot$plot
得到热点图，颜色深浅代表z-score值，点的大小代表RSS评分。
然后就可以挑选各个细胞类型代表性regulon绘制热图等等啦。
注：帖子里面大家给出的建议是:根据regulon list提取到norm的AUCscore之后，再scale（计算z-score）比直接计算所有regulons的z-score再提取，绘制效果要好。
