regsplice 1.25.0
的regsplice包实现了用于检测RNA测序(RNA-seq)和外显子微阵列数据集中差异外显子使用(差异剪接)的统计方法。
的regsplice方法是基于使用套索(l1 -正则化)来提高标准广义线性模型的功率。的一个主要优势regsplice就是运行时比其他主流方法快。我们预计类似的基于正则化的方法也可能在其他设置中有应用。
详细的统计方法和与其他方法的性能比较将在接下来的一篇论文中描述。
的示例工作流regsplice软件包，使用小型模拟RNA-seq数据集。
跑步有两种选择regsplice:您可以使用包装器功能一步运行一个完整的工作流regsplice ()；或者您可以按顺序为每个步骤运行单个函数，这提供了额外的灵活性和对方法的深入了解。下面演示了这两个选项。
用于示例工作流的数据集包括来自模拟人类RNA-seq数据集的100个基因子集的外显子级读取计数，该数据集由6个生物样本组成，2种条件下各有3个样本。
原始数据集来自论文:
Soneson, Matthes等人(2016)，Isoform预过滤提高了基于计数的分析差异转录本使用方法的性能，基因组生物学，可以在这里
本文中的原始数据文件包含模拟的RNA-seq读取(FASTQ和BAM文件)，可以从ArrayExpress的登录代码中获得e - mtab - 3766．
类提供的Python计数脚本生成外显子bin计数DEXSeq包装，使用选项从重叠基因中排除外显子，而不是将它们聚集到多基因复合体中(参见Soneson et al. 2016，补充材料)。
对于这个示例工作流，我们从这个模拟数据集中选择了一个由前100个基因组成的子集。这100个基因的外显子级读取计数和真正的差异剪接状态标签保存在文本文件中vignette_counts.txt而且vignette_truth.txt在extdata /目录中的regsplice包源代码。
的regsplice方法设计用于RNA-seq读取计数和外显子微阵列强度。
如果使用外显子微阵列数据，工作流程中的主要步骤与RNA-seq数据相同，如下图所示。但是，需要对工作流程进行以下调整:
而不是RNA-seq读取计数，外显子微阵列强度的矩阵或数据帧提供给计数输入参数。这个参数的名字是静止的计数，不管输入数据类型是什么。
外显子微阵列的强度在被提供之前，应该在外部进行log2转换regsplice．这通常在微阵列数据的预处理过程中完成，也可能根据你的软件自动完成。
应该通过设置参数禁用零计数和低计数外显子箱的过滤filter_zero = FALSE而且filter_low_counts = FALSE．
应该通过设置禁用归一化因子的计算normalize = FALSE．
计算limma-voom应该通过设置禁用转换和权重voom = FALSE．
加载小插图示例数据文件，其中包含6个生物样本中100个基因的模拟RNA-seq读取计数。从原始数据中，提取计数表(计数)、基因id (gene_IDs)，以及每个基因的外显子仓数量(n_exons）.
然后创建条件向量，它指定了每个生物样本的实验条件或处理组。
# load data file_counts <- system.file("extdata/vignette_counts.txt"， package = "regsplice") data <- read。表(file_counts头= TRUE, 9 =“t \”,stringsAsFactors = FALSE)头(数据)# #外显子sample1 sample2 sample3 sample4 sample5 sample6 ENSG00000000003:001 # # 576 506 526 643 482 826 # # 2 ENSG00000000003:002 141 122 126 157 121 191 # # 3 ENSG00000000003:003 123 102 106 133 99 156 # # 4 ENSG00000000003:004 86 76 77 98 72 112 # # 5 ENSG00000000003:005 97 83 87 113 76 126 # # 6 ENSG00000000003:006 133 107 116 155 97 170 #提取计数,gene_IDs, n_exons计数(< -数据,2:7] tbl_exons <- table(sapply(strsplit(data$exon， ":")， function(s) s[[1]])) gene_IDs <- names(tbl_exons) n_exons <- unname(tbl_exons) dim(counts) ## [1] 3191 6 length(gene_IDs) ## [1] 100 head(gene_IDs) ## [1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ## [5] "ENSG00000000460" "ENSG00000000938" length(n_exons) ## [1] 100 sum(n_exons) ## [1] 3191 # create condition vector condition <- rep(c("未处理"，"已处理")，每个= 3)条件##[1]“未处理”“未处理”“未处理”“已处理”“已处理”“已处理”的regsplice ()包装器函数用一个命令运行分析管道。包装器函数所需的输入格式是RegspliceData对象创建的RegspliceData ()构造函数。
研究的结果regsplice分析包括一组多重测试调整的p值(Benjamini-Hochberg错误发现率，FDR)，量化每个基因的差异外显子使用(DEU)的统计证据。调整后的p值被用来根据DEU的证据对数据集中的基因进行排序，并且可以指定一个显著性阈值来生成具有统计显著证据的DEU的基因列表。
包装器函数返回基因名称、拟合模型结果、原始p值、多次测试调整的p值、似然比(LR)测试统计信息以及LR测试的自由度。
所需的输入RegspliceData ()构造函数为计数(RNA-seq读取计数或外显子微阵列强度的矩阵或数据帧)，gene_IDs(基因id向量)，n_exons(外显子长度向量，即每个基因的外显子箱数)和条件(每个生物样本的实验条件向量)。
或者，输入可以作为SummarizedExperiment对象，该对象将被解析以提取这些组件。这在运行时可能很有用regsplice作为与其他包连接的管道的一部分。
注意，我们已经使用了suppressWarnings ()为了隐藏与该示例数据集中每个基因的少量观察结果相关的警告信息。对于大多数数据集，不应该出现这些警告消息。
看到RegspliceData ?而且regsplice ?有关其他详细信息，包括其他可用输入和选项。的种子参数用于生成可重复的结果。
库(regsplice) rs_data <- RegspliceData(counts, gene_IDs, n_exons, condition) rs_results <- suppressWarnings(regsplice(rs_data, seed = 123)) ##删除936个外显子(s)，零计数##删除1个剩余单外显子基因##删除240个低计数外显子(s) ##删除4个剩余单外显子基因##装正则化(lasso)模型…##拟合空模型…这个函数summaryTable ()用于生成结果的汇总表。
结果显示为顶部的数据帧n大多数高显著性基因，根据错误发现率(FDR)或原始p值排序，直到指定的显著性阈值(例如FDR < 0.05)。
这个论点rank_by选择是否按FDR或原始p值排序。
要显示达到显著性阈值的所有基因的结果，请设置参数n = Inf．若要显示数据集中所有基因的结果，请同时设置两者n = Inf而且阈值= 1．
有关更多细节，请参见summaryTable ?．
摘要表(rs_results) ## # gene_IDs p_vals p_adj LR_stats df_tests ## 1 ENSG00000004766 1.963457e-17 6.086718e-16 137.05717 25 ## 2 ENSG00000001461 8.635485e-06 1.338500e-04 26.20603 3 ## 3 ENSG00000003436 1.316879e-04 1.360775e-03 17.87015 2 . ENSG00000004766 8.663457e -17 6.086718e-16 137.05717 25或者,regsplice分析管道可以使用每个步骤的单独函数运行，这提供了额外的灵活性和对统计方法的洞察。下面描述这些步骤。
第一步是创建一个RegspliceData对象中函数所需的格式的数据regsplice分析管道。
的RegspliceData格式是基于SummarizedExperimentBioconductor的容器。这种格式的主要优点是，子设置操作自动保持行和列的数据和元数据同步，这有助于避免由于选择不正确的行或列索引而导致的错误。
最初,RegspliceData对象包含原始数据以及行(基因和外显子箱)和列(生物样本)的元数据。的后续步骤中regsplice分析管道，修改数据值，并添加额外的数据和元数据。最终结果存储在RegspliceResults对象。
所需的输入是计数(RNA-seq读取计数或外显子微阵列强度的矩阵或数据帧)，gene_IDs(基因id向量)，n_exons(外显子长度向量，即每个基因的外显子箱数)和条件(每个生物样本的实验条件向量)。
或者，输入可以作为SummarizedExperiment对象，该对象将被解析以提取这些组件。这在运行时可能很有用regsplice作为与其他包连接的管道的一部分。
注意，警告消息是由于在这个小插图中使用的数据集中，每个基因的观察数量很少。在大多数数据集中，这些警告消息不应该出现。
有关更多细节，请参见RegspliceData ?．
库(regsplice) rs_data <- RegspliceData(计数，gene_IDs, n_exons，条件)接下来，使用函数filterZeros ()过滤所有生物样品(列)中计数为零的外显子箱(行)。
任何剩余的单外显子基因也被移除(因为差异剪接需要多个外显子)。
如果您使用的是外显子微阵列数据，则应跳过此步骤。
有关更多细节，请参见filterZeros ?．
rs_data <- filterZeros(rs_data) ## removed 936外显子(s) with zero counts ## removed 1剩1个单外显子基因(s)过滤低计数外显子箱用filterLowCounts ()．
的参数filter_min_per_exon而且filter_min_per_sample控制过滤量。提供了默认值;但是，这些应该根据样本总数和每个条件的样本数量进行调整。
任何剩余的单外显子基因也被移除。
如果您使用的是外显子微阵列数据，则应跳过此步骤。
有关更多细节，请参见filterLowCounts ?．
rs_data <- filterLowCounts(rs_data) ## removed 240 low-count exon(s) ## removed 4个剩余的单外显子基因(s)这个函数runNormalization ()计算用于缩放库大小的归一化因子。
默认情况下,runNormalization ()采用TMM (m值的修剪平均值)归一化方法(Robinson and Oshlack, 2010)，实现于刨边机包中。有关详细信息，请参阅相关文档calcNormFactors ()在刨边机包中。
这一步应该在过滤之后完成。然后使用归一化因子limma-voom在下一步。
如果您使用的是外显子微阵列数据，则应跳过此步骤。
有关更多细节，请参见runNormalization ?．
rs_data <- runNormalization(rs_data)下一步是使用limma-voom转换计数并计算外显子级权值。这是用runVoom ()函数。
的limma-voom方法将计数转换为每百万log2计数(logCPM)，并根据观察到的均值-方差关系计算外显子级别的权重。这是必需的，因为原始或对数转换计数不满足线性建模所需的统计假设(即方差相等)。后limma-voom经过变换和权值计算，可以采用线性建模方法。
更多细节，请参见下面的论文，其中介绍轰；或者是limma用户指南(“差分拼接”部分)可在Bioconductor上获得。
请注意,轰假设计数为零或低的外显子箱(行)已经被删除，因此这一步应该在过滤后执行filterZeros ()而且filterLowCounts ()．
如果标准化因子可用(从上一步与runNormalization ())，它们将被使用轰计算标准化库的大小。如果没有，轰将改为使用非规范化库大小(按列的总计数)。
如果您使用的是外显子微阵列数据，则应跳过此步骤。
有关更多细节，请参见runVoom ?．
rs_data <- runVoom(rs_data) #视图列元数据包括规范化因子和规范化库大小colData(rs_data) ## 6行4列的DataFrame ## sample_names条件norm_factors lib_sizes ## <字符> <字符> <数字> <数字> ## 1 sample1未处理0.872019 348678 ## 2 sample2未处理1.003102 327180 ## 3 sample3未处理1.082822 315879 ## 4 sample4处理1.033876 305823 ## 5 sample5处理0.983314 359353 ## 6 sample6处理1.038511 300455的initializeResults ()函数创建一个RegspliceResults对象，该对象将包含分析的结果。该对象将在后续步骤中填充。
有关更多细节，请参见initializeResults ?．
rs_data <- initializeResults(rs_data)模型拟合函数有三个:
fitRegMultiple ()拟合包含最优外显子子集的正则化(套索)模型:每个基因的条件相互作用项。模型拟合过程只考虑相互作用项，因此外显子和样品的主要效应项总是包含在内。这确保了空模型是嵌套的，允许计算似然比测试。
fitNullMultiple ()适合不包含任何交互项的空模型。
fitFullMultiple ()适合“完整”模型，其中包含每个基因的所有外显子:条件相互作用术语。
拟合函数拟合数据集中所有基因的模型。
注意，我们已经使用了suppressWarnings ()为了隐藏与该示例数据集中每个基因的少量观察结果相关的警告信息。对于大多数数据集，不应该出现这些警告消息。
有关更多细节，请参见fitRegMultiple ?，fitNullMultiple ?,或fitFullMultiple ?．
#拟合正则化模型rs_results <- suppressWarnings(fitRegMultiple(rs_results, rs_data, seed = seed)) ##拟合正则化(lasso)模型…# fit null models rs_results <- fitNullMultiple(rs_results, rs_data, seed = seed) ## fit null models…# fit "full" models(不需要，如果'when_null_selected = "ones"'在下一步)rs_results <- fitFullMultiple(rs_results, rs_data, seed = seed) ## fit full models…这个函数LRTests ()计算拟合模型与零模型之间的似然比(LR)检验。
如果拟合的正则化(套索)模型包含至少一个外显子:条件相互作用项，LR检验将套索模型与嵌套的空模型进行比较。但是，如果套索模型包含零项相互作用，则套索模型和零模型是相同的，因此LR检验无法计算。的when_null_selected参数允许用户选择在这些情况下做什么:要么将p-value设置为1 (When_null_selected = "ones"）;或者使用包含所有外显子的“完整”模型计算LR检验:条件相互作用项(When_null_selected = "full")，由于术语数量较多，这降低了功率，但允许区分这些基因之间的差异外显子使用的证据。你也可以返回这些基因的NAs (when_null_selected = "NA"）.
默认选项是When_null_selected = "ones"．这就简单地称所有这些基因都是不显著的，在大多数情况下这就足够了，因为我们更感兴趣的是有强烈证据证明差异外显子使用的基因。但是，如果对数据集中证据较低的基因进行排序很重要，则使用When_null_selected = "full"选择。如果When_null_selected = "ones"或when_null_selected = "NA"，“完整”安装的型号是不需要的。
结果对象包含基因名称、拟合模型结果、原始p值、多次测试调整的p值(Benjamini-Hochberg错误发现率，FDR)、似然比(LR)检验统计量和LR检验的自由度。
有关更多细节，请参见LRTests ?．
rs_results <- LRTests(rs_results)这个函数summaryTable ()用于生成结果的汇总表。
结果显示为顶部的数据帧n大多数高显著性基因，根据错误发现率(FDR)或原始p值排序，直到指定的显著性阈值(例如FDR < 0.05)。
这个论点rank_by选择是否按FDR或原始p值排序。
要显示达到显著性阈值的所有基因的结果，请设置参数n = Inf．若要显示数据集中所有基因的结果，请同时设置两者n = Inf而且阈值= 1．
有关更多细节，请参见summaryTable ?．
摘要表(rs_results) ## # gene_IDs p_vals p_adj LR_stats df_tests ## 1 ENSG00000004766 1.963457e-17 6.086718e-16 137.05717 25 ## 2 ENSG00000001461 8.635485e-06 1.338500e-04 26.20603 3 ## 3 ENSG00000003436 1.316879e-04 1.360775e-03 17.87015 2 . ENSG00000004766 8.663457e -17 6.086718e-16 137.05717 25对于这个小插图中的模拟数据集，每个基因的真实差异剪接状态是已知的。在本节中，我们将展示如何分析结果并计算列联表，显示真阳性、真阴性、假阳性和假阴性的数量。
如上面的工作流程所示，我们可以使用summaryTable ()带参数的函数n = Inf显示所有具有差异外显子使用(DEU)显著证据的基因列表。
summaryTable(rs_results, n = Inf) ## gene_IDs p_vals p_adj LR_stats df_tests ## 1 ENSG00000004766 1.963457e-17 6.086718e-16 137.05717 25 ## 2 ENSG00000001461 8.635485e-06 1.338500e-04 26.20603 3 ## 3 ENSG00000003436 1.316879e-04 1.360775e-03 17.87015 2 . ENSG00000004766 1.963457e-17 6.086718e-16 137.05717 25在给定的阈值下，具有显著证据的DEU基因的总数也可以计算出来。
注意，我们使用多重测试调整的p值(Benjamini-Hochberg错误发现率，fdr)进行计算。FDR < 0.05的标准阈值意味着列表中5%的基因可能是错误发现。
sum(p_adj(rs_results) < 0.05) ## [1] 3 table(p_adj(rs_results) < 0.05) ## ## FALSE TRUE ## 78如上所述，每个基因的真差异剪接(DS)状态是已知的，因为这是一个模拟数据集。因此，我们可以计算列联表，比较每个基因在给定显著性阈值下的真实和预测DS状态。提高显著性阈值返回更多的基因，代价是更多的假阳性。
# load true DS状态标签file_truth <- system.file("extdata/vignette_truth.txt"， package = "regsplice") data_truth <- read. txttable(file_truth, header = TRUE, sep = "\t"， stringsAsFactors = FALSE) str(data_truth) ## 'data.frame': 100 obs。## $ gene: chr "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"…## $ ds_status: int 0 0 0 0 0 0 0 0 1 0…#删除基因筛选在regsplice分析data_truth < - data_truth [data_truth基因% % gene_IDs美元(rs_results)]暗(data_truth) # # 81年[1]2长度(gene_IDs (rs_results)) # #[1] 81 #真正的DS基因数量在模拟数据集和(data_truth ds_status = = 1美元)# #表(data_truth ds_status美元)[1]6 0 1 # # 75 # # # # 6 #列联表比较真实和预测DS地位对于每个基因#(意义阈值:FDR < 0.05)表(true = data_truth$ds_status，预言= p_adj(rs_results) < 0.05) ##预言## true FALSE true ## 0 73 2 ## 1 5 1 #增加阈值检测到更多的基因，以牺牲更多的假阳性表(true = data_truth$ds_status，预言= p_adj(rs_results) < 0.99) ##预言##真假真## 0 69 6 ## 1 4 2上面的工作流程中没有讨论的其他用户选项包括:
α:弹性网参数glmnet模型拟合函数。的价值α必须介于0(脊回归)和1(套索)之间。默认值为1，符合套索模型。看到glmnet包文档以获取更多详细信息。lambda_choice:选择哪个最优λ值来选择cv.glmnet交叉验证适合。可用的选择有“lambda.min”(具有最小交叉验证误差的模型)和“lambda.1se”(最正则化的模型，交叉验证误差在一个最小标准误差内)。默认值为“lambda.min”．看到glmnet包文档以获取更多详细信息。有关进一步的详细信息，包括所有可用用户选项的完整列表和描述，请参阅regsplice ()包装器函数，可以通过regsplice ?或帮助(regsplice)．
这个函数createDesignMatrix ()为每个基因创建模型设计矩阵。该函数由模型拟合函数自动调用，因此不需要直接使用。在本节中，我们将演示它如何对单个基因起作用，并展示一个示例设计矩阵，以便进一步深入了解统计方法。
设计矩阵包括每个外显子和每个样品的主要效应项以及外显子与条件之间的相互作用项。
注意，设计矩阵不包括条件的主效应项，因为这些被吸收到样本的主效应项中。此外，设计矩阵不包括拦截列，因为让模型拟合函数稍后添加拦截项更简单。
有关更多细节，请参见createDesignMatrix ?．
# 3外显子基因# 4个生物样本2在每个样品2条件design_example < - createDesignMatrix(条件=代表(c(0, 1),每个= 2),n_exons = 3) design_example # # Exon2 Exon3 Samp2 Samp3 Samp4 Exon2: Cond1 Exon3: Cond1 # # 1 0 0 0 0 0 0 0 # # 2 1 0 0 0 0 0 0 # # 3 0 1 0 0 0 0 0 # # 4 0 0 1 0 0 0 0 # # 5 1 0 1 0 0 0 0 # # 6 0 1 1 0 0 0 0 # # 7 0 0 0 1 0 0 0 # # 8 1 0 0 1 0 1 0 # # 9 0 1 0 1 0 0 1 # # 10 0 0 0 0 1 0 0 # # 11 1 0 0 0 1 1 0 # # 12 0 1 0 0 1 0 1