用黑羊法分析离群值
04/28/2022
摘要
Blacksheep在两个种群之间进行差异比较的情况下进行极值分析。其基本机制是检验两个总体之间的异常值比例，并评估异常值比例之间的统计差异。包
blacksheepr 1.11.0
内容
0．1简介
离群值分析是一种成熟的方法，用于在其余人群的背景下探索极端值。在生物学背景下，探索这些异常值比例的变化可以阐明亚种群之间的差异，表明可以进一步探索的潜在差异特征。Blacksheep是一个致力于将这种分析改进为生物数据环境中异常值分析的功能工具的项目。这些数据可以有多种形式:蛋白质、磷蛋白、RNA、CNA等。黑羊的输入是某种形式的计数数据，以及用于比较的指示种群的注释数据。主要的函数调用是提婆(微分极值分析)。本小插图将通过一些标准用例，说明blacksheep的功能，并希望回答问题，并显示其在生物探索中的有用性。
0.2安装
安装类似于所有的Bioconductor包。启动R并运行以下代码行进行安装:
如果(!requireNamespace("BiocManager"， quiet = TRUE)) install.packages("BiocManager") BiocManager::install("blacksheep ")方法加载库图书馆调用
库(blacksheepr)0.3输入数据
0.3.1计数数据
计数数据应该是一个计数表，沿x轴是样本，沿y轴是特征，标签是行名和冒号。
数据(“sample_phosphodata”)sample_phosphodata (1:5, 1:5)NP_116026-S179s NA 0.97943811 0.4830637 NA NA NA NA NA NA NA NA NA ## tcga - a2 - a154 ## NP_149131-S462s T469t -1.512795 ## np_1182345 - s192s NA ## np_1182345 - s413s 2.658696 ## NP_116026-S179s NA ## NP_002388-S192s NA ##0.3.2注释数据
Annotation数据应该是一个与计数数据中包含的相同样本的表，沿y轴列出。然后，每一列都是要执行的比较，列的值是某种形式的二进制系统，表示属于每个子组的样本。注意-这些应该是字符串。如果您的列实际上包含有用的信息，如“高”/“低”，“突变体”/“WT”等，而不是“1/0”或任何其他非描述性二进制，则信息更丰富。
数据(“sample_annotationdata”)sample_annotationdata [1:5]## TCGA-AO-A12D Her2 not_基础not_LumA not_LumB阴性阴性## TCGA-BH-A18Q not_Her2 not_基础not_LumA not_LumB阴性阴性## TCGA-C8-A130 Her2 not_基础not_LumA not_LumB阳性阳性## TCGA-C8-A138 Her2 not_基础not_LumA not_LumB阳性阴性## TCGA-E2-A154 not_Her2 not_基础LumA not_LumB阳性阳性0.3.3SummarizedExperiment数据
blacksheeper将summarizeexperiment对象带入主程序提婆函数调用(下文将详细介绍)。因此blacksheepr可以从任何summarizeexperiment对象开始，只要summarizeexperiment的colData包含与上面描述的数据中的列名相同的行名，并且元数据由BINARY分类组成。注意，这可能需要一些额外的格式。附录中有关于辅助函数的说明，这些函数可以帮助进行这种格式化。
1磷酸蛋白
在下一节中，我们将通过一个使用磷蛋白数据的异常值分析的例子。输入的数据是从Github最初使用的是来自TCGA和CPTAC的乳腺癌数据。doi: 10.1038 / nature18003
1．1输入数据
1.1.1阅读注释
阅读注释表。
数据(sample_annotationdata) <- sample_annotationdata comptable[1:5，]## TCGA-AO-A12D Her2 not_基础not_LumA not_LumB阴性阴性## TCGA-BH-A18Q not_Her2 not_基础not_LumA not_LumB阴性阴性## TCGA-C8-A130 Her2 not_基础not_LumA not_LumB阳性阳性## TCGA-C8-A138 Her2 not_基础not_LumA not_LumB阳性阴性## TCGA-E2-A154 not_Her2 not_基础LumA not_LumB阳性阳性暗(comptable)## [1] 76 61.1.2读入磷数据
接下来，我们需要有我们正在分析的实际计数数据。请注意，这些数据是由特征-蛋白质-组成的，然后有许多亚特征-磷酸化位点。这方面的数据将直接与后面的分析相关，因为blacksheep将对该数据进行聚合，并将数据折叠到主要特征上。稍后再详细介绍。还要注意，这些数据已经被预规范化了。对于更多关于数据输入的处理提婆，详见附录。
Data (sample_phosphodata) phosphotable <- sample_phosphodata phosphotable[1:5,1:5]NP_116026-S179s NA 0.97943811 0.4830637 NA NA NA NA NA NA NA NA NA ## tcga - a2 - a154 ## NP_149131-S462s T469t -1.512795 ## np_1182345 - s192s NA ## np_1182345 - s413s 2.658696 ## NP_116026-S179s NA ## NP_002388-S192s NA ##暗(phosphotable)## [1] 15532 761.1.3从我们的数据创建一个summary experiment。
blacksheep -作为Bioconductor宇宙的一部分-从包含感兴趣的分析和底层元数据的单个对象开始。一个SummarizedExperiment对象易于创建-并有助于确保与每个示例关联的数据和元数据之间没有错位。注:如下所示，计数数据必须格式化为MATRIX，注释数据必须格式化为DATAFRAME。
suppressPackageStartupMessages(library(summarizeexperiment)) blacksheep_SE <- summarizeexperiment (assays=list(counts=as.matrix(phosphotable))， colData=DataFrame(comptable)) blacksheep_SE##类:summarize实验## dim: 15532 76 ##元数据(0):## assays(1):计数## rownames(15532): NP_149131-S462s T469t NP_061893-S25s…## NP_542417-S1276s NP_001136254-S832s ## rowData names(0): ## colnames(76): TCGA-AO-A12D tcga - bhh - a18q…TCGA-BH-A0C7 TCGA-A2-A0SX ## colData名称(6):PAM50_Her2 PAM50_Basal…ER_Status PR_Status1．2运行deva(差分极值分析)
的提婆函数有许多步骤，下面分别描述这些步骤。但是请注意，单个步骤只需要用于特定的查询或更改。在一般情况下，提婆函数本身对于所需的分析应该是足够的。
deva_out <- deva(se = blacksheep_SE, analyze_negative_outliers = FALSE, aggregate_features = TRUE, feature_delineator = "-"， fraction_samples_cutoff = 0.3, fdrcutoffvalue = 0.1)1.2.1 "天神参数
- se -作为deva输入的summary experiment
- analyze_negative_outliers——分析是查看正异常值还是负异常值。
- aggregate_features——分析是否将特征分解为一个主要术语。适用于蛋白质上的磷酸基、蛋白质同型等。
- 类时，描述主要特征的字符标记
参数 - Fraction_samples_cutoff -一组中有多少样本必须包含一个异常值，才能认为该特征是显著的。该参数用于避免单个样本的结果过于偏倚。(见附录)
- fdrcutoffvalue -用于重要输出的FDR值。
1.3探索deva输出
blacksheep的输出是带有期望结果的列表。有一个outlier_analysis和一个significant_heatmaps进行所需的分析。在本例中，我们只对正离群值进行分析。因此输出将是2项长-热图和表的积极异常值分析。
名(deva_out)## [1] "outlier_analysis_out" "significant_heatmaps" "fraction_table" ## [4] "median_value_table" "outlier_boundary_table"的outlier_analysis是一个嵌套列表的分析-与一个分析每个比较列
名称(deva_out pos_outlier_analysis美元)# #空的significant_heatmapsSection是热图遵从的嵌套列表，同样每个比较列有一个分析
名称(deva_out significant_pos_heatmaps美元)# #空1.3.1deva_results()函数
deva_results ()是一个帮助探索数据的实用函数。如果你使用deva_results在deva_out对象，它将返回返回重要结果的已执行分析的列表。
deva_results (deva_out)## [1] "outlier_analysis_out" "significant_heatmaps" "fraction_table" ## [4] "median_value_table" "outlier_boundary_table"deva_out的其他参数是ID而且类型．ID是一个关键字来抓取所需的分析。类型是以下选项之一:表格或的热图指定您想要获取的分析对象，后面跟着特定分析的ID。fraction_table是输出表，它显示了每个特征的每个样本的离群值的比例(注意，如果没有进行聚合，这将与离群值表相同)。中位数而且边界返回中值列表，以及每个特征的离群值边界值。
subanalysis_Her2 <- deva_results(deva_out, ID = "Her2"， type = "table") head(subanalysis_Her2)# # 1 # #基因pval_more_pos_outliers_in_PAM50_Her2__Her2 NP_000215 0.000876167768317827 # # 2 NP_000393 0.0104700645400028 # # 3 NP_001002814 0.312790178663637 # # 4 NP_001003408 0.540506340033899 # # 5 NP_001005751 # # 6 NP_001020262 0.0519534804138539 # # pval_more_pos_outliers_in_PAM50_Her2__not_Her2 # # 1 # # 2 # 3 # 4 # 5 # 1 # 6 # # # fdr_more_pos_outliers_in_PAM50_Her2__Her2 # 0.0055210906126839 # 1 # 2 # 0.0260165240084919 # 0.371721661600264 # 3 # 4 # 0.575604154321814 # # 0.0774579162533822 # # 5 # # 6fdr_more_pos_outliers_in_PAM50_Her2__not_Her2 PAM50_Her2__Her2_nonposoutlier ## 1 32 ## 2 12 ## 3 47 ## 4 65 ## 5 1 100 ## 6 31 ## PAM50_Her2__Her2_posoutlier PAM50_Her2__not_Her2_nonposoutlier ## 1 6 186 ## 2 4 67 ## 3 4 285 ## 4 4 392 ## 5 5 562 ## 6 5 165 ## PAM50_Her2__not_Her2_posoutlier ## 1 3 ## 2 2 ## 3 14 ## 4 18 ## 5 33 ## 6 81.3.2结果
对于您输入的比较表的每一列(占用colData在summarizeexperimental对象中)，提婆如果有任何适用的结果，将输出一个分析表。表的一个例子如下:
subanalysis_Her2 <- deva_results(deva_out, ID = "Her2"， type = "table") head(subanalysis_Her2)# # 1 # #基因pval_more_pos_outliers_in_PAM50_Her2__Her2 NP_000215 0.000876167768317827 # # 2 NP_000393 0.0104700645400028 # # 3 NP_001002814 0.312790178663637 # # 4 NP_001003408 0.540506340033899 # # 5 NP_001005751 # # 6 NP_001020262 0.0519534804138539 # # pval_more_pos_outliers_in_PAM50_Her2__not_Her2 # # 1 # # 2 # 3 # 4 # 5 # 1 # 6 # # # fdr_more_pos_outliers_in_PAM50_Her2__Her2 # 0.0055210906126839 # 1 # 2 # 0.0260165240084919 # 0.371721661600264 # 3 # 4 # 0.575604154321814 # # 0.0774579162533822 # # 5 # # 6fdr_more_pos_outliers_in_PAM50_Her2__not_Her2 PAM50_Her2__Her2_nonposoutlier ## 1 32 ## 2 12 ## 3 47 ## 4 65 ## 5 1 100 ## 6 31 ## PAM50_Her2__Her2_posoutlier PAM50_Her2__not_Her2_nonposoutlier ## 1 6 186 ## 2 4 67 ## 3 4 285 ## 4 4 392 ## 5 5 562 ## 6 5 165 ## PAM50_Her2__not_Her2_posoutlier ## 1 3 ## 2 2 ## 3 14 ## 4 18 ## 5 33 ## 6 8输出顺序为:* col1 -基因列表* col2-3 -基因的p值显著不同，它出现的列表示该基因在group1或group2中较高。* col4-5 - pvalue的FDR值* col6-9 -统计检验背后的原始数据，显示每个组有多少异常值和非异常值。
热图作为重要基因的快照，满足参数化fdr截止在提婆函数调用。可以以与分析表类似的方式访问它们。不过要注意的是，只有存在显著基因时才会有热图。如果没有满足fdr截断的基因，那么就不会产生热图。
subanalysis_Her2_HM <- deva_results(deva_out, ID = "Her2"， type = "heatmap") subanalysis_Her2_HM
热图对象本身就是ComplexHeatmap对象——并且在调用时将被绘制，可以在默认的Rplot中，也可以保存到外部文件中。
##不运行##单独输出到pdf pdf("outfile.pdf") draw(subanalysis_Her2_HM) dev.off()2分段分析
本节假设已经完成了前面读入数据和注释的步骤下一节将演示deva的各个步骤。注意用户可以从来没有需要调用特定的步骤，但如果需要进行特定的调整，或者需要提取中间步骤，则可以遵循此工作流以查看如何生成分析。
2.1创建分组
根据元数据创建子组。注意comparison_groupings函数通过遍历比较列来创建组，并为所有比较创建第一个子组，然后为所有比较创建第二个子组。顺序取决于向下移动时遇到的第一个子类别。当您稍后查看比较时，这个顺序很重要，但是这个信息将包含在输出表中以避免混淆，稍后将对此进行详细介绍。
<- comparison_groupings(comptable) ##打印出我们的前5个组中的前6个样本lapply(groupings, head)[1:5]# # # # $ PAM50_Her2__Her2[1]“TCGA-AO-A12D”“TCGA-C8-A130”“TCGA-C8-A138”“TCGA-C8-A12L”“TCGA-C8-A12T”# #[6]“TCGA-A2-A0EQ”PAM50_Basal__not_Basal美元# # # # # #[1]“TCGA-AO-A12D”“TCGA-C8-A130”“TCGA-C8-A138”“TCGA-E2-A154”“TCGA-C8-A12L”# #[6]“TCGA-A2-A0EX”PAM50_LumA__not_LumA美元# # # # # #[1]“TCGA-AO-A12D”“TCGA-BH-A18Q”“TCGA-C8-A130”“TCGA-C8-A138”“TCGA-C8-A12L”# #[6]“TCGA-BH-A0AV”PAM50_LumB__not_LumB美元# # # # # #[1]“TCGA-AO-A12D”“TCGA-BH-A18Q”“TCGA-C8-A130”“TCGA-C8-A138”“TCGA-E2-A154”## [6] "TCGA-C8-A12L" ## ## $ER_Status__Negative ## [1] "TCGA-AO-A12D" "TCGA-BH-A18Q" "TCGA-C8-A12L" "TCGA-BH-A0AV" "TCGA-A2-A0CM" ## [6] "TCGA-A2-A0EQ"2．2制作离群值表
下一个函数make_outlier_table将接收countdata并输出一个已转换为显示异常值的表。值为0表示它不是异常值，1表示它是正方向上的异常值，如果参数设置为分析负异常值，则-1表示负方向上的异常值。
这个函数的输出是一个对象列表，它取决于输入和指定的参数。也就是说，只有当analyze_negative_outliers的参数被关闭时，它才会输出$upperboundtab，如果analyze_negative_outliers参数被打开，它才会输出$lowerboundtab
##执行函数reftable_function_out <- make_outlier_table(phosphotable, analyze_negative_outliers = FALSE) ##查看输出对象的名称(reftable_function_out)## [1] "outliertab" "sampmedtab" "upperboundtab"##将它们分配给单个变量outliertab <- reftable_function_out$outliertab upperboundtab <- reftable_function_out$upperboundtab sampmedtab <- reftable_function_out$sampmedtab ##注意，我们将只使用异常值表-现在看起来像这样outliertab[1:5,1:5]TCGA-AO-A12D TCGA-BH-A18Q TCGA-C8-A130 tcga - c8 - s49s T469t 0 NA NA 00 NA ## NP_001182345-S413s 00 0 NP_116026-S179s NA 00 0 NA ## np_001182345 - s192s NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ## tcga - b2 - a154 ## np_149131 - s42s T469t 0 ## NP_116026-S179s NA ## NP_002388-S192s NA2．3汇总异常值
对于我们的每个组，运行离群值表并计算离群值和非离群值的总数。对于磷(以及本例)，我们将使用AGGREGATION函数以每个蛋白质为基础将我们的计数聚合在一起
这个函数的输出是一个对象列表，它取决于输入和指定的参数。对于磷酸化蛋白-你可以有数据，每个蛋白质有几个磷酸化位点。作为分析的一部分，一种选择是聚合该蛋白质的数据，并将异常值信息分解为单个特征。把aggregate_features = TRUE将执行此功能，而feature_delineator在ex) Feature1-1和Feature1-2-1上使用“-”描述符的字符串是否会折叠到Feature1上
的输出aggregate_features = TRUE将包含两个额外的对象。它将输出正常的异常值表和边界表，以及“aggoutlierstab”和“fractiontab”。“aggoutlierstab”将是折叠表，汇总每个特征的异常值数量。“fractiontab”将返回每个样本每个特征的异常值百分比，给定可用的信息——这对以后的进一步过滤很重要。ex) 1,0,NA >> 1/2
count_outliers_out <- count_outliers(groupings, outliertab, aggregate_features = TRUE, feature_delineator = "-") grouptablist <- count_outlierout $grouptablist aggoutliertab <- count_outlierout $aggoutliertab fractiontab <- count_outliers_out$fractiontab names(grouptablist)##[1]“PAM50_Her2__Her2”“PAM50_Basal__not_Basal”“PAM50_LumA__not_LumA”##[4]“PAM50_LumB__not_LumB”“ER_Status__Negative”“PR_Status__Negative”##[7]“PAM50_Her2__not_Her2”“PAM50_Basal__Basal”“PAM50_LumA__LumA”##[10]“PAM50_LumB__LumB”“ER_Status__Positive”“PR_Status__Positive”每个表都有特征计数，并存储在样本中，进入计数
名称(grouptablist PAM50_Her2__Her2美元)## [1] "feature_counts" "samples"特性计数的例子如下:
头(grouptablist PAM50_Her2__Her2 feature_counts美元)## np_000011 9 1 ## np_000012 00 ## np_000019 1 0 ## np_000034 4 0 ## np_000035 7 1 ## np_000042 00进行分析的样本是什么:
grouptablist PAM50_Her2__Her2美元样本[1] " tcga-ao-a12d " " tcga-c8-a130 " " tcga-c8-a138 " " tcga-c8-a12l " " tcga-c8-a12t " ## [6] " tcga-a2-a0eq " " tcga-ar-a0tx " " tcga-a8-a076 " " tcga-c8-a135 " " tcga-c8-a12z " ## [11] " tcga-ao-a0je " " tcga - a8 - a12g "2.4运行异常值分析
使用表格，运行异常值分析，以寻找组间异常值的丰富。注意，这个函数有一个内置的功能，可以将表写入外部文件，我们现在不使用这个参数，但是可以通过转换to参数来设置write_out_tables = TRUE
在这个异常值分析中，我们将包含一个额外的过滤器。的fractiontab上一步中的输出有一个度量，度量每个组中有一个离群值的样本数量。在下一步中，我们将使用该信息来过滤只有至少的特征0.3(30%)有异常值的样本。这个过滤器很重要，因为我们的聚合步骤分解了每个样本中的所有站点，然后我们计算了总数。如果一个样本的所有站点都是异常值——虽然有趣——但这并不意味着我们整个群体都有过高的代表性——只是一个样本。该过滤器仅通过挑选具有多个具有异常值的样本的特征来实现更清晰的分析。注意，可以通过只留下fraction_table参数out或作为零
Outlier_analysis_out <- outlier_analysis(grouptablist = grouptablist, fraction_table = fractiontab, fraction_samples_cutoff = 0.3) names(Outlier_analysis_out)[1]“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”“[6]”头(outlier_analysis_out outlieranalysis_for_PAM50_Her2__Her2_vs_PAM50_Her2__not_Her2美元)# # 1 # #基因pval_more_pos_outliers_in_PAM50_Her2__Her2 NP_000215 0.000876167768317827 # # 2 NP_000393 0.0104700645400028 # # 3 NP_001002814 0.312790178663637 # # 4 NP_001003408 0.540506340033899 # # 5 NP_001005751 # # 6 NP_001020262 0.0519534804138539 # # pval_more_pos_outliers_in_PAM50_Her2__not_Her2 # # 1 # # 2 # 3 # 4 # 5 # 1 # 6 # # # fdr_more_pos_outliers_in_PAM50_Her2__Her2 # 0.0055210906126839 # 1 # 2 # 0.0260165240084919 # 0.371721661600264 # 3 # 4 # 0.575604154321814 # # 0.0774579162533822 # # 5 # # 6fdr_more_pos_outliers_in_PAM50_Her2__not_Her2 PAM50_Her2__Her2_nonposoutlier ## 1 32 ## 2 12 ## 3 47 ## 4 65 ## 5 1 100 ## 6 31 ## PAM50_Her2__Her2_posoutlier PAM50_Her2__not_Her2_nonposoutlier ## 1 6 186 ## 2 4 67 ## 3 4 285 ## 4 4 392 ## 5 5 562 ## 6 5 165 ## PAM50_Her2__not_Her2_posoutlier ## 1 3 ## 2 2 ## 3 14 ## 4 18 ## 5 33 ## 6 82.5使用热图生成函数绘制结果
在获得结果之后，在图形中获得结果的快照是很有用的。blacksheeper包含一个实用函数，用于输出带有自定义注释和数据的热图。使用原始注释数据的绘图函数，并使用您想要表示的任何信息填充热图。在这种情况下，我们将用每个特征的异常值的分数填充表。outlier_analysis对象用于选择差异基因。注意，您可以使用给定的参数直接在函数中写出绘图write_out_plot或者从函数中保存的对象是热图，因此您可以打开自己的PDF并使用下面的注释代码打印出来。注意，元表必须在所需的ORDER中，并且这个函数不会为您排序。
可绘图<- comptable[做。调用(order, c(decreasing = TRUE, data.frame(comptable[,1:ncol(comptable)]))),] hm1 <- outlier_heatmap(outlier_analysis_out = outlier_analysis_out, counttab = fractiontab, metatable = plottable, fdrcutoffvalue = 0.1) ## To output heatmap to pdf outside of the function #pdf(paste0(outfilepath, "test_hm1.pdf")) #hm1 #junk<-dev.off()hm1 print_outlieranalysis_for_PAM50_Her2__Her2_vs_PAM50_Her2__not_Her2美元
图1:输出Heatmap
3.附录
3．1格式化注释表
格式化注释表是确保blacksheepr平稳运行的关键步骤。关键在于注释的行名与分析数据的列名完全匹配，并且注释列是BINARY，因此每列只有2个子类别。
3.1.1使用make_comparison_columns实用函数
有一个内置的效用函数make_comparison_columns帮助将一个多因素列转换为几个二进制列。用户可以输入任意多列，该函数将输出一个表，其中每个多阶乘列都转换为若干二进制列
Dummyannotations <- data.frame(comp1 = c(1,1,2,2,3,3)， comp2 = c("red"， "blue"， "red"， "blue"， "blue"， "green"， "green")， row.names = paste0("sample"， seq_len(6))## comp1 comp2 ## sample1 1红## sample2 1蓝## sample3 2红## sample4 2蓝## sample5 3绿## sample6 3绿使用make_comparison_columns函数创建二进制列##抽样抽样1“1”“not_2”“not_3”“not_blue”“not_green”##抽样抽样2“1”“not_2”“not_3”“not_red”“blue”“not_green”##抽样抽样3“not_1”“2”“not_3”“not_red”“not_blue”“not_green”##抽样抽样4“not_1”“2”“not_3”“not_red”“blue”“not_green”##抽样抽样5“not_1”“not_2”“3”“not_red”“not_blue”“绿色”##抽样抽样6“not_1”“not_2”“3”“not_red”“not_blue”“绿色”3.2运行blacksheep函数为其他组学数据
Blacksheepr是一套具有算法的通用工具，可用于前面描述的分析和注释格式的任何类型的组学数据。沿着这些路线，对于其他组学数据类型有一些常见的实践。
3.2.1之上使用RNAseq数据运行deva
原理是一样的，唯一的主要区别是你不会聚合到一个主要的特性上，所以aggregate_features参数应保留为默认值假．还请注意，如果希望输出所有基因，而不考虑其重要性，则可以通过设置FDRcutoff到1，和fraction_value_cutoff为0。
3．3处理使用deva运行的数据
提婆它的核心是基于与中值的差异来探索异常值。严重倾斜的数据需要确保准确的分析。Blacksheepr包括一个归一化函数，它将执行一个中位数的比率归一化和一个log2变换。方法的示例请参见下面pasilla数据集:
库(pasilla)##加载所需的包:DEXSeq##加载所需的包:BiocParallel##加载所需的包:DESeq2##加载所需的包:AnnotationDbi##加载所需的包:RColorBrewerpasCts <- system。文件(“extdata”、“pasilla_gene_counts。tsv"， package="pasilla") cts <- as.matrix(read.csv(pasts,sep="\t"，row.names="gene_id")) norm_cts <- deva_normalization(cts, method =" MoR-log")