稀疏的偏相关基因表达与海绵

马库斯列表,Azim Dehghani Amirabad,丹尼斯Kostka,马塞尔·h·舒尔茨

2022-04-10

目的

海绵的标志

海绵是第一个方法解决计算要求的任务报告重要的龙头、有效地在全基因组范围内的相互作用。除了龙头,这种方法非常适合推断出其他类型的监管如转录因子调节的相互作用。

介绍

小分子核糖核酸(microrna)小19 - 22日核苷酸长分子,促进退化的信使核糖核酸(mRNA)记录目标通过匹配种子序列。内源性RNA(龙头)竞争假说表明,拥有mrna结合位点的microrna在竞争。这存在的所谓的海绵,即。,genes that exert strong regulatory control via miRNA binding in a ceRNA interaction network. It is currently an unsolved problem how to estimate miRNA-mediated ceRNA interactions genome-wide. The most widely used approach considers miRNA and mRNA expression jointly measured for the same cell state. Several partial association methods have been proposed for determining ceRNA interaction strength using条件互信息偏相关为例。

然而,我们发现三个关键现有方法的局限性,防止建设一个精确的全基因组龙头、交互网络:(i)没有一个现有的方法考虑几个microrna的组合效应;(2)由于计算需求,所有公认的龙头、交互的推理gene-miRNA-gene交互在人类基因组中是禁止的;(3)一个有效的策略来确定推断龙头、互动的意义缺失,因此系统的重要参数是被忽视的。

为了克服这些挑战,我们开发了一个新颖的方法称为稀疏的偏相关基因表达(海绵)。我们减少了计算复杂度的构造一个全基因组龙头、交互网络几个步骤。首先,我们只考虑miRNA-gene交互或预测的实验验证。第二,我们仅保留miRNA-gene交互正规化回归模型系数为负。第三,而不是每个gene-miRNA-gene三个一组,我们计算一个敏感性的相关性(相关-偏相关)对于每一对基因基因给定所有共享microrna通上述滤波器作为假定的监管机构。最后,我们得到第一个数学公式来模拟不同参数的过程的零分布的系统:许多microrna基因之间的相关性和样本数量。我们的配方使经验假定值的计算统计的一个非常有效的方式,一个数量级的速度比先前的方法。分析表明,先前的研究已经低估了这些参数的影响在他们的网络推理。网络中心措施可以应用到海绵推断龙头、网络小说揭示已知和海绵,其中许多是潜在的生物标记物。

海绵进一步的细节展示如何提高对艺术的状态以及海绵推断电抗器,网络可以用于生物标志物的发现将会在我们的论文(包括提交的手稿,链接将在稍后的点)。

一般工作流程

海绵工作流

海绵工作流的概述。(A)预测和/或实验验证gene-miRNA交互进行正规化回归基因和microrna的表达数据。与负相互作用系数microrna的监管机构被保留,因为它们表明microrna的诱导基因表达的抑制作用。(B)计算灵敏度相关系数基于共享microrna基因对中标识(A)。(C)给定的样本数量,我们计算经验零模型对各种基因基因相关系数(k,列)和microrna (m,行)的数量。敏感性相关系数被分配到最匹配零模型和假定值是被推断出来的。(D)多个测试校正后,重要的电抗器,可以用来构造一个全基因组的相互作用,疾病或组织龙头、交互网络。在接下来的实践教程中,我们将重点介绍如何实现这些步骤与海绵R包。

(一)gene-miRNA交互

我们开始加载包及其依赖项

图书馆(海绵)

海绵有一个很小的例子和microrna的表达基因数据集用于说明功能。我们建议获得较大的表达数据集通过GEOquery或TCGAbiolinks R包,例如。加载包后,可以访问数据集的例子:

基因表达:

