3．1数据类型和输入格式

理论上，OVESEG-test可以应用于任何分子表达数据类型，只要这些数据类型经过转换后可以被正态分布建模，例如微阵列和蛋白质组数据的log2-transformation, RNAseq计数的logCPM, DNA甲基化β值的logit2-transformation。如果变换后存在均值-方差关系，limma: vooma() /轰()需要首先执行，得到的权重矩阵作为的输入OVESEG．

对输入表达式数据的要求:

• (log2/logCPM/logit2)。
• 已经过归一化和批效应去除处理。
• 没有缺失值;含有缺失值的分子应在oveseg测试前去除。
• 无低表达分子;否则会影响测试结果。

输入表达式数据应该存储在一个矩阵中。数据帧或总结实验对象也被接受，并将在分析之前在内部强制转换为矩阵格式。行对应探测，列对应样本。组织/细胞标签必须与表达矩阵中的每一列相匹配。如果提供了权重矩阵，则必须匹配表达式矩阵中的每个点。

3.2微阵列数据

我们使用纯化的B细胞，CD4+ T细胞和CD8+ T细胞的数据集GSE28490)为例，对微阵列数据进行oveseg测试。

library(OVESEG) #导入数据data(RocheBT) y <- RocheBT$y #5000*15 matrix group <- RocheBT$group #15 labels #filter low- expression probes groupMean <- sapply(levels(group)， function(x) rowMeans(y[，group == x])) groupMeanMax <- apply(groupMean, 1, max) keep <- (groupMeanMax > log2(64)) y <- y[keep，] # OVESEGtest rtest1 <- OVESEGtest(y, group, NumPerm = 999，BPPARAM=BiocParallel::SnowParam()) #>计算零假设的后后概率#>排列前2个组#>排列前3个组#>计算p值#>计算每个组的p值

注意，有许多低表达的探测过滤方法。在这里，我们过滤了即使在最高表达的组中，其平均表达量也小于阈值的探测。如果均值-方差关系明显，可以考虑应用limma: vooma ()首先。

#OVESEG-test yvooma <- limma::vooma(y, model.matrix(~0+factor(group))) rtest2 <- OVESEGtest(yvooma$E, group, weights=yvooma$weights, NumPerm = 999, BPPARAM=BiocParallel::SnowParam()) #>计算零假设的后验概率#>排列前2个组#>排列前3个组#>计算p值#>具体计算每个组的p值

3．3RNAseq计数数据

我们使用纯化的B细胞，CD4+ T细胞和CD8+ T细胞的数据集GSE60424)为例，对RNAseq计数数据进行OVESEG-test。

library(OVESEG) #导入数据data(countBT) count <- countBT$count #10000*60矩阵组<- countBT$group #60标签#过滤低表达探针groupMean <- sapply(levels(group)， function(x) rowMeans(count[，group == x])) groupMeanMax <- apply(groupMean, 1, max) keep <- (groupMeanMax > 2) count <- count[keep，] #OVESEG-test lib。size <- max(colsum (count)) yvoom <- limma::voom(count, model.matrix(~0+factor(group)))， lib。size = library .size) rtest3 <- OVESEGtest(yvoom$E, group, weights=yvoom$weights, NumPerm = 999, BPPARAM=BiocParallel::SnowParam()) #>计算零假设的后验概率#>排列前2个组#>排列前3个组#>计算p值#>具体计算每个组的p值

注意，有许多低表达的探测过滤方法。在这里，我们过滤了即使在最高表达组中，其平均计数也小于阈值的探针。lib.size和其他参数limma:轰()可以手动设置根据limma包中。

3.4假定值

OVESEGtest返回三个p值:pv.overall，pv.oneside而且pv.multiside．pv.overall由所有排列计算而不考虑上调的子类型。pv.oneside子类型特定的p值是专门为每个探针的上调子类型计算的。pv.multiside是pv.oneside乘以K (K比较校正)并截断为1。更多细节可以在论文中找到(Chen et al. 2019)．

4．1测试统计数据

oveseg检验统计量定义为\ [\ mathbf {\ max_ {k = 1,……，K}\left\{min_{l \neq k}\left\{ \frac{\mu_k(j)-\mu_l(j)}{\sigma(j)\sqrt{\frac{1}{N_k}+\frac{1}{N_l}}} \right\}\right\}}\]在哪里\ (\ mathbf {\ mu_k (j)} \)为k亚型中基因j对数表达量的均值，而进行p值估计的排列时间较长，OVESEGtstat ()如果用户只需要测试统计数据来对基因进行排名，则此方法很有用:

