GOCompare可以安装如下
# CRAN
install.packages( " GOCompare " )
# Alternative: GitHub
library( devtools )
remotes :: install_github( " ccsosa/GOCompare " )
下面包含该包所需的库的完整列表。
依赖项: R (>= 4.0.0)
导入: base, utils, methods, stats, grDevices, ape, vegan, ggplot2, ggrepel, igraph, parallel, stringr
建议: testthat
该 R 包提供了六个函数,以提供简单的工作流程来比较功能富集分析的结果:
功能: mostFrequentGOs. graphGOspecies
旨在提供对一种物种的分析。
功能: compareGOspecies graph_two_GOspecies evaluateCAT_species evaluateGO_species
允许比较属于用户所需类别的两个物种GO术语列表。
最后,包中提供了一组用于测试的四个数据集: A_thaliana, A_thaliana_compress, H_sapiens, H_sapiens_compress, comparison_example
功能富集分析结果
作为主要输入,您将需要两个数据框,其中包含来自您最喜欢的资源的功能丰富分析结果,例如:BinGO、AmiGO、ShinnyGO 或 TopGO。每个文件必须具有以下结构:
功能富集分析结果的data.frame
其中一列包含要分析的 GO 术语,一列包含要比较的类别
根据功能,您需要指定物种名称:species1 =“H_sapiens”和species2 =“A_thaliana”
必须提供包含要分析的 GO 术语所在的列名称的字段(例如:GOterm_field <- "Functional_Category")
功能_类别 | 特征 |
---|---|
对压力的反应 | 援助 |
防御反应 | 援助 |
细胞大小的调节 | 援助 |
防御反应 | 目的 |
对外部生物刺激的反应 | 大连商品交易所 |
gorth
函数的示例。require( gprofiler2 );require( stringr );require( GOCompare )
url_file = " https://raw.githubusercontent.com/ccsosa/R_Examples/master/Hallmarks_of_Cancer_AT.csv "
x <- read.csv( url_file )
x [, 1 ] <- NULL
CH <- c( " AID " , " AIM " , " DCE " , " ERI " , " EGS " , " GIM " , " IA " , " RCD " , " SPS " , " TPI " )
x_Hsap <- lapply(seq_len(length( CH )), function ( i ){
x_unique <- unique(na.omit( x [, i ]))
x_unique <- x_unique [which( x_unique != " " )]
x_unique <- as.list( x_unique )
return ( x_unique )
})
names( x_Hsap ) <- CH
# Using as background the unique genes for the ten CH.
GOterm_field <- " term_name "
x_s <- gprofiler2 :: gost( query = x_Hsap ,
organism = " hsapiens " , ordered_query = FALSE ,
multi_query = FALSE , significant = TRUE , exclude_iea = FALSE ,
measure_underrepresentation = FALSE , evcodes = FALSE ,
user_threshold = 0.05 , correction_method = " g_SCS " ,
domain_scope = " annotated " , custom_bg = unique(unlist( x_Hsap )),
numeric_ns = " " , sources = " GO:BP " , as_short_link = FALSE )
colnames( x_s $ result )[ 1 ] <- " feature "
# Check number of enriched terms per category
tapply( x_s $ result $ feature , x_s $ result $ feature , length )
# Running function to get graph of a list of features and GO terms
x <- graphGOspecies( df = x_s $ result ,
GOterm_field = GOterm_field ,
option = " Categories " ,
numCores = 1 ,
saveGraph = FALSE ,
outdir = NULL ,
filename = NULL )
# visualize nodes
View( x $ nodes )
# Get nodes with values greater than 95%
perc <- x $ nodes [which( x $ nodes $ WEIGHT > quantile( x $ nodes $ WEIGHT , probs = 0.95 )),]
# visualize nodes filtered
View( perc )
# ########
# Running function to get graph of a list of GO terms and categories
x_GO <- graphGOspecies( df = x_s $ result ,
GOterm_field = GOterm_field ,
option = " GO " ,
numCores = 1 ,
saveGraph = FALSE ,
outdir = NULL ,
filename = NULL )
# visualize nodes
View( x_GO $ nodes )
# Get GO terms nodes with values greater than 95%
perc_GO <- x_GO $ nodes [which( x_GO $ nodes $ GO_WEIGHT > quantile( x_GO $ nodes $ GO_WEIGHT , probs = 0.95 )),]
# visualize GO terms nodes filtered
View( perc_GO )
# #######################################################################################################
# two species comparison assuming they are the same genes in Drosophila melanogaster
orth_genes <- gprofiler2 :: gorth( query = unique(unlist( x_Hsap )), source_organism = " hsapiens " , target_organism = " dmelanogaster " )
# assigning genes
x_Dmap <- list ()
for ( i in 1 : length( x_Hsap )){
D_list <- list ()
for ( j in 1 : length( x_Hsap [[ i ]])){
x_orth <- orth_genes [ orth_genes $ input == x_Hsap [[ i ]][ j ],]
if (nrow( x_orth ) > 0 ){
D_list [[ j ]] <- data.frame ( orth = x_orth $ ortholog_name )
} else {
D_list [[ j ]] <- NULL
}
rm( x_orth )
};rm( j )
D_list <- unique(do.call( rbind , D_list ))
D_list <- D_list [which( ! is.null( D_list ))]
x_Dmap [[ i ]] <- D_list
rm( D_list )
};rm( i )
names( x_Dmap ) <- CH
GOterm_field <- " term_name "
x_s2 <- gprofiler2 :: gost( query = x_Dmap ,
organism = " dmelanogaster " , ordered_query = FALSE ,
multi_query = FALSE , significant = TRUE , exclude_iea = FALSE ,
measure_underrepresentation = FALSE , evcodes = FALSE ,
user_threshold = 0.05 , correction_method = " g_SCS " ,
domain_scope = " annotated " , custom_bg = unique(unlist( x_Dmap )),
numeric_ns = " " , sources = " GO:BP " , as_short_link = FALSE )
colnames( x_s2 $ result )[ 1 ] <- " feature "
# preparing input for compare two species
x_input <- GOCompare :: compareGOspecies( x_s $ result , x_s2 $ result , GOterm_field , species1 = " H. sapiens " , species2 = " D. melanogaster " , paired_lists = T )
# try to test similarities using clustering
plot(hclust( x_input $ distance , method = " ward.D " ))
# Comparing species results
comp_species_graph <- GOCompare :: graph_two_GOspecies( x_input , species1 = " H. sapiens " , species2 = " D. melanogaster " , option = " Categories " )
# View nodes order by combined weight (SPS and GIM categories have more frequent GO terms co-occurring)
View( comp_species_graph $ nodes [order( comp_species_graph $ nodes $ COMBINED_WEIGHT , decreasing = T ),])
comp_species_graph_GO <- GOCompare :: graph_two_GOspecies( x_input , species1 = " H. sapiens " , species2 = " D. melanogaster " , option = " GO " )
# Get GO terms nodes with values greater than 95%
perc_GO_two <- comp_species_graph_GO $ nodes [which( comp_species_graph_GO $ nodes $ GO_WEIGHT > quantile( comp_species_graph_GO $ nodes $ GO_WEIGHT , probs = 0.95 )),]
# visualize GO terms nodes filtered and ordered (more frequent GO terms in both species and categories)
View( perc_GO_two [order( perc_GO_two $ GO_WEIGHT , decreasing = T ),])
# evaluating if there are different in proportions of GO terms for each category
x_CAT <- GOCompare :: evaluateCAT_species( x_s $ result , x_s2 $ result , species1 = " H. sapiens " , species2 = " D. melanogaster " , GOterm_field = " term_name " , test = " prop " )
x_CAT <- x_CAT [which( x_CAT $ FDR < = 0.05 ),]
# View Categories with FDR <0.05 (RCD,SPS,GIM, AIM,ERI,DCE)
View( x_CAT )
# evaluating if there are different in proportions of categories for GO terms
x_GO <- GOCompare :: evaluateGO_species( x_s $ result , x_s2 $ result , species1 = " H. sapiens " , species2 = " D. melanogaster " , GOterm_field = " term_name " , test = " prop " )
x_GO <- x_GO [which( x_GO $ FDR < = 0.05 ),]
# View Categories with FDR <0.05 (No significant results in proportions)
View( x_GO )
# #Optional plots (omit # symbol and run)
# source("https://raw.githubusercontent.com/ccsosa/Supplementary-information/refs/heads/main/CHAPTER3/PLOT_TWO_SP_GRAPH_CAT.R")
# source("https://raw.githubusercontent.com/ccsosa/Supplementary-information/refs/heads/main/CHAPTER3/PLOT_TWO_SP_GRAPH_GO.R")
# plot_twosp_CAT("D:/",comp_species_graph)
# plot_twosp_GO("D:/",comp_species_graph_GO)
主要成员:克里斯蒂安·C·索萨、戴安娜·卡罗莱纳·克拉维霍-布里蒂卡、毛里西奥·金巴亚、维克多·雨果·加西亚-梅尔查
其他贡献者:Nicolas Lopéz-Rozo、Camila Riccio Rengifo、David Arango Londoño、Maria Victoria Diaz
索萨、克里斯蒂安·C.、戴安娜·卡罗莱纳·克拉维霍-布里蒂卡、维克多·雨果·加西亚-梅尔钱、尼古拉斯·洛佩斯-罗佐、卡米拉·里奇奥-伦吉福、玛丽亚·维多利亚·迪亚兹、大卫·阿兰戈·隆多尼奥、毛里西奥·阿尔贝托·金巴亚。 «GOCompare:用于比较两个物种之间功能富集分析的 R 包»。基因组学 115,第 1 期(2023 年 1 月):110528。https://doi.org/10.1016/j.ygeno.2022.110528。
GNU 通用公共许可证版本 3