(gene_expr)
科大 FBXO32 TMEM132A CTTNBP2 MITF RANBP3L DNMT3A NTN4
TCGA-CC-A7IH 9.763889 8.664702 8.143475 2.004897 5.449300 5.7478519 9.482231 8.355215
TCGA-DD-A3A9 4.186326 5.754722 9.618885 2.150949 5.330171 0.7554858 9.985951 6.665077
TCGA-CC-A5UE 2.394404 5.514267 5.041414 1.045723 5.688147 3.1671973 9.742771 8.601548
TCGA-BC-A3KF 3.378858 5.423568 7.644005 7.085019 5.983949 1.1212144 8.863079 8.595803
TCGA-DD-A3A6 8.016376 10.199888 9.349632 3.828246 8.824847 5.0691178 7.681065 8.423871

microrna的表达:

(mir_expr)
MIMAT0000070 MIMAT0000680 MIMAT0000222 MIMAT0000460 MIMAT0000243
TCGA-CC-A7IH 8.110843 7.637696 14.397147 13.337312 17.77343
TCGA-DD-A3A9 9.425956 7.141432 14.627978 13.753494 13.84755
TCGA-CC-A5UE 10.017687 8.191404 16.030031 14.579990 16.12510
TCGA-BC-A3KF 11.260798 7.721752 17.036269 15.497766 16.15116
TCGA-DD-A3A6 7.717456 7.700341 7.605197 7.191296 13.79406

注意,预期格式的表达矩阵行和列对应的基因/ microrna与样品相对应。请确保命令是完全相同的矩阵和样本,每个样本有一个入口在两个矩阵(成对基因和microrna的表达数据)。

海绵使用双层过滤的方法来识别gene-miRNA交互。

首先,用户可以结合先验知识来缩小gene-miRNA交互的数量。例如,这可能的输出序列microrna的目标预测等软件TargetScanmiRcode或从数据库等获得的实验证据miRTarBaseLncBase。用户免费使用他或她的信任资源(s)和只需要确保miRNA-gene交互的海绵预计:

海绵允许任意数量的miRNA-gene交互矩阵我\ \ ()需要格式化的基因\ (G \)行和microrna\ \(米)在列。为每个元素\ (i_ {g、m}我\ \),我们认为如果miRNA-gene交互\ (i_ {g、m} > 0 \)。而海绵仅检查是否一个条目是零,这种方法允许各自的数量等附加信息结合位点或得分矩阵的编码后的分析。我们的示例包括TargetScan:

(targetscan_symbol)
MIMAT0000062 MIMAT0000063 MIMAT0000064 MIMAT0000065 MIMAT0000066
A1BG 0 0 0 0 0
A1CF 1 1 1 1 1
A4GNT 0 0 0 0 0
aac格式 0 0 0 0 0
AAED1 0 0 0 0 0

注意,基因和microrna的标识符的交互矩阵需要相同的类型和格式的表达式中使用矩阵。否则他们不能匹配的海绵。

海绵将用户提供的矩阵作为总结如上所示的条目和随后标识gene-miRNA交互与非零项候选人。在第二步中,我们使用稀疏正规化的回归与基因表达的依赖和microrna的表达作为解释变量。更具体地说,我们使用elasticnet(通过R包glmnet)套索和岭回归平衡。microrna的套索驱动器系数与微不足道的对基因表达的影响为零,而岭回归阻止我们踢出相关的microrna有类似的影响。Elasticnet有两个参数,\α(\ \)是用来平衡套索和岭回归的影响。\λ(\ \)是正则化参数。我们选择最优参数弹性网络如下:\ \(α= (0.1,0.2,…,1.0)\)我们执行10倍交叉验证来获得最优参数\λ(\ \)。最好的\α(\ \)然后选择这样的残差平方和误差最小。我们只是感兴趣的microrna抑制基因表达。因此,对于每一个基因,我们只保留那些有负面的microrna互动伙伴系数。而不是检查如果系数是小于零,海绵还允许更严格的阈值设置。缺省值是-0.05。另外,海绵提供了可能性进行野生识别重要的microrna。有关详细信息,请参阅文档。