# ovesegg -test statistics rtstat1 <- OVESEGtstat(y, RocheBT$group) rtstat2 <- OVESEGtstat(yvooma$E, RocheBT$group, weights=yvooma$weights) rtstat3 <- OVESEGtstat(yvoom$E, countBT$group, weights=yvoom$weights)

4.2零假设分量的后验概率

oveseg检验的零假设被建模为多个组件的混合，其中组件的权重由所有探针的后验概率估计。这些后验概率也由OVESEGtest ()．如果用户只想观察探测后验概率，postProbNull ()是有益的。

ppnull1 <- postProbNull(y, RocheBT$group) ppnull2 <- postProbNull(yvooma$E, RocheBT$group, weights=yvooma$weights) ppnull3 <- postProbNull(yvom $E, countBT$group, weights= yvom $weights)

利用函数进行后验概率的累加和归一化nullDistri ()，在零假设条件下，我们也可以得到一个亚型被上调的概率。概率越高的亚型，假阳性mg值越高。

nullDistri(ppnull1) #> B CD4 CD8 #> 0.2544411 0.3699901 0.3755688 nullDistri(ppnull2) #> B CD4 CD8 #> 0.25551930 0.3689916 0.3758154 nullDistri(ppnull3) # bbb CD4 CD8 #> 0.2555459 0.3664029 0.3780513

完全有\ (\ mathbf {2 ^ k - 1} \)在真实数据集中的表达式模式类型:仅用1、2、…或K个子类型表示的探测。未在任何子类型中表示的探针应该在预处理期间被过滤。利用函数积累和归一化零假设和备择假设的后验概率patternDistri ()可以向我们展示所有亚型之间可能的上调模式的比率:

pattern <- patternDistri(ppnull1) #或pattern <- patternDistri(rtest1)

图案的比例可以说明如下。

库(gridExtra)库(网格)图书馆(reshape2)图书馆(ggplot2) gridpatterns < -函数(ppnull) {K < -长度(ppnull标签美元)cellcomp < - patternDistri (ppnull) cellcomp < - cellcomp[秩序(cellcomp [K + 1],减少= TRUE),] #条形图显示模式比例DF1 < - data.frame(排名= seq_len (2 ^ K - 1), cellcomp) p1 < - ggplot (DF1, aes (x =排名,y = Wpattern)) + geom_bar(统计=“身份”)+ scale_y_continuous(标签=尺度::百分比)+ theme_bw (base_size = 16) + (axis.text主题。X = element_blank()， axis.title.x = element_blank()， plot.margin=unit(c(1,1，-0.2,1)， "cm")， panel.grid.minor = element_line(size = 1)， panel.grid.major = element_line(size = 1)， panel.grid.major = element_line(size = 1)，边境= element_blank (), axis.ticks.x = element_blank ()) + ylab(“的百分比在某些亚型(s)”)#模式作为轴DF2 < - data.frame(排名= seq_len (2 ^ K - 1), cellcomp[,——(K + 1)]) DF2 < -融化(DF2 id.var =“排名”)p2 < - ggplot (DF2, aes (x =, y =价值,填补=变量))+ geom_bar(统计=“身份”)+ scale_fill_brewer(面板=“关于我校”)+主题(线= element_blank (), axis.text.x = element_blank (), axis.text.y = element_blank()、title = element_blank(),面板。Background = element_rect(fill = "white")，图例。position="bottom"， legend.key.size =unit(1.2，"line")， plot.margin=unit(c(-0.2,1,1,1)， "cm")) + scale_y_reverse() + guides(fill = guide_legend(nrow = 1, byrow = TRUE)) g1 <- ggplotGrob(p1) g2 <- ggplotGrob(p2) colnames(g1) <- paste0(seq_len(ncol(g1))) colnames(g2)) g <- gtable_combine(g1, g2, along=2) panels <- g$layout$t[grepl("panel"，G $layout$name)] n1 <- length(G $heights[panels[1]]) n2 <- length(G $heights[panels[2]]) #为面板分配新的(相对)高度#注意这里的*4以获得不同的高度G $heights[panels[1]] <- unit(n1*4，“null”)G $heights[panels[2]] <- unit(n2，“null”)return(G)} grid.newpage() grid.draw(gridpatterns(ppnull1)) #或grid.draw(gridpatterns(rtest1))