R包SCENIC：分析流程

2019-05-10 本文已影响4人阿糖胞苷_SYSU

GRN：基因调控网络

分为四步进行GRN分析

1.GENIE3/GRNBoost：基于共表达情况鉴定每个TF的潜在靶点;

过滤表达矩阵并运行R包GENIE3/GRNBoost；

将得到的靶点格式化为共表达模块；

2.RcisTarget：基于DNA-motif 分析选择潜在的直接结合靶点（调控原件）

鉴定细胞状态和细胞的调控分子

3. AUCell：分析每个细胞的网络活动度

计算AUC对细胞的调控元件进行评分；

将网络活动度转为二分类矩阵（on、off）

4.细胞聚类：基于GRN的活动度鉴定稳定的细胞状态并对结果进行探索

一.input数据准备：

exprMatPath <- "data/sceMouseBrain.RData"

library(SingleCellExperiment)

sceMouseBrain <- readRDS(exprMatPath)

exprMat <- counts(sceMouseBrain)

dim(exprMat)

提供cell info或者表型相关的信息,用于step3~4的比较分析

cellInfo <- get_cellAnnotation(loom)

close_loom(loom)# cellInfo$nGene <- colSums(exprMat>0)# cellInfo <- data.frame(cellInfo)head(cellInfo)

dir.create("int")

saveRDS(cellInfo, file="int/cellInfo.Rds")# Color to assign to the variables (same format as for NMF::aheatmap)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))

二.初始化SCENIC的设置参数

org，dbDir，dbs，datasetTitle

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)

根据需要进行修正

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"

saveRDS(scenicOptions, file="int/scenicOptions.Rds")

三.共表达网络

1.利用R包GENIE3/GRNBoost推断潜在转录因子靶点

input数据形式：表达矩阵（过滤）或者转录因子的列表

output数据形式：相关矩阵，用于构建共表达模块（runSCENIC_1_coexNetwork2modules()）

关于R包GENIE3/GRNBoost的选择：

GENIE3：推荐3-5k细胞数据集，耗时较长

GRNBoost：较大的数据集，耗时稍短

细胞的二次抽样：

如果input的数据集存在较多低质量的细胞或者耗时较长，可以选择对数据集的细胞先进行二次抽样（可以是随机抽样或者是选高质量细胞），作为共表达的训练集，用于AUCell所有细胞评估。注意二次抽样的细胞必须具有代表性。

基因的过滤、筛选：推荐较松的过滤标准，移除低表达或者少数细胞表达的基因

（1）total number of reads per gene 平均每个基因的读序数

目的：是移除疑似“噪声”的基因

默认保存所有样本≥6UMI counts的基因，根据数据集的表达量单位（UMI,TPM等)校正（minCountsPerGene）。

（2） number of cells in which the gene is detected表达该基因的细胞数

目的：移除少量“噪声”细胞，即该基因在少数细胞高表达，导致噪声产生

默认保存至少在1%的细胞表达的基因(minSamples)。

为避免移除有生物学意义的罕见细胞，推荐设置低于细胞最小分群占比的百分数作为筛选标

准

（3）仅保留RcisTarget包中数据库有的基因

目的：移除下游分析不需要的基因，节省GENIE3/GRNBoost运行的时间

# (Adjust minimum values according to your dataset)genesKept <- geneFiltering(exprMat, scenicOptions=scenicOptions,

minCountsPerGene=3*.01*ncol(exprMat),

minSamples=ncol(exprMat)*.01)

## Maximum value in the expression matrix: 421

## Ratio of detected vs non-detected: 0.24

## Number of counts (in the dataset units) per gene:

## Min. 1st Qu. Median Mean 3rd Qu. Max.

## 0.0 17.0 48.0 146.8 137.0 5868.0

## Number of cells in which each gene is detected:

## Min. 1st Qu. Median Mean 3rd Qu. Max.

## 0.00 10.00 25.00 39.01 60.00 180.00

##

## Number of genes left after applying the following filters (sequential):

## 771 genes with counts per gene > 6

## 770 genes detected in more than 2 cells

## 770 genes available in RcisTarget database

## Gene list saved in int/1.1_genesKept.Rds

在进行共表达网络推断之前，先检查是否把相关的TF筛除了，以判断过滤标准是否设置正确。

interestingGenes <- c("Sox9", "Sox10", "Dlx5")

# any missing?

interestingGenes[which(!interestingGenes %in% genesKept)]

## character(0)