我们应用步骤海绵的工作流程如下:

genes_miRNA_candidates < -sponge_gene_miRNA_interaction_filter(gene_expr =gene_expr,mir_expr =mir_expr,mir_predicted_targets =targetscan_symbol)

结果是一个基因列表,每个条目包含一个数据帧与保留microrna及其系数(或假定值,以防使用野生)。注意,mir_predicted_targets参数也可能是零。在这种情况下,没有先验知识和所有公认的miRNA-gene交互用于elasticnet回归。

genes_miRNA_candidates [1:2]
# # $ ASAP2 # # 1 # # microrna系数MIMAT0000074 CBR4美元-0.2459728 # # # # # # NULL

(B)电抗器,交互

第二步的海绵工作流,我们确定龙头、行动的候选人。海绵从(A)使用信息,并检查每一对基因如果他们共享一个或多个microrna与监管的效果。这减少了公认的龙头、交互的数量\ (n \选择2 \)考虑显著。接下来,海绵使用ppcor R包来计算多个microrna的敏感性相关()值。注意,这是一个泛化的敏感性相关定义的奶嘴等。。这些值捕获的共同贡献几个microrna的龙头、监管两个基因同时占他们的互相关:

\[开始\ {eqnarray} \标签{formula_scor} mscor (g_1里面,g_2, M) & = &软木(g_1里面,g_2) - pc机(g_1里面,g_2 |米)\ {eqnarray}结束\]

在哪里\ (M = {1,…,m_i} \)我\ \ ()之间共享microrna的数量\ (g_1里面\)\ (g_2 \)

我们应用步骤B的海绵工作流程如下:

ceRNA_interactions < -海绵(gene_expr =gene_expr,mir_expr =mir_expr,mir_interactions =genes_miRNA_candidates)

结果是一个数据表与基因(a和B), microrna的数量(),基因基因相关性(),偏相关microrna()和的值,即:\ (mscor =和pc机\):

(ceRNA_interactions)
日内瓦 geneB df 天哪 pc机 mscor
FBXO32 CTTNBP2 1 0.3036081 0.2268930 0.0767151
FBXO32 MITF 3 0.1595507 0.0697054 0.0898452
FBXO32 DCLK2 2 0.2046023 0.1147868 0.0898155
FBXO32 ASAP2 1 0.2666062 0.2118974 0.0547088
FBXO32 SRGAP1 1 0.1223776 0.1237864 -0.0014088
FBXO32 LIMCH1 2 0.1854762 0.1264297 0.0590465

(C) Null-model-based假定值计算

值的随机分布的理论并不存在和先前的研究使用灵敏度相关使用简单的阈值或随机背景分布推断从随机的三胞胎。仿真实验表明,这种方法不考虑事实的分布参数强烈影响基因基因相关性和样本的数量被认为。

制定一项战略,这让我们推断值的分布在零假设下,microrna不会影响两个基因之间的相关性,即。\ (mscor = 0 \)。这是一个基础的随机采样协方差矩阵方法在这个零假设我们控制我)基因基因相关性和ii)的microrna的数量对应于协方差矩阵的列数。有关详细信息,请参阅[额外计算协方差矩阵]如何构造这些矩阵。因为这个抽样程序计算量,我们包括预先计算的集协方差矩阵的参数组合,即10协方差矩阵为每一对基因基因相关性\(软木= (0.2,……,0.9)\)和microrna的数量\ (m = (1、2、…, 8) \)。使用这些协方差矩阵和包mvrnorm允许我们样本值的随机分布对可用的数量的样品(我们表达矩阵的行数):

mscor_null_model < -sponge_build_null_model(number_of_datasets =One hundred.,number_of_samples =nrow(gene_expr))

我们可以画出这些模拟的结果的分布是受基因基因相关性和microrna的数量(如果ggplot2安装):

sponge_plot_simulation_results(mscor_null_model)

显然,这是有点粗糙,因为我们只有样品100数据点为每个分区保持装饰图案的编译时检查。默认的数据集采样是muche更合理的价值1 e6。

