【SCENIC】安装测试1（R版本）

2023-09-23 本文已影响0人 jjjscuedu

今天开始测试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。