从exprMat取过滤后的基因的部分矩阵作为过滤后的矩阵exprMat_filtered

exprMat_filtered <- exprMat[genesKept, ]

dim(exprMat_filtered)

## [1] 770 200

移除exprMat

rm(exprMat)

计算spearman相关系数

runCorrelation(exprMat_filtered, scenicOptions)

推断TF靶点

函数 runGenie3 以默认参数运行r包GENIE3，将RcisTarget的数据库作为候选的调控因子，对于大多数数据集是足够的。

GENIE3运用了随机森林的算法，每次运行的结果都有些许不同，ntrees越多，可变性就越低。推荐设置set.seed从而产生较为准确的结果。因为这一步耗时较长（数小时或数天），推荐可以在另外的r console或者r server或者HPC进行下面的分析。

# log转换（如果未进行归一化或者log转换）

exprMat_filtered <- log2(exprMat_filtered+1)

# 运行GENIE3

runGenie3(exprMat_filtered, scenicOptions)

四.构建并对GRN进行评分

构建基因调控网络：

1.获得共表达模块

2.r 包RcisTarget进行TFmotif分析获得调控元件

鉴定细胞状态：

3.r包AUCell对细胞的GRN进行评分

4.根据GRN活动度对细胞进行聚类

loom <- open_loom(loomPath, mode="r")

exprMat <- get_dgem(loom)

close_loom(loom)

# Optional: log expression (for TF expression plot, it does not affect any other calculation)

logMat <- log2(exprMat+1)

dim(exprMat)

scenicOptions@settings$verbose <- TRUE

scenicOptions@settings$nCores <- 10

scenicOptions@settings$seed <- 123

#用的测试数据

scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"]

runSCENIC_1_coexNetwork2modules(scenicOptions)

#用top5perTarget进行测试

runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget"))

runSCENIC_3_scoreCells(scenicOptions, logMat)

对网络活动度进行AUC评分，对细胞进行二分类

output是一个二进制矩阵，行名代表调控元件，列名代表细胞，“1”代表激活调控元件，“0”代表其他。r包AUCell可以自动确定AUCscore的阈值，但是一般比较保守，建议人为进行调整从而选择最佳阈值。

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)

五.根据GRN评分对细胞进行聚类

聚类方法可以是tsne，umap或者其他HC方法均可；建议聚类时尝试不同参数看细胞状态是否稳定；

#选择PC的数目

nPcs <- c(5,15,50)

scenicOptions@settings$seed <- 123

#用不同参数进行t-SNE降维聚类

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")

#保存为PDF文件

fileNames <- paste0("int/",grep(".Rds", grep("tSNE_", list.files("int"), value=T), value=T))

#对t-SNE结果进行可视化，并进行比较

par(mfrow=c(length(nPcs), 3))

fileNames <- paste0("int/",grep(".Rds", grep("tSNE_AUC", list.files("int"), value=T, perl = T), value=T))

plotTsne_compareSettings(fileNames, scenicOptions, showLegend=FALSE, varName="CellType", cex=.5)

#选择高度可信的调控元件

par(mfrow=c(3,3))

fileNames <- paste0("int/",grep(".Rds", grep("tSNE_oHC_AUC", list.files("int"), value=T, perl = T), value=T))

plotTsne_compareSettings(fileNames, scenicOptions, showLegend=FALSE, varName="CellType", cex=.5)

上述选择的t-SNE用于绘图，并保存为.rds文件

scenicOptions@settings$defaultTsne$aucType <- "AUC"

scenicOptions@settings$defaultTsne$dims <- 5

scenicOptions@settings$defaultTsne$perpl <- 15

saveRDS(scenicOptions, file="int/scenicOptions.Rds")

六.输出为.loom文件

需要r包：SCopeLoomR；运用函数export2scope()

# Export:

scenicOptions@fileNames$output["loomFile",] <- "output/mouseBrain_SCENIC.loom"

export2scope(scenicOptions, exprMat)

从.loom读取结果：regulons, AUC, embeddings 这3个结果；motif富集分析和共表达模块分析的结果以txt.file形式另外保存。

library(SCopeLoomR)

scenicLoomPath <- getOutName(scenicOptions, "loomFile")

loom <- open_loom(scenicLoomPath)# Read information from loom file:regulons_incidMat <- get_regulons(loom)

regulons <- regulonsToGeneLists(regulons_incidMat)

regulonsAUC <- get_regulonsAuc(loom)

regulonsAucThresholds <- get_regulonThresholds(loom)

embeddings <- get_embeddings(loom)