使用上面的零模型、海绵推断实证假定值的系数在步骤(B)。这是通过第一选择最匹配零模型,例如最近的匹配和许多microrna基因基因相关性。在这个分区的随机分布,我们现在排名我们计算的值。这个排名可以用来计算实证假设机率。注意,精度,即最小的假定值获得数据集采样的数量是有限的。例如,100年抽样数据集是这里最好的假定值我们可以是手段\ (p = 0.01 \)。如果我们样品1 e6的数据集,最好的精度提高,假定值\ (p = 1 e-6 \)

后计算假定值系数,我们应用多个测试修正所有的值在一个分区中,占空模型基因基因相关性较低这一事实将经常往往大大超过那些高基因基因相关性。推断出我们所做的假定值如下:

ceRNA_interactions_sign < -sponge_compute_p_values(sponge_result =ceRNA_interactions,null_model =mscor_null_model)

上面的方法将两个附加列添加到龙头的结果,即\ (p.val \)\ (p.adj \)

(ceRNA_interactions_sign)
日内瓦 geneB df 天哪 pc机 mscor p.val p.adj
FBXO32 SRGAP1 1 0.1223776 0.1237864 -0.0014088 0.54 0.540
SRGAP1 ESR1 1 0.1234821 -0.0011236 0.1246057 0.01 0.020
FBXO32 CTTNBP2 1 0.3036081 0.2268930 0.0767151 0.02 0.070
FBXO32 ASAP2 1 0.2666062 0.2118974 0.0547088 0.06 0.098
TMEM132A TET1 1 0.3128573 0.2308158 0.0820415 0.02 0.070
CTTNBP2 DCLK2 1 0.3186208 0.2666053 0.0520155 0.07 0.098

(D)龙头、交互网络

最后一步是相当直接的。我们决定一个罗斯福截止选择重要的龙头、交互,例如罗斯福< 0.01和提取的结果。

ceRNA_interactions_fdr < -ceRNA_interactions_sign [哪一个(ceRNA_interactions_sign美元p.adj<0.2),)
日内瓦 geneB df 天哪 pc机 mscor p.val p.adj
2 SRGAP1 ESR1 1 0.1234821 -0.0011236 0.1246057 0.01 0.020
3 FBXO32 CTTNBP2 1 0.3036081 0.2268930 0.0767151 0.02 0.070
4 FBXO32 ASAP2 1 0.2666062 0.2118974 0.0547088 0.06 0.098
5 TMEM132A TET1 1 0.3128573 0.2308158 0.0820415 0.02 0.070
6 CTTNBP2 DCLK2 1 0.3186208 0.2666053 0.0520155 0.07 0.098
9 ASAP2 LIMCH1 1 0.2887748 0.2370206 0.0517542 0.07 0.098

由此产生的网络也可以绘制

sponge_plot_network(ceRNA_interactions_fdr genes_miRNA_candidates)
# #加载需要名称空间:visNetwork

网络分析

重要的龙头、交互从上面可以用来构造一个基因(龙头)的龙头、网络节点和边是重要的龙头、交互时预测。海绵构造这样一个网络使用igraph R包和计算单个基因的不同的网络统计值。特征向量中间状态中心可以用来评估基因在电抗器,网络的重要性。在这个示例中,我们只选择很少基因使得这个练习有点徒劳的。然而,如果应用到更大的基因集,中间性中心显示强大的龙头,可能造成生物标志物的microrna的监管。

network_centralities < -sponge_node_centralities(ceRNA_interactions_fdr)

请注意,关于龙头效应强度的附加信息可以通过加权特征向量合并或betweeness中心。在这种情况下,我们添加一个列权重结果表如下:

ceRNA_interactions_fdr_weight < -ceRNA_interactions_fdrceRNA_interactions_fdr_weight美元重量< -- - - - - -log10(ceRNA_interactions_fdr美元p.adj)weighted_network_centralities < -sponge_node_centralities(ceRNA_interactions_fdr)

