| Title: | Generalized Reporter Score-Based Enrichment Analysis for Omics Data |
|---|---|
| Description: | Inspired by the classic 'RSA', we developed the improved 'Generalized Reporter Score-based Analysis (GRSA)' method, implemented in the R package 'ReporterScore', along with comprehensive visualization methods and pathway databases. 'GRSA' is a threshold-free method that works well with all types of biomedical features, such as genes, chemical compounds, and microbial species. Importantly, the 'GRSA' supports multi-group and longitudinal experimental designs, because of the included multi-group-compatible statistical methods. |
| Authors: | Chen Peng [aut, cre] |
| Maintainer: | Chen Peng <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.9 |
| Built: | 2024-11-28 16:53:42 UTC |
| Source: | CRAN |
Test the proper clusters k for c_means
C-means cluster
cm_test_k(otu_group, filter_var, fast = TRUE) c_means(otu_group, k_num, filter_var)cm_test_k(otu_group, filter_var, fast = TRUE) c_means(otu_group, k_num, filter_var)
otu_group |
standardize data |
filter_var |
filter the highest var |
fast |
whether do the gap_stat? |
k_num |
cluster number |
ggplot
ggplot
Other C_means:
RSA_by_cm()
if (requireNamespace("e1071") && requireNamespace("factoextra")) { data(otutab, package = "pcutils") pcutils::hebing(otutab, metadata$Group) -> otu_group cm_test_k(otu_group, filter_var = 0.7) cm_res <- c_means(otu_group, k_num = 3, filter_var = 0.7) plot(cm_res, 0.8) }if (requireNamespace("e1071") && requireNamespace("factoextra")) { data(otutab, package = "pcutils") pcutils::hebing(otutab, metadata$Group) -> otu_group cm_test_k(otu_group, filter_var = 0.7) cm_res <- c_means(otu_group, k_num = 3, filter_var = 0.7) plot(cm_res, 0.8) }
Combine the results of 'step by step GRSA'
combine_rs_res(kodf, group, metadata, ko_stat, reporter_s, modulelist = NULL)combine_rs_res(kodf, group, metadata, ko_stat, reporter_s, modulelist = NULL)
kodf |
KO_abundance table, rowname are feature ids (e.g. K00001 if feature="ko"; PEX11A if feature="gene"; C00024 if feature="compound"), colnames are samples. |
group |
The comparison groups (at least two categories) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. And you can use factor levels to change order. |
metadata |
sample information data.frame contains group |
ko_stat |
result of |
reporter_s |
result of |
modulelist |
NULL or customized modulelist dataframe, must contain 'id','K_num','KOs','Description' columns. Take the 'KOlist' as example, use |
reporter_score object
Other GRSA:
get_reporter_score(),
ko.test(),
pvalue2zs(),
reporter_score()
data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed") reporter_s1 <- get_reporter_score(ko_stat, perm = 499) reporter_res <- combine_rs_res(KO_abundance, "Group", metadata, ko_stat, reporter_s1)data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed") reporter_s1 <- get_reporter_score(ko_stat, perm = 499) reporter_res <- combine_rs_res(KO_abundance, "Group", metadata, ko_stat, reporter_s1)
Compound htable from 'KEGG'
Other data:
CPDlist,
GOlist,
KO_htable,
KOlist,
Module_htable,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
an list contains two data.frame named pathway and module.
four columns in each data.frame.
"map0010" or "M00001"
contians how many Compounds in this pathway or module
Compounds name
the description of this pathway or module
Other data:
Compound_htable,
GOlist,
KO_htable,
KOlist,
Module_htable,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
Build a custom modulelist
Transform a modulelist to a list
custom_modulelist(pathway2ko, pathway2desc = NULL, verbose = TRUE) transform_modulelist(mymodulelist, mode = 1)custom_modulelist(pathway2ko, pathway2desc = NULL, verbose = TRUE) transform_modulelist(mymodulelist, mode = 1)
pathway2ko |
user input annotation of Pathway to KO mapping, a data.frame of 2 column with pathway and ko. |
pathway2desc |
user input of Pathway TO Description mapping, a data.frame of 2 column with pathway and description. |
verbose |
verbose |
mymodulelist |
mymodulelist |
mode |
1~2 |
a custom modulelist
modulelist
Other modulelist:
custom_modulelist_from_org(),
get_features()
Other modulelist:
custom_modulelist_from_org(),
get_features()
mydat <- data.frame(pathway = paste0("PATHWAY", rep(seq_len(2), each = 5)), ko = paste0("K", 1:10)) mymodulelist <- custom_modulelist(mydat) print(mymodulelist) transform_modulelist(mymodulelist)mydat <- data.frame(pathway = paste0("PATHWAY", rep(seq_len(2), each = 5)), ko = paste0("K", 1:10)) mymodulelist <- custom_modulelist(mydat) print(mymodulelist) transform_modulelist(mymodulelist)
Custom modulelist from a specific organism
custom_modulelist_from_org( org = "hsa", feature = "ko", gene = "symbol", verbose = TRUE )custom_modulelist_from_org( org = "hsa", feature = "ko", gene = "symbol", verbose = TRUE )
org |
kegg organism, listed in https://www.genome.jp/kegg/catalog/org_list.html, default, "hsa" |
feature |
one of "ko", "gene", "compound" |
gene |
one of "symbol","id" |
verbose |
logical |
modulelist
Other modulelist:
custom_modulelist(),
get_features()
hsa_pathway <- custom_modulelist_from_org(org = "hsa", feature = "gene")hsa_pathway <- custom_modulelist_from_org(org = "hsa", feature = "gene")
Export report score result tables
export_report_table(reporter_res, dir_name, overwrite = FALSE)export_report_table(reporter_res, dir_name, overwrite = FALSE)
reporter_res |
a reporter_score object or rs_by_cm object |
dir_name |
the directory to save the report tables |
overwrite |
overwrite the existed files or not, default is FALSE. |
No return value
You can use 'clusterProfiler::bitr()' to transfer your table from other gene_id to gene_symbol.
gene2ko(genedf, org = "hsa")gene2ko(genedf, org = "hsa")
genedf |
,rowname is gene symbol (e.g. PFKM), colnames is samples |
org |
kegg organism, listed in 'https://www.genome.jp/kegg/catalog/org_list.html', default, 'hsa' |
kodf
data("genedf") KOdf <- gene2ko(genedf, org = "hsa")data("genedf") KOdf <- gene2ko(genedf, org = "hsa")
human gene table
Other test_data:
KO_abundance,
reporter_score_res
get features in a modulelist
get_features(map_id = "map00010", ko_stat = NULL, modulelist = NULL)get_features(map_id = "map00010", ko_stat = NULL, modulelist = NULL)
map_id |
map_id in modulelist |
ko_stat |
NULL or ko_stat result from |
modulelist |
NULL or customized modulelist dataframe, must contain 'id','K_num','KOs','Description' columns. Take the 'KOlist' as example, use |
KOids, or data.frame with these KOids.
Other modulelist:
custom_modulelist_from_org(),
custom_modulelist()
get_features(map_id = "map00010")get_features(map_id = "map00010")
Calculate reporter score
get_reporter_score( ko_stat, type = c("pathway", "module")[1], feature = "ko", threads = 1, modulelist = NULL, perm = 4999, verbose = TRUE, p.adjust.method2 = "BH", min_exist_KO = 3, max_exist_KO = 600 )get_reporter_score( ko_stat, type = c("pathway", "module")[1], feature = "ko", threads = 1, modulelist = NULL, perm = 4999, verbose = TRUE, p.adjust.method2 = "BH", min_exist_KO = 3, max_exist_KO = 600 )
ko_stat |
ko_stat result from |
type |
'pathway' or 'module' for default KOlist for microbiome, 'CC', 'MF', 'BP', 'ALL' for default GOlist for homo sapiens. And org in listed in 'https://www.genome.jp/kegg/catalog/org_list.html' such as 'hsa' (if your kodf is come from a specific organism, you should specify type here). |
feature |
one of 'ko', 'gene', 'compound' |
threads |
default 1 |
modulelist |
NULL or customized modulelist dataframe, must contain 'id','K_num','KOs','Description' columns. Take the 'KOlist' as example, use |
perm |
permutation number, default: 4999. |
verbose |
logical |
p.adjust.method2 |
p.adjust.method for the correction of ReporterScore, see |
min_exist_KO |
min exist KO number in a pathway (default, 3, when a pathway contains KOs less than 3, there will be no RS) |
max_exist_KO |
max exist KO number in a pathway (default, 600, when a pathway contains KOs more than 600, there will be no RS) |
reporter_res data.frame
Other GRSA:
combine_rs_res(),
ko.test(),
pvalue2zs(),
reporter_score()
data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed") reporter_s1 <- get_reporter_score(ko_stat, perm = 499)data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed") reporter_s1 <- get_reporter_score(ko_stat, perm = 499)
an list contains three data.frame named BP, CC, MF.
four columns in each data.frame.
"map0010" or "M00001"
contians how many Genes in this GO term
Genes name
the description of this GO term
Other data:
CPDlist,
Compound_htable,
KO_htable,
KOlist,
Module_htable,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
pathway information for "hsa"
Other data:
CPDlist,
Compound_htable,
GOlist,
KO_htable,
KOlist,
Module_htable,
Pathway_htable,
mmu_kegg_pathway
The KOs abundance table and group table.
The KOs abundance table and group table.
Other test_data:
genedf,
reporter_score_res
This function performs KO enrichment analysis using the 'clusterProfiler' package.
KO_enrich( ko_stat, padj_threshold = 0.05, logFC_threshold = NULL, add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE ) as.enrich_res(gsea_res)KO_enrich( ko_stat, padj_threshold = 0.05, logFC_threshold = NULL, add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE ) as.enrich_res(gsea_res)
ko_stat |
ko_stat dataframe from |
padj_threshold |
p.adjust threshold to determine whether a feature significant or not. p.adjust < padj_threshold, default: 0.05 |
logFC_threshold |
logFC threshold to determine whether a feature significant or not. abs(logFC)>logFC_threshold, default: NULL |
add_mini |
add_mini when calculate the logFC. e.g (10+0.1)/(0+0.1), default 0.05*min(avg_abundance) |
p.adjust.method |
The method used for p-value adjustment (default: "BH"). |
type |
"pathway" or "module" for default KOlist_file. |
feature |
one of "ko", "gene", "compound" |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
verbose |
logical |
gsea_res |
gsea_res from KO_gsea |
A data frame containing the enrichment results.
enrich_res object
Other common_enrich:
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_safe(),
KO_sea(),
plot_enrich_res()
Perform fisher's exact enrichment analysis
KO_fisher( ko_stat, padj_threshold = 0.05, logFC_threshold = NULL, add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE )KO_fisher( ko_stat, padj_threshold = 0.05, logFC_threshold = NULL, add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE )
ko_stat |
ko_stat dataframe from |
padj_threshold |
p.adjust threshold to determine whether a feature significant or not. p.adjust < padj_threshold, default: 0.05 |
logFC_threshold |
logFC threshold to determine whether a feature significant or not. abs(logFC)>logFC_threshold, default: NULL |
add_mini |
add_mini when calculate the logFC. e.g (10+0.1)/(0+0.1), default 0.05*min(avg_abundance) |
p.adjust.method |
The method used for p-value adjustment (default: "BH"). |
type |
"pathway" or "module" for default KOlist_file. |
feature |
one of "ko", "gene", "compound" |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
verbose |
logical |
data.frame
Other common_enrich:
KO_enrich(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_safe(),
KO_sea(),
plot_enrich_res()
## use `fisher.test` from the `stats` package. data("reporter_score_res") fisher_res <- KO_fisher(reporter_score_res)## use `fisher.test` from the `stats` package. data("reporter_score_res") fisher_res <- KO_fisher(reporter_score_res)
Perform gene set analysis
KO_gsa( reporter_res, method = "Two class unpaired", p.adjust.method = "BH", verbose = TRUE, perm = 1000, ... )KO_gsa( reporter_res, method = "Two class unpaired", p.adjust.method = "BH", verbose = TRUE, perm = 1000, ... )
reporter_res |
reporter_res |
method |
Problem type: "quantitative" for a continuous parameter; "Two class unpaired" ; "Survival" for censored survival outcome; "Multiclass" : more than 2 groups, coded 1,2,3...; "Two class paired" for paired outcomes, coded -1,1 (first pair), -2,2 (second pair), etc |
p.adjust.method |
"BH" |
verbose |
TRUE |
perm |
1000 |
... |
additional parameters to |
enrich_res object
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_safe(),
KO_sea(),
plot_enrich_res()
## use `GSA` from the `GSA` package. if (requireNamespace("GSA")) { data("reporter_score_res") gsa_res <- KO_gsa(reporter_score_res, p.adjust.method = "none", perm = 200) plot(gsa_res) }## use `GSA` from the `GSA` package. if (requireNamespace("GSA")) { data("reporter_score_res") gsa_res <- KO_gsa(reporter_score_res, p.adjust.method = "none", perm = 200) plot(gsa_res) }
Perform gene set enrichment analysis
KO_gsea( ko_stat, weight = "logFC", add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE )KO_gsea( ko_stat, weight = "logFC", add_mini = NULL, p.adjust.method = "BH", type = c("pathway", "module")[1], feature = "ko", modulelist = NULL, verbose = TRUE )
ko_stat |
ko_stat dataframe from |
weight |
the metric used for ranking, default: logFC |
add_mini |
add_mini when calculate the logFC. e.g (10+0.1)/(0+0.1), default 0.05*min(avg_abundance) |
p.adjust.method |
The method used for p-value adjustment (default: "BH"). |
type |
"pathway" or "module" for default KOlist_file. |
feature |
one of "ko", "gene", "compound" |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
verbose |
logical |
DOSE object
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsva(),
KO_padog(),
KO_safe(),
KO_sea(),
plot_enrich_res()
message("The following example require some time to run:") ## use `GSEA` from the `clusterProfiler` package. if (requireNamespace("clusterProfiler")) { data("reporter_score_res") gsea_res <- KO_gsea(reporter_score_res, p.adjust.method = "none") enrichplot::gseaplot(gsea_res, geneSetID = data.frame(gsea_res)$ID[1]) gsea_res_df <- as.enrich_res(gsea_res) plot(gsea_res_df) }message("The following example require some time to run:") ## use `GSEA` from the `clusterProfiler` package. if (requireNamespace("clusterProfiler")) { data("reporter_score_res") gsea_res <- KO_gsea(reporter_score_res, p.adjust.method = "none") enrichplot::gseaplot(gsea_res, geneSetID = data.frame(gsea_res)$ID[1]) gsea_res_df <- as.enrich_res(gsea_res) plot(gsea_res_df) }
Perform Gene Set Variation Analysis
KO_gsva( reporter_res, verbose = TRUE, method = "wilcox.test", p.adjust.method = "BH", ... )KO_gsva( reporter_res, verbose = TRUE, method = "wilcox.test", p.adjust.method = "BH", ... )
reporter_res |
reporter_res |
verbose |
verbose |
method |
see |
p.adjust.method |
p.adjust.method |
... |
additional parameters to |
enrich_res
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_padog(),
KO_safe(),
KO_sea(),
plot_enrich_res()
## use `gsva` from the `GSVA` package. if (requireNamespace("GSVA")) { data("reporter_score_res") gsva_res <- KO_gsva(reporter_score_res, p.adjust.method = "none") }## use `gsva` from the `GSVA` package. if (requireNamespace("GSVA")) { data("reporter_score_res") gsva_res <- KO_gsva(reporter_score_res, p.adjust.method = "none") }
KO htable from 'KEGG'
Other data:
CPDlist,
Compound_htable,
GOlist,
KOlist,
Module_htable,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
Perform Pathway Analysis with Down-weighting of Overlapping Genes (PADOG)
KO_padog( reporter_res, verbose = TRUE, perm = 1000, p.adjust.method = "BH", ... )KO_padog( reporter_res, verbose = TRUE, perm = 1000, p.adjust.method = "BH", ... )
reporter_res |
The input reporter result. |
verbose |
If TRUE, print verbose messages. Default is TRUE. |
perm |
The number of permutations. Default is 1000. |
p.adjust.method |
Method for p-value adjustment. Default is "BH". |
... |
Additional parameters to be passed to |
A data frame containing PADOG results for KO enrichment.
A data frame with columns "ID," "Description," "K_num," "Exist_K_num," "p.value," and "p.adjust."
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_safe(),
KO_sea(),
plot_enrich_res()
## use `PADOG` from the `PADOG` package. if (requireNamespace("PADOG")) { data("reporter_score_res") padog_res <- KO_padog(reporter_score_res, verbose = TRUE, perm = 200, p.adjust.method = "none" ) }## use `PADOG` from the `PADOG` package. if (requireNamespace("PADOG")) { data("reporter_score_res") padog_res <- KO_padog(reporter_score_res, verbose = TRUE, perm = 200, p.adjust.method = "none" ) }
Perform Significance Analysis of Function and Expression
KO_safe( reporter_res, verbose = TRUE, perm = 1000, C.matrix = NULL, p.adjust.method = "BH", ... )KO_safe( reporter_res, verbose = TRUE, perm = 1000, C.matrix = NULL, p.adjust.method = "BH", ... )
reporter_res |
The input reporter result. |
verbose |
If TRUE, print verbose messages. Default is TRUE. |
perm |
The number of permutations. Default is 1000. |
C.matrix |
The contrast matrix. Default is NULL, and it will be generated from the module list. |
p.adjust.method |
Method for p-value adjustment. Default is "BH". |
... |
Additional parameters to be passed to |
A data frame containing SAFE results for KO enrichment.
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_sea(),
plot_enrich_res()
## use `safe` from the `safe` package. if (requireNamespace("safe")) { data("reporter_score_res") safe_res <- KO_safe(reporter_score_res, verbose = TRUE, perm = 200, p.adjust.method = "none" ) }## use `safe` from the `safe` package. if (requireNamespace("safe")) { data("reporter_score_res") safe_res <- KO_safe(reporter_score_res, verbose = TRUE, perm = 200, p.adjust.method = "none" ) }
Perform Simultaneous Enrichment Analysis
KO_sea(reporter_res, verbose = TRUE, ...)KO_sea(reporter_res, verbose = TRUE, ...)
reporter_res |
The input reporter result. |
verbose |
If TRUE, print verbose messages. Default is TRUE. |
... |
Additional parameters to be passed to |
enrich_res
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_safe(),
plot_enrich_res()
## use `SEA` from the `rSEA` package. if (requireNamespace("rSEA")) { data("reporter_score_res") sea_res <- KO_sea(reporter_score_res, verbose = TRUE) }## use `SEA` from the `rSEA` package. if (requireNamespace("rSEA")) { data("reporter_score_res") sea_res <- KO_sea(reporter_score_res, verbose = TRUE) }
Differential analysis or Correlation analysis for KO-abundance table
ko.test( kodf, group, metadata = NULL, method = "wilcox.test", pattern = NULL, p.adjust.method1 = "none", threads = 1, verbose = TRUE )ko.test( kodf, group, metadata = NULL, method = "wilcox.test", pattern = NULL, p.adjust.method1 = "none", threads = 1, verbose = TRUE )
kodf |
KO_abundance table, rowname are feature ids (e.g. K00001 if feature="ko"; PEX11A if feature="gene"; C00024 if feature="compound"), colnames are samples. |
group |
The comparison groups (at least two categories) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. And you can use factor levels to change order. |
metadata |
sample information data.frame contains group |
method |
the type of test. Default is 'wilcox.test'. Allowed values include:
|
pattern |
a named vector matching the group, e.g. c('G1'=1,'G2'=3,'G3'=2), use the correlation analysis with specific pattern to calculate p-value. |
p.adjust.method1 |
p.adjust.method for 'ko.test', see |
threads |
default 1 |
verbose |
logical |
ko_pvalue data.frame
Other GRSA:
combine_rs_res(),
get_reporter_score(),
pvalue2zs(),
reporter_score()
data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata)data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata)
an list contains two data.frame named pathway and module.
four columns in each data.frame.
"map0010" or "M00001"
contians how many KOs in this pathway or module
KOs name
the description of this pathway or module
Other data:
CPDlist,
Compound_htable,
GOlist,
KO_htable,
Module_htable,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
Load the CARDinfo (from CARD database)
load_CARDinfo(verbose = TRUE)load_CARDinfo(verbose = TRUE)
verbose |
logical |
CARDinfo
Load the GOlist (from 'GO' database)
Load the GOinfo (from GO)
load_GOlist(verbose = TRUE) load_GOinfo(verbose = TRUE)load_GOlist(verbose = TRUE) load_GOinfo(verbose = TRUE)
verbose |
logical |
GOlist
GOinfo
Load the specific table (from 'KEGG')
Load the KOlist (from 'KEGG')
Load the CPDlist (from 'KEGG')
Load the KO description (from 'KEGG')
Load the KO_htable (from 'KEGG')
Load the Pathway_htable (from 'KEGG')
Load the Module_htable (from 'KEGG')
Load the Compound_htable (from 'KEGG')
Load the pathway information for an organism (from 'KEGG')
load_htable(type, verbose = TRUE) load_KOlist(verbose = TRUE) load_CPDlist(verbose = TRUE) load_KO_desc(verbose = TRUE) load_KO_htable(verbose = TRUE) load_Pathway_htable(verbose = TRUE) load_Module_htable(verbose = TRUE) load_Compound_htable(verbose = TRUE) load_org_pathway(org = "hsa", verbose = TRUE)load_htable(type, verbose = TRUE) load_KOlist(verbose = TRUE) load_CPDlist(verbose = TRUE) load_KO_desc(verbose = TRUE) load_KO_htable(verbose = TRUE) load_Pathway_htable(verbose = TRUE) load_Module_htable(verbose = TRUE) load_Compound_htable(verbose = TRUE) load_org_pathway(org = "hsa", verbose = TRUE)
type |
"ko", "module", "pathway", "compound" ... |
verbose |
logical |
org |
kegg organism, listed in https://www.genome.jp/kegg/catalog/org_list.html, default, "hsa" |
KO_htable
KOlist
CPDlist
KO description
KO_htable
Pathway_htable
Module_htable
Compound_htable
KOlist
Pathway_htable <- load_htable("pathway") head(Pathway_htable)Pathway_htable <- load_htable("pathway") head(Pathway_htable)
pathway information for "mmu"
Other data:
CPDlist,
Compound_htable,
GOlist,
KO_htable,
KOlist,
Module_htable,
Pathway_htable,
hsa_kegg_pathway
Modify the pathway description before plotting
modify_description( reporter_res, pattern = " - Homo sapiens (human)", replacement = "" )modify_description( reporter_res, pattern = " - Homo sapiens (human)", replacement = "" )
reporter_res |
reporter_res |
pattern |
str, like " - Homo sapiens (human)" |
replacement |
str, like "" |
reporter_res
data("reporter_score_res") modify_description(reporter_score_res, pattern = " - Homo sapiens (human)")data("reporter_score_res") modify_description(reporter_score_res, pattern = " - Homo sapiens (human)")
Module htable from 'KEGG'
Other data:
CPDlist,
Compound_htable,
GOlist,
KO_htable,
KOlist,
Pathway_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
Pathway htable from 'KEGG'
Other data:
CPDlist,
Compound_htable,
GOlist,
KO_htable,
KOlist,
Module_htable,
hsa_kegg_pathway,
mmu_kegg_pathway
Plot enrich_res
Plot enrich_res
plot_enrich_res( enrich_res, mode = 1, padj_threshold = 0.05, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, str_width = 50, facet_str_width = 15, ... ) ## S3 method for class 'enrich_res' plot( x, mode = 1, padj_threshold = 0.05, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, str_width = 50, facet_str_width = 15, ... )plot_enrich_res( enrich_res, mode = 1, padj_threshold = 0.05, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, str_width = 50, facet_str_width = 15, ... ) ## S3 method for class 'enrich_res' plot( x, mode = 1, padj_threshold = 0.05, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, str_width = 50, facet_str_width = 15, ... )
enrich_res |
enrich_res object |
mode |
plot style: 1~2 |
padj_threshold |
p.adjust threshold |
show_ID |
show pathway id |
Pathway_description |
show KO description rather than KO id. |
facet_level |
facet plot if the type is "pathway" or "module" |
facet_anno |
annotation table for facet, two columns, first is level summary, second is pathway id. |
str_width |
default: 50 |
facet_str_width |
str width for facet label |
... |
add |
x |
enrich_res object |
ggplot
ggplot
Other common_enrich:
KO_enrich(),
KO_fisher(),
KO_gsa(),
KO_gsea(),
KO_gsva(),
KO_padog(),
KO_safe(),
KO_sea()
Plot features boxplot
plot_features_box( kodf, group = NULL, metadata = NULL, map_id = "map00780", select_ko = NULL, only_sig = FALSE, box_param = NULL, modulelist = NULL, KO_description = FALSE, str_width = 50 )plot_features_box( kodf, group = NULL, metadata = NULL, map_id = "map00780", select_ko = NULL, only_sig = FALSE, box_param = NULL, modulelist = NULL, KO_description = FALSE, str_width = 50 )
kodf |
KO_abundance table, rowname is ko id (e.g. K00001),colnames is samples. or result of 'get_reporter_score' |
group |
The compare group (two category) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. |
metadata |
metadata |
map_id |
the pathway or module id |
select_ko |
select which ko |
only_sig |
only show the significant features |
box_param |
parameters pass to |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
KO_description |
show KO description rather than KO id. |
str_width |
str_width to wrap |
ggplot
data("reporter_score_res") plot_features_box(reporter_score_res, select_ko = c("K00059", "K00208", "K00647", "K00652", "K00833", "K01012"), box_param = list(p_value1 = FALSE, trend_line = TRUE) ) plot_features_box(reporter_score_res, select_ko = "K00059", KO_description = TRUE, box_param = list(p_value1 = FALSE, trend_line = TRUE) )data("reporter_score_res") plot_features_box(reporter_score_res, select_ko = c("K00059", "K00208", "K00647", "K00652", "K00833", "K01012"), box_param = list(p_value1 = FALSE, trend_line = TRUE) ) plot_features_box(reporter_score_res, select_ko = "K00059", KO_description = TRUE, box_param = list(p_value1 = FALSE, trend_line = TRUE) )
plot the Z-score of features distribution
plot_features_distribution( reporter_res, map_id, text_size = 4, text_position = NULL, rug_length = 0.04 )plot_features_distribution( reporter_res, map_id, text_size = 4, text_position = NULL, rug_length = 0.04 )
reporter_res |
result of 'reporter_score' |
map_id |
the pathway or module id |
text_size |
text_size=4 |
text_position |
text_position, e.g. c(x=3,y=0.4) |
rug_length |
rug_length=0.04 |
ggplot
data("reporter_score_res") plot_features_distribution(reporter_score_res, map_id = c("map05230", "map03010"))data("reporter_score_res") plot_features_distribution(reporter_score_res, map_id = c("map05230", "map03010"))
Plot features heatmap
plot_features_heatmap( kodf, group = NULL, metadata = NULL, map_id = "map00780", select_ko = NULL, only_sig = FALSE, columns = NULL, modulelist = NULL, KO_description = FALSE, str_width = 50, heatmap_param = list() )plot_features_heatmap( kodf, group = NULL, metadata = NULL, map_id = "map00780", select_ko = NULL, only_sig = FALSE, columns = NULL, modulelist = NULL, KO_description = FALSE, str_width = 50, heatmap_param = list() )
kodf |
KO_abundance table, rowname is ko id (e.g. K00001),colnames is samples. or result of 'get_reporter_score' |
group |
The compare group (two category) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. |
metadata |
metadata |
map_id |
the pathway or module id |
select_ko |
select which ko |
only_sig |
only show the significant KOs |
columns |
change columns |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
KO_description |
show KO description rather than KO id. |
str_width |
str_width to wrap |
heatmap_param |
parameters pass to |
ggplot
if (requireNamespace("pheatmap")) { data("reporter_score_res") plot_features_heatmap(reporter_score_res, map_id = "map00780") }if (requireNamespace("pheatmap")) { data("reporter_score_res") plot_features_heatmap(reporter_score_res, map_id = "map00780") }
Plot features trend in one pathway or module
plot_features_in_pathway( ko_stat, map_id = "map00780", modulelist = NULL, select_ko = NULL, box_color = reporter_color, show_number = TRUE, scale = FALSE, feature_type = "KOs", line_color = c(Depleted = "seagreen", Enriched = "orange", None = "grey", Significant = "red2") )plot_features_in_pathway( ko_stat, map_id = "map00780", modulelist = NULL, select_ko = NULL, box_color = reporter_color, show_number = TRUE, scale = FALSE, feature_type = "KOs", line_color = c(Depleted = "seagreen", Enriched = "orange", None = "grey", Significant = "red2") )
ko_stat |
ko_stat result from |
map_id |
the pathway or module id |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
select_ko |
select which ko |
box_color |
box and point color, default: c("#e31a1c","#1f78b4") |
show_number |
show the numbers. |
scale |
scale the data by row. |
feature_type |
show in the title ,default: KOs |
line_color |
line color, default: c("Depleted"="seagreen","Enriched"="orange","None"="grey") |
ggplot
data("reporter_score_res") plot_features_in_pathway(ko_stat = reporter_score_res, map_id = "map00860")data("reporter_score_res") plot_features_in_pathway(ko_stat = reporter_score_res, map_id = "map00860")
Plot features network
plot_features_network( ko_stat, map_id = "map00780", near_pathway = FALSE, modulelist = NULL, kos_color = c(Depleted = "seagreen", Enriched = "orange", None = "grey", Significant = "red2", Pathway = "#80b1d3"), pathway_label = TRUE, kos_label = TRUE, pathway_description = FALSE, kos_description = FALSE, str_width = 50, mark_module = FALSE, mark_color = NULL, return_net = FALSE, ... )plot_features_network( ko_stat, map_id = "map00780", near_pathway = FALSE, modulelist = NULL, kos_color = c(Depleted = "seagreen", Enriched = "orange", None = "grey", Significant = "red2", Pathway = "#80b1d3"), pathway_label = TRUE, kos_label = TRUE, pathway_description = FALSE, kos_description = FALSE, str_width = 50, mark_module = FALSE, mark_color = NULL, return_net = FALSE, ... )
ko_stat |
ko_stat result from |
map_id |
the pathway or module id |
near_pathway |
show the near_pathway if any features exist. |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
kos_color |
default, c("Depleted"="seagreen","Enriched"="orange","None"="grey","Significant"="red2") |
pathway_label |
show pathway_label? |
kos_label |
show kos_label? |
pathway_description |
show the pathway description? |
kos_description |
show the kos description? |
str_width |
str width |
mark_module |
mark the modules? |
mark_color |
mark colors, default, c("Depleted"="seagreen","Enriched"="orange","None"="grey","Significant"="red2") |
return_net |
return the network |
... |
additional arguments for |
network plot
if (requireNamespace("MetaNet")) { data("reporter_score_res") plot_features_network(reporter_score_res, map_id = "map05230") plot_features_network(reporter_score_res, map_id = "map00780", near_pathway = TRUE) }if (requireNamespace("MetaNet")) { data("reporter_score_res") plot_features_network(reporter_score_res, map_id = "map05230") plot_features_network(reporter_score_res, map_id = "map00780", near_pathway = TRUE) }
Plot htable levels
plot_htable(type = "ko", select = NULL, htable = NULL)plot_htable(type = "ko", select = NULL, htable = NULL)
type |
"ko", "module", "pathway", "compound" |
select |
select ids |
htable |
custom a htable |
ggplot
data("KO_abundance_test") plot_htable(select = rownames(KO_abundance))data("KO_abundance_test") plot_htable(select = rownames(KO_abundance))
plot_KEGG_map
plot_KEGG_map( ko_stat, map_id = "map00780", modulelist = NULL, type = "pathway", feature = "ko", color_var = "Z_score", save_dir, color = c("seagreen", "grey", "orange") )plot_KEGG_map( ko_stat, map_id = "map00780", modulelist = NULL, type = "pathway", feature = "ko", color_var = "Z_score", save_dir, color = c("seagreen", "grey", "orange") )
ko_stat |
ko_stat result from |
map_id |
the pathway or module id |
modulelist |
NULL or customized modulelist dataframe, must contain "id","K_num","KOs","Description" columns. Take the 'KOlist' as example, use |
type |
"pathway" or "module" for default KOlist for microbiome, "CC", "MF", "BP", "ALL" for default GOlist for homo sapiens. And org in listed in 'https://www.genome.jp/kegg/catalog/org_list.html' such as "hsa" (if your kodf is come from a specific organism, you should specify type here). |
feature |
one of "ko", "gene", "compound" |
color_var |
use which variable to color |
save_dir |
where to save the png files |
color |
color |
png files
https://zhuanlan.zhihu.com/p/357687076
message("The following example will download some files, run yourself:") if (requireNamespace("pathview")) { output_dir <- tempdir() data("reporter_score_res") plot_KEGG_map(reporter_score_res$ko_stat, map_id = "map00780", type = "pathway", feature = "ko", color_var = "Z_score", save_dir = output_dir ) }message("The following example will download some files, run yourself:") if (requireNamespace("pathview")) { output_dir <- tempdir() data("reporter_score_res") plot_KEGG_map(reporter_score_res$ko_stat, map_id = "map00780", type = "pathway", feature = "ko", color_var = "Z_score", save_dir = output_dir ) }
Plot the reporter_res
plot_report( reporter_res, rs_threshold = 1.64, mode = 1, y_text_size = 13, str_width = 100, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, facet_str_width = 15, plot_line = TRUE, reorder = FALSE )plot_report( reporter_res, rs_threshold = 1.64, mode = 1, y_text_size = 13, str_width = 100, show_ID = FALSE, Pathway_description = TRUE, facet_level = FALSE, facet_anno = NULL, facet_str_width = 15, plot_line = TRUE, reorder = FALSE )
reporter_res |
result of 'get_reporter_score' or 'reporter_score' |
rs_threshold |
plot threshold vector, default:1.64 |
mode |
1~3 plot style. |
y_text_size |
y_text_size |
str_width |
str_width to wrap |
show_ID |
show pathway id |
Pathway_description |
show KO description rather than KO id. |
facet_level |
facet plot if the type is "pathway" or "module" |
facet_anno |
annotation table for facet, two columns, first is level summary, second is pathway id. |
facet_str_width |
str width for facet label |
plot_line |
plot line or not |
reorder |
reorder the order of the pathways |
ggplot
data("reporter_score_res") plot_report(reporter_score_res, rs_threshold = c(2.5, -2.5), y_text_size = 10, str_width = 40)data("reporter_score_res") plot_report(reporter_score_res, rs_threshold = c(2.5, -2.5), y_text_size = 10, str_width = 40)
Plot the reporter_res as circle_packing
plot_report_circle_packing( reporter_res, rs_threshold = 1.64, mode = 2, facet_anno = NULL, show_ID = FALSE, Pathway_description = TRUE, str_width = 10, show_level_name = "all", show_tip_label = TRUE )plot_report_circle_packing( reporter_res, rs_threshold = 1.64, mode = 2, facet_anno = NULL, show_ID = FALSE, Pathway_description = TRUE, str_width = 10, show_level_name = "all", show_tip_label = TRUE )
reporter_res |
result of 'get_reporter_score' |
rs_threshold |
plot threshold vector, default:1.64 |
mode |
1~2 plot style. |
facet_anno |
annotation table for facet, more two columns, last is pathway name, last second is pathway id. |
show_ID |
show pathway id |
Pathway_description |
show KO description rather than KO id. |
str_width |
str_width to wrap |
show_level_name |
show the level name? |
show_tip_label |
show the tip label? |
ggplot
data("reporter_score_res") if (requireNamespace("igraph") && requireNamespace("ggraph")) { plot_report_circle_packing(reporter_score_res, rs_threshold = c(2, -2), str_width = 40) }data("reporter_score_res") if (requireNamespace("igraph") && requireNamespace("ggraph")) { plot_report_circle_packing(reporter_score_res, rs_threshold = c(2, -2), str_width = 40) }
Plot the significance of pathway
plot_significance(reporter_res, map_id)plot_significance(reporter_res, map_id)
reporter_res |
result of 'get_reporter_score' or 'reporter_score' |
map_id |
the pathway or module id |
ggplot
data("reporter_score_res") plot_significance(reporter_score_res, map_id = c("map05230", "map03010"))data("reporter_score_res") plot_significance(reporter_score_res, map_id = c("map05230", "map03010"))
Plot c_means result
## S3 method for class 'cm_res' plot( x, filter_membership, mode = 1, show.clust.cent = TRUE, show_num = TRUE, ... )## S3 method for class 'cm_res' plot( x, filter_membership, mode = 1, show.clust.cent = TRUE, show_num = TRUE, ... )
x |
a cm_res object |
filter_membership |
filter membership |
mode |
1~2 |
show.clust.cent |
show cluster center? |
show_num |
show number of each cluster? |
... |
additional |
ggplot
Print reporter_score
## S3 method for class 'reporter_score' print(x, ...)## S3 method for class 'reporter_score' print(x, ...)
x |
reporter_score |
... |
add |
No value
Print rs_by_cm
## S3 method for class 'rs_by_cm' print(x, ...)## S3 method for class 'rs_by_cm' print(x, ...)
x |
rs_by_cm |
... |
add |
No value
Transfer p-value of KOs to Z-score
pvalue2zs( ko_pvalue, mode = c("directed", "mixed")[1], p.adjust.method1 = "none" )pvalue2zs( ko_pvalue, mode = c("directed", "mixed")[1], p.adjust.method1 = "none" )
ko_pvalue |
data.frame from |
mode |
'mixed' or 'directed' (default, only for two groups differential analysis or multi-groups correlation analysis.), see details in |
p.adjust.method1 |
p.adjust.method for 'ko.test', see |
'mixed' mode is the original reporter-score method from Patil, K. R. et al. PNAS 2005. In this mode, the reporter score is undirected, and the larger the reporter score, the more significant the enrichment, but it cannot indicate the up-and-down regulation information of the pathway! (Liu, L. et al. iMeta 2023.)
steps:
1. Use the Wilcoxon rank sum test to obtain the P value of the significance of each KO difference between the two groups (ie , i represents a certain KO);
2. Using an inverse normal distribution, convert the P value of each KO into a Z value (), the formula:
3. 'Upgrade' KO to pathway: , calculate the Z value of the pathway, the formula:
where k means A total of k KOs were annotated to the corresponding pathway;
4. Evaluate the degree of significance: permutation (permutation) 1000 times, get the random distribution of , the formula:
is The mean of the random distribution, is the standard deviation of the random distribution.
Instead, 'directed' mode is a derived version of 'mixed', referenced from https://github.com/wangpeng407/ReporterScore.
This approach is based on the same assumption of many differential analysis methods: the expression of most genes has no significant change.
steps:
1. Use the Wilcoxon rank sum test to obtain the P value of the significance of each KO difference between the two groups (ie , i represents a certain KO), and then divide the P value by 2, that is, the range of (0,1] becomes (0,0.5], ;
2. Using an inverse normal distribution, convert the P value of each KO into a Z value (), the formula:
since the above P value is less than 0.5, all Z values will be greater than 0;
3. Considering whether each KO is up-regulated or down-regulated, calculate ,
,
so is greater than 0 Up-regulation, less than 0 is down-regulation;
4. 'Upgrade' KO to pathway: , calculate the Z value of the pathway, the formula:
where k means A total of k KOs were annotated to the corresponding pathway;
5. Evaluate the degree of significance: permutation (permutation) 1000 times, get the random distribution of , the formula:
is The mean of the random distribution, is the standard deviation of the random distribution.
The finally obtained is the Reporter score value enriched for each pathway.
In this mode, the Reporter score is directed, and a larger positive value represents a significant up-regulation enrichment, and a smaller negative values represent significant down-regulation enrichment.
However, the disadvantage of this mode is that when a pathway contains about the same number of significantly up-regulates KOs and significantly down-regulates KOs, the final absolute value of Reporter score may approach 0, becoming a pathway that has not been significantly enriched.
ko_stat data.frame
1. Patil, K. R. & Nielsen, J. Uncovering transcriptional regulation of metabolism by using metabolic network topology. Proc Natl Acad Sci U S A 102, 2685–2689 (2005).
2. Liu, L., Zhu, R. & Wu, D. Misuse of reporter score in microbial enrichment analysis. iMeta n/a, e95.
3. https://github.com/wangpeng407/ReporterScore
Other GRSA:
combine_rs_res(),
get_reporter_score(),
ko.test(),
reporter_score()
data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed")data("KO_abundance_test") ko_pvalue <- ko.test(KO_abundance, "Group", metadata) ko_stat <- pvalue2zs(ko_pvalue, mode = "directed")
One step to get the reporter score of your KO abundance table.
reporter_score( kodf, group, metadata = NULL, method = "wilcox.test", pattern = NULL, p.adjust.method1 = "none", mode = c("directed", "mixed")[1], verbose = TRUE, feature = "ko", type = c("pathway", "module")[1], p.adjust.method2 = "BH", modulelist = NULL, threads = 1, perm = 4999, min_exist_KO = 3, max_exist_KO = 600 )reporter_score( kodf, group, metadata = NULL, method = "wilcox.test", pattern = NULL, p.adjust.method1 = "none", mode = c("directed", "mixed")[1], verbose = TRUE, feature = "ko", type = c("pathway", "module")[1], p.adjust.method2 = "BH", modulelist = NULL, threads = 1, perm = 4999, min_exist_KO = 3, max_exist_KO = 600 )
kodf |
KO_abundance table, rowname are feature ids (e.g. K00001 if feature="ko"; PEX11A if feature="gene"; C00024 if feature="compound"), colnames are samples. |
group |
The comparison groups (at least two categories) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. And you can use factor levels to change order. |
metadata |
sample information data.frame contains group |
method |
the type of test. Default is 'wilcox.test'. Allowed values include:
|
pattern |
a named vector matching the group, e.g. c('G1'=1,'G2'=3,'G3'=2), use the correlation analysis with specific pattern to calculate p-value. |
p.adjust.method1 |
p.adjust.method for 'ko.test', see |
mode |
'mixed' or 'directed' (default, only for two groups differential analysis or multi-groups correlation analysis.), see details in |
verbose |
logical |
feature |
one of 'ko', 'gene', 'compound' |
type |
'pathway' or 'module' for default KOlist for microbiome, 'CC', 'MF', 'BP', 'ALL' for default GOlist for homo sapiens. And org in listed in 'https://www.genome.jp/kegg/catalog/org_list.html' such as 'hsa' (if your kodf is come from a specific organism, you should specify type here). |
p.adjust.method2 |
p.adjust.method for the correction of ReporterScore, see |
modulelist |
NULL or customized modulelist dataframe, must contain 'id','K_num','KOs','Description' columns. Take the 'KOlist' as example, use |
threads |
default 1 |
perm |
permutation number, default: 4999. |
min_exist_KO |
min exist KO number in a pathway (default, 3, when a pathway contains KOs less than 3, there will be no RS) |
max_exist_KO |
max exist KO number in a pathway (default, 600, when a pathway contains KOs more than 600, there will be no RS) |
reporter_score object:
kodf |
your input KO_abundance table |
ko_stat |
ko statistics result contains p.value and z_score |
reporter_s |
the reporter score in each pathway |
modulelist |
default KOlist or customized modulelist dataframe |
group |
The comparison groups in your data |
metadata |
sample information dataframe contains group |
for the 'reporter_s' in result, whose columns represent:
ID |
pathway id |
Description |
pathway description |
K_num |
total number of KOs/genes in the pathway |
Exist_K_num |
number of KOs/genes in your inputdata that exist in the pathway |
Significant_K_num |
number of kos/genes in your inputdata that are significant in the pathway |
Z_score |
|
BG_Mean |
Background mean, |
BG_Sd |
Background standard deviation, |
ReporterScore |
ReporterScore of the pathway, |
p.value |
p.value of the ReporterScore |
p.adjust |
adjusted p.value by p.adjust.method2 |
Other GRSA:
combine_rs_res(),
get_reporter_score(),
ko.test(),
pvalue2zs()
message("The following example require some time to run:") data("KO_abundance_test") reporter_score_res <- reporter_score(KO_abundance, "Group", metadata, mode = "directed", perm = 499 ) head(reporter_score_res$reporter_s) reporter_score_res2 <- reporter_score(KO_abundance, "Group2", metadata, mode = "mixed", method = "kruskal.test", p.adjust.method1 = "none", perm = 499 ) reporter_score_res3 <- reporter_score(KO_abundance, "Group2", metadata, mode = "directed", method = "pearson", pattern = c("G1" = 1, "G2" = 3, "G3" = 2), perm = 499 )message("The following example require some time to run:") data("KO_abundance_test") reporter_score_res <- reporter_score(KO_abundance, "Group", metadata, mode = "directed", perm = 499 ) head(reporter_score_res$reporter_s) reporter_score_res2 <- reporter_score(KO_abundance, "Group2", metadata, mode = "mixed", method = "kruskal.test", p.adjust.method1 = "none", perm = 499 ) reporter_score_res3 <- reporter_score(KO_abundance, "Group2", metadata, mode = "directed", method = "pearson", pattern = c("G1" = 1, "G2" = 3, "G3" = 2), perm = 499 )
'reporter_score()' result from KO_abundance_test
'reporter_score()' result from KO_abundance_test
a list contain 7 elements.
your input KO_abundance table
ko statistics result contains p.value and z_score
the reporter score in each pathway
default KOlist or customized modulelist dataframe
The compare group (two category) in your data
sample information dataframe contains group
Other test_data:
KO_abundance,
genedf
Reporter score analysis after C-means clustering
Extract one cluster from rs_by_cm object
Plot c_means result
RSA_by_cm( kodf, group, metadata = NULL, k_num = NULL, filter_var = 0.7, verbose = TRUE, method = "pearson", ... ) extract_cluster(rsa_cm_res, cluster = 1) plot_c_means( rsa_cm_res, filter_membership, mode = 1, show.clust.cent = TRUE, show_num = TRUE, ... )RSA_by_cm( kodf, group, metadata = NULL, k_num = NULL, filter_var = 0.7, verbose = TRUE, method = "pearson", ... ) extract_cluster(rsa_cm_res, cluster = 1) plot_c_means( rsa_cm_res, filter_membership, mode = 1, show.clust.cent = TRUE, show_num = TRUE, ... )
kodf |
KO_abundance table, rowname is ko id (e.g. K00001),colnames is samples. |
group |
The comparison groups (at least two categories) in your data, one column name of metadata when metadata exist or a vector whose length equal to columns number of kodf. And you can use factor levels to change order. |
metadata |
sample information data.frame contains group |
k_num |
if NULL, perform the cm_test_k, else an integer |
filter_var |
see c_means |
verbose |
verbose |
method |
method from |
... |
additional |
rsa_cm_res |
a cm_res object |
cluster |
integer |
filter_membership |
filter membership 0~1. |
mode |
1~2 |
show.clust.cent |
show cluster center? |
show_num |
show number of each cluster? |
rs_by_cm
reporter_score object
ggplot
Other C_means:
cm_test_k()
message("The following example require some time to run:") if (requireNamespace("e1071") && requireNamespace("factoextra")) { data("KO_abundance_test") rsa_cm_res <- RSA_by_cm(KO_abundance, "Group2", metadata, k_num = 3, filter_var = 0.7, method = "pearson", perm = 199 ) extract_cluster(rsa_cm_res, cluster = 1) }message("The following example require some time to run:") if (requireNamespace("e1071") && requireNamespace("factoextra")) { data("KO_abundance_test") rsa_cm_res <- RSA_by_cm(KO_abundance, "Group2", metadata, k_num = 3, filter_var = 0.7, method = "pearson", perm = 199 ) extract_cluster(rsa_cm_res, cluster = 1) }
Upgrade the KO level
up_level_KO( KO_abundance, level = "pathway", show_name = FALSE, modulelist = NULL, verbose = TRUE )up_level_KO( KO_abundance, level = "pathway", show_name = FALSE, modulelist = NULL, verbose = TRUE )
KO_abundance |
KO_abundance |
level |
one of 'pathway', 'module', 'level1', 'level2', 'level3', 'module1', 'module2', 'module3'. |
show_name |
logical |
modulelist |
NULL or customized modulelist dataframe, must contain 'id','K_num','KOs','Description' columns. Take the 'KOlist' as example, use |
verbose |
logical |
data.frame
data("KO_abundance_test") KO_level1 <- up_level_KO(KO_abundance, level = "level1", show_name = TRUE)data("KO_abundance_test") KO_level1 <- up_level_KO(KO_abundance, level = "level1", show_name = TRUE)
update CARDinfo from (from 'CARD' database)
update_CARDinfo(download_dir = NULL, card_data = NULL)update_CARDinfo(download_dir = NULL, card_data = NULL)
download_dir |
download_dir |
card_data |
card_data from https://card.mcmaster.ca/download/0/broadstreet-v3.2.8.tar.bz2 |
No value
Download links:
http://geneontology.org/docs/download-ontology/
https://asa12138.github.io/FileList/GOlist.rda
update_GOlist(download_dir = NULL, GO_file = NULL) update_GOinfo(download_dir = NULL, obo_file = NULL)update_GOlist(download_dir = NULL, GO_file = NULL) update_GOinfo(download_dir = NULL, obo_file = NULL)
download_dir |
download_dir |
GO_file |
GO_file |
obo_file |
obo_file from http://current.geneontology.org/ontology/go.obo |
No value
Download links:
https://rest.kegg.jp/list/pathway
https://rest.kegg.jp/link/pathway/ko
https://rest.kegg.jp/link/pathway/compound
https://rest.kegg.jp/list/module
https://rest.kegg.jp/link/module/ko
https://rest.kegg.jp/link/module/compound
update_KEGG(download_dir) update_KO_file(download_dir, RDSfile = NULL) update_htable(type, keg_file = NULL, download = FALSE, download_dir = NULL) update_org_pathway( org = "hsa", RDS_file = NULL, download = TRUE, download_dir = NULL )update_KEGG(download_dir) update_KO_file(download_dir, RDSfile = NULL) update_htable(type, keg_file = NULL, download = FALSE, download_dir = NULL) update_org_pathway( org = "hsa", RDS_file = NULL, download = TRUE, download_dir = NULL )
download_dir |
where to save the .keg file? |
RDSfile |
saved KO_files.RDS file |
type |
"ko", "module", "pathway", "compound" ... |
keg_file |
path of a .keg file, such as ko00001.keg from https://www.genome.jp/kegg-bin/download_htext?htext=ko00001&format=htext. |
download |
save the .keg file? |
org |
kegg organism, listed in https://www.genome.jp/kegg/catalog/org_list.html, default, "hsa" |
RDS_file |
path of a org.RDS file if you saved before. |
No value