fgga 1.5.0
随着基因组数据量的增长,计算方法对于提供对基因注释的初步了解变得至关重要。当注释结果的可解释性是一个主要问题时,基于层次集成分类技术的自动基因本体(GO)注释方法特别有趣。在这些方法中,通过检查预定义GO关系的一致性,利用基本二元分类器计算的原始GO术语预测。FGGA是一种因子图方法,用于自动GO注释问题的分层集成公式。在这种形式化方法中,首先基于GO结构构建核心因子图,然后对其进行充实,以考虑GO项预测的噪声性质。因此,从原始GO-term预测开始,使用因子图节点之间的迭代消息传递算法来计算目标GO-term的边际概率。
该方法在原文中有详细说明出版物[1,2],但这个简短装饰图案说明如何使用FGGA从给定生物体的基因产物序列中预测一组GO-term (GO节点标签)。目的是提高现有电子标注的质量,为Blast等传统方法无法分类的未知序列提供新的标注方法。
因此,FGGA是一种自动标注蛋白质序列的工具,然而,其他类型的醒目功能基因产物的标注也是可能的,例如长链非编码rna。FGGA通过对长链非编码RNA序列的一级、二级和三级结构进行表征,提高碱基二元分类器的置信度,从而为改进长链非编码RNA序列的标注提供了机会。沿着小插图,我们使用的蛋白质编码基因的数据从canis familiaris生物体获得与UniProt。
的fggaR源码包可以从Bioconductor库或GitHub库。这个R包包含一个实验数据集作为示例,一个预运行的R对象和运行FGGA所需的所有函数。
##从Bioconductor存储库如果(!)requireNamespace("BiocManager",悄悄地= TRUE) {install.packages("BiocManager")} BiocManager::install("fgga") ##或使用devtools从GitHub存储库BiocManager::install("devtools") devtools::install_github("fspetale/fgga")
安装包后,可以通过以下命令将其加载到R工作区中
库(fgga)
目前,该方法直接支持基因产物序列表征的数据。这种特征可以根据专家的标准生成。例如,可能的特征可以是:主题包含了,物理化学性质,信号肽等等。允许的值可以是数字、布尔或字符类型。然而,为了提高二元分类算法的效率,建议使用规范化的数值。这些数据集必须每个go术语至少有50个注释才能训练FGGA模型。
在本节中,我们应用FGGA来预测的生物学功能,go术语犬狼疮蛋白质。我们下载6962犬狼疮, Cf,蛋白质序列及其GO注释来自UniProt然后通过R包进行理化性质表征肽。
让我们在工作区中导入玩具数据集:
#加载Canis狼疮数据集和示例R对象数据(CfData)
#看总结实验对象总结(CfData) # >长度类模式# > dxCf 501264那个没有数字# > graphCfGO 1296那个没有数字# > indexGO 2那个没有-列表# > tableCfGO 250632那个没有数字# > nodesGO 36那个没有字# > varianceGOs 36那个没有——数字#数据的特征信息(CfData dxCf美元)# > [1]6962 72 colnames (CfData dxCf美元)[1:20]# >[1]“Numb_Amino”“A”“C”“D”“E”“F”“G”# >[8]“H”“我”“K”“L”“M”“N”“P”# >[15]“Q”“R”“S”“T”“V”“W”rownames(CfData$dxCf)[1:10] #> [1] "E2QTB2" "F1PQ93" "E2RQU4" "E2RN67" "E2RGI4" "E2RFJ5" "J9P0Z9" "E2RL61" "F1PAF1" "J9P1G2" head。矩阵(CfData$dxCf[, 51:61], n = 10) #> pkc。po_charged_res2 Pos_changed_res3 #> E2QTB2 -0.0790 0.0474 58225.55 4.9492 0.1175 0.1368 #> F1PQ93 -0.0953 0.0494 83933.13 4.6956 0.1250 0.1371 #> E2RQU4 -0.0117 0.1234 83933.57 9.0031 0.1344 0.1694 #> E2RN67 0.0185 0.0424 73821.11 7.9705 0.0848 0.1161 #> E2RGI4 0.0715 0.0368 28713.98 5.3223 0.1304 0.1423 #> E2RFJ5 -0.0194 0.0066 30101.85 4.6895 0.1049 0.1199 #> J9P0Z9 0.0010 0.0619 67510.79 6.3768 0.1267 0.1467 #> E2RL61 -0.0559 0.0226 163422.56 5.6757 0.1291 0.1463 #>F1PAF1 -0.0407 -0.0550 107808.45 5.4839 0.0932 - 0.1087 # > J9P1G2 -0.1119 0.0330 26295.98 9.1588 0.1040 - 0.1200 # > Neg_changed_res碳氢氮氧# > E2QTB2 0.1676 2518 3935 739 816 # > F1PQ93 0.1976 1209 1922 330 402 # > E2RQU4 0.1183 3779 5963 1033 1083 # > E2RN67 0.0818 3374 5180 910 914 # > E2RGI4 0.1542 1289 2035 331 387 # > E2RFJ5 0.1723 1327 2072 348 426 # > J9P0Z9 0.1333 3003 4694 808 908 # > E2RL61 0.1470 7269 11378 2012 2209 # > F1PAF1 0.1118 4755 7406 1294 1451 # > J9P1G2 0.08001166 1842 328 347 #查看GO data dim(CfData$tableCfGO) #> >962 36个colnames(CfData$tableCfGO)[1:10] #> [1] "GO:0000166" "GO:0003674" "GO:0003676" "GO:0003677" "GO:0003824" "GO:0004888" "GO:0005102" #> [8] "GO:0005215" "GO:0005488" "GO:0005515" rownames(CfData$tableCfGO)[1:10] #> [1] "E2QTB2" "F1PQ93" "E2RQU4" " e2rg4 " "E2RFJ5" "J9P0Z9" "E2RL61" "F1PAF1" "J9P1G2" head(CfData$tableCfGO)[,][1:8] #> go:0000166 go:0003674 go:0003676 go:0003677 go:0003824 go:0004888 go:0005102 go:0005215 #> e2qtb2 01 0000 0000 000 #> f1pq93 01 0000 0000 0 #> e2rqu4 1 10 0000 0000 0000 0000 #> e2rn67 01 0000 0000 01 #> e2rgi4 01 0000 001 #> e2rfj5 01 0000 0000 001 10 0000 00001 10 0000 0000 01 10 0000 0000 0000 0001 10 0000 000
tableCfGO一个二元矩阵是谁的\ ((i, j) _ {th} \)成分为1if蛋白质我\ \ ()注释为\ (GO-term_j \),否则为0。注意,dxCf和tablefgo的行名是相同的。我们的二元分类器每个go术语至少需要50个注释。因此,检查此条件,并消除不符合的go条款。
#检查GO-term的注释量apply(CfData$tableCfGO, MARGIN=2,sum) #> GO:0000166 GO:0003674 GO:0003676 GO:0003677 GO: 000582 GO:0005215 GO:0005488 #> 995 6962 1558 857 2679 516 520 565 5183 #> GO:0005515 GO:0005524 GO:0008144 GO:0016740 GO:0016787 GO:0030554 GO:0032553 #> 2738 682 764 1040 1097 878 870 707 882 #> GO:0032555 GO:0032559 GO:0035639 GO:0036094 GO:0043168 GO:0043168 GO:0043169 #> 871 702 847 1147 593 656 2273 1232 1307 #> GO:0046872 GO:0060089 GO:0097367 GO:0098772Go:0140096 Go:0140110 Go:1901265 Go:1901363 #> 1274 610 2523 992 650 978 615 995 2498
如果我们想要预测单个子域中的go项,BP, MF或CC,我们必须建立与这些go术语相关的GO-DAG。
库(GO.db)库(GOstats) mygraph <- GOGraph(CfData$nodesGO, GOMFPARENTS) #删除所有名为mygraph <- subGraph(CfData$nodesGO, mygraph)的根节点#我们将图调整为FGGA使用的格式mygraph <- t(as(mygraph, "matrix")) mygraphGO <- as(mygraph, "graphNEL") #我们搜索根GO项rootGO <- leaves(mygraphGO, "in") rootGO #> [1] "GO:0003674" "GO:0008144" plot(mygraphGO)
另一方面,如果您希望预测两个或三个子域,则应该使用preCoreFG函数。该函数建立了一个尊重推理GO约束的图,并将所选子域的GO项链接起来。
#我们在元胞成分子域myGOs <- c(CfData[['nodesGO']]], "GO:1902494", "GO:0032991", "GO:1990234", "GO:0005575")中添加GO-terms, #我们建立了一个尊重MF和CC子域推断的GO约束的图mygraphGO <- preCoreFG(myGOs, domains="MF-CC") plot(mygraphGO)
我们来看一个GO子图,GO项去:我映射到二值变量节点\ (x_i \)而go项之间的关系被映射到逻辑因子节点\ (f_k \)描述有效去:我真路径图约束下的配置。
模型fgga <- fgga2bipartite(mygraphGO)
现在,让我们使用MF子图用Cf数据的训练集构建二元SVM分类器模型。
#我们取Cf数据的一个子集来训练我们的模型idsTrain <- CfData$indexGO[["indexTrain"]][1:750] #我们建立我们的二元支持向量机分类器模型SVM <- lapply(CfData[["nodesGO"]]], FUN = svmTrain, tableGOs = CfData[["tableCfGO"]]][idsTrain,], dxfeatured = CfData[["dxCf"]][idsTrain,], graphGO = mygraphGO, kernelSVM = "radial")
#我们计算每个GO-term varianceGOs <- varianceGO(tableGOs = CfData[["tableCfGO"]][idsTrain,], dxcharacterated= CfData[["dxCf"]][idsTrain,], kFold = 5, graphGO = mygraphGO, rootNode = rootGO, kernelSVM = "radial")的可靠性
#> go:0000166 go:0003674 go:0003676 go:0003677 go:0004888 go:0005102 go: 000515 go:0005524 go: 000524 go:0008144 go:0016740 go:0016787 go:0017076 go:0030554 go:0032553 #> 3.153 0.911 1.039 2.026 2.428 1.128 2.171 0.891 0.960 #> go:0032555 go:0032559 go:0035639 go:0036094 go:0042802 go:0043168 go:0043169 #> 1.014 0.949 0.975 1.243 1.770 1.383 1.682 1.654 2.973 #> go:0046872Go:0060089 Go:0097159 Go:0097367 Go:0098772 Go:0140096 Go:0140110 Go:1901265 Go:1901363 #> 3.070 1.548 3.112 1.431 1.839 2.205 1.745 1.429 3.135
接下来,我们使用二元支持向量机分类器模型从测试集中预测每个go项。
dxtestcharacterization <- CfData[["dxCf"]][CfData$indexGO[["indexTest"]][1:50],] matrixGOTest <- svmGO(svmMoldel = modelSVMs, dxcharacteristic = dxtestcharacterization, rootNode = rootGO,varianceSVM = varianceGOs) head(matrixGOTest)[,1:8] #> GO:0000166 GO:0003674 GO:0003676 GO:0003677 GO: 0003877 GO: 0003882 GO:0005102 GO:0005215 #> J9P8X1 0.36261032 0.9999 0.5673337 0.29859880 0.23412172 0.53385507 0.5726663 #> F1PVD9 0.18383584 0.9999 0.6810226 0.70831700 0.1875907 0.18733458 0.50163000 0.3016141 #> F1PTQ9 0.12416758 0.9999 0.5674606 0.59928405 0.49460696 0.45442303 0.4494959 #> E2R8P1 0.62299290 0.9999 0.2439744 0.09448527 0.5969643 0.26052513 0.342110040.2724795 #> f6uwv4 0.03045405 0.9999 0.1485072 0.05638052 0.5342228 0.07593993 0.24838076 0.5000584 #> e2rqr5 0.83101925 0.9999 0.1871167 0.01917385 0.5083126 0.14971871 0.09746675 0.5850919
一旦因子图模型成品针对自动GO标注问题,定义了节点间的迭代消息传递算法成品可以用来计算变量节点的最大后验(MAP)估计\ (x_i \)实际建模去:我注释。
这个函数msgFGGA返回一个矩阵k行和米列分别对应Cf蛋白和MF go项。
matrixFGGATest <- t(apply(matrixGOTest, MARGIN = 1, FUN = msgFGGA, matrixFGGA = modelFGGA, graphGO = mygraphGO, tmax = 50,) head(matrixFGGATest)[,1:8] #> GO:0000166 GO:0003674 GO:0003676 GO: 0003877 GO:0004888 GO:0005215 #> J9P8X1 3.19485e -01 1.0000000 0.45560604 0.1359551814 0.6067942 0.107944377 0.45269695 0.5726663 #> F1PVD9 1.054230e-02 1.0000000 0.5511484659 0.2933632 0.012031746 0.45355209 0.3016141 #> F1PTQ9 1.585561e-02 1.0000000 0.55523576 0.3327302849 0.256601774 0.42717439 0.4494959 #> E2R8P1 1.000000e+00 1.0000000 0.26274295 0.02482534030.9457410 0.021207564 0.26828099 0.2724795 #> F6UWV4 1.482777e-06 0.9999998 0.00260364 0.0001465912 0.9495346 0.006502711 0.07246064 0.5000584 #> E2RQR5 9.610457e-01 1.0000000 0.18775918 0.0035986341 0.8871151 0.002950320 0.07612951 0.5850919
我们现在评估FGGA在等级f分数方面。分层分类性能度量,如分层精度(HP)、分层召回率(HR)和分层F-score (HF)度量出版物[3]适当地识别部分正确的分类,并相应地惩罚更遥远或更肤浅的错误。层次指标的公式如下:
\[惠普(s) = \压裂{1}{\ l中期(P_} {G (s))中期\}\水平间距{1.25毫米}\ sum_{问\水平间距{0.5毫米}\ \的水平间距{0.5毫米}l (P_} {G (s))} \水平间距{1.5毫米}\ max_ {c \水平间距{0.5毫米}\ \的水平间距{0.5毫米}l (C_} {G (s))} \压裂{\ \向上光标键c \中期水平间距{1毫米}\帽\向上光标键q \中期}{\向上光标键q} \]
\ [HR (s) = \压裂{1}{\ l中期(C_} {G (s))中期\}\水平间距{1.25毫米}\ sum_ {c \水平间距{0.5毫米}\ \的水平间距{0.5毫米}(C_} {G (s))} \水平间距{1.5毫米}\ max_{问\水平间距{0.5毫米}\ \的水平间距{0.5毫米}l (P_} {G (s))} \压裂{\ \向上光标键c \中期水平间距{1毫米}\帽\向上光标键q \中期}{\向上光标键c} \]
\[HF(s) = \frac{2 \cdot HP \cdot HR}{HP + HR} \]
HP HR HF样品0.5949843 0.7047748 0.6178593 2399
最后,我们对性能进行了评估FGGA在ROC曲线下的面积(AUC)和在x召回水平上的精度(PxR)。我们使用R包PerfMeas,它提供了计算我们需要的性能度量的函数。
library(PerfMeas) for (i in 1:dim(matrixFGGATest)[1]){matrixFGGATest[which(matrixFGGATest[i,] >= 0.5),] <- 1 matrixFGGATest[which(matrixFGGATest[i,] < 0.5),] <- 0} #计算F-score Fs <- f .度量。类(CfData$tableCfGO[rownames(matrixFGGATest),], matrixFGGATest) #平均F-score F $ Average [4] #> f# > 0.3070424 #计算AUC AUC <- AUC.single. overs .classes(CfData$tableCfGO[rownames(matrixFGGATest),], matrixFGGATest);#在不同召回水平下的计算精度PXR <- precision.at.multiple.recall.level.over。类(CfData$tableCfGO[rownames(matrixFGGATest),], matrixFGGATest, seq(from = 0.1, to = 1, by = 0.1));#平均PxR PxR $avgPXR #> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 #> 0.3546985 0.3417015 0.3515627 0.3419220 0.3232068 0.3124189 0.3129685 0.3144420 0.3009597 0.3030871
[1]张建军,李建军,李建军,张建军。基于因子图的围棋自动标注方法研究[j]。科学通报,2016,36(1):444 - 444。
[2]黄晓明,黄晓明,黄晓明,等。氧化石墨烯蛋白定位的研究进展[j]。科学通报,2018,77 (8)
[3]王晓明,王晓明,王晓明,等。一种面向本体函数标注的分类方法[J]。蛋白质科学,2006;15:1544-1549