我们还可以绘制这些结果与上一个节点度分布和中心与节点下面的程度。参数前控制多少前x样品标签。

sponge_plot_network_centralities(weighted_network_centralities顶级=1)

您也可以选择个人情节

sponge_plot_network_centralities(weighted_network_centralities测量=“顺便说一句”,顶级=1)

并行化和共享内存

虽然海绵很快在推断龙头、交互仍然需要几天处理单个CPU与几百的样本数据集。我们因此建立海绵的foreach R包:

foreach R包将允许您在多个cpu并行处理海绵。这种方法的好处是,它是由用户来决定使用哪个后端。我们把foreach方案细节。的后端用于foreach doParallel R包。一旦建立了并行计算环境中,海绵可以使用它自动(不需要额外参数):

图书馆(doParallel)图书馆(foreach)num.of。天哪es <-2#更多的计算集群上#如果你想使用日志记录#日志记录。文件< -"where_my_log_file_should_go.log"日志记录。文件< -cl < -makeCluster(num.of.cores输出文件=logging.file)registerDoParallel(cl)genes_miRNA_candidates < -sponge_gene_miRNA_interaction_filter(gene_expr =gene_expr,mir_expr =mir_expr,mir_predicted_targets =targetscan_symbol)ceRNA_interactions < -海绵(gene_expr =gene_expr,mir_expr =mir_expr,mir_interactions =genes_miRNA_candidates)stopCluster(cl)

请注意,如果你想使用BiocParallel包只需注册任何并行foreach端如上面使用如下:

图书馆(“BiocParallel”)注册(DoparParam(),默认=真正的)

如果你有访问集群计算,你可以注册它作为并行端使用clustermq包。有关详细信息,请参阅clustermq包。

图书馆(clustermq)register_dopar_cmq(n_jobs =2)

经常遇到警告的foreach就是每一个平行的工人需要复制的数据操作。大型数据集这成为一个问题甚至与大内存计算服务器。我们因此支持bigmemory R包允许每个并行工作进程/线程访问共享内存中的数据。利用大内存包只是gene_expr和mir_expr参数替换bigmemory描述对象。我们把文档bigmemory包的详细信息。如果你有访问计算服务器与奖金的核心,通常你可以计算全基因组电抗器,网络在几个小时内消失。

计算协方差矩阵为零模型

海绵允许采样协方差矩阵得到一个和一个特定数量的microrna基因基因相关性。预先计算的矩阵的包包含许多不同参数组合。然而,它也遭遇的计算额外的矩阵:

more_covariance_matrices < -sample_zero_mscor_cov(m =1,number_of_solutions =10,gene_gene_correlation =0.5)

这些可以用来模拟数据样本的随机分布系数给定样本中可用的数量表达式的数据。

mscor_coefficients < -sample_zero_mscor_data(cov_matrices =more_covariance_matrices,number_of_samples =200年,number_of_datasets =One hundred.)(unlist(mscor_coefficients),主要=“mscor系数的随机分布”)

更广泛的应用

海绵主要是用来推断龙头、交互效率。我们指出,然而,海绵,探讨出部分关联背后的原则是适用于其他类型的调节分子生物学。海绵并不局限于消极监管microrna的监管的情况,也可以是有用的学习,例如,监管RNA结合蛋白或转录因子调控。此外,我们可以预见,这种方法也适用于生物信息学之外感兴趣的部分协会在哪里。据我们所知,海绵了第一个方法随机样本协方差矩阵的零假设下偏相关,相关性是相等的,这可能是有用的作为起点的理论更深入的了解在未来的局部相关性属性。

引用

列表中,M。,Dehghani Amirabad, A., Kostka, D., & Schulz, M. H. (2019). Large-scale inference of competing endogenous RNA networks with sparse partial correlation. Bioinformatics, 35(14), i596-i604.doi: 10.1093 /生物信息学/ btz314