cafferychen777 / ggpicrust2

Make Picrust2 Output Analysis and Visualization Easier
https://cafferychen777.github.io/ggpicrust2/
MIT License
102 stars 11 forks source link

Error after pathway_errorbar #25

Closed andressamv closed 1 year ago

andressamv commented 1 year ago

Hi (and thanks for creating this tool!)

I am running ggpicrust2 one function at a time but struggling with the pathway_errorbar function.

My code kegg_abundance <- ko2kegg_abundance(abundance_file)

daa_results <- pathway_daa(kegg_abundance, meta_table, 'Day', daa_method = 'ALDEx2', select = NULL, p.adjust = 'bonferroni', reference = NULL)

daa_results_method <- daa_results[daa_results$method == 'ALDEx2_Kruskal-Wallace test', ]

daa_results_method1 <- daa_results_method[1:100,] daa_results_method2 <- daa_results_method[101:200,] daa_results_method3 <- daa_results_method[201:253,]

pa_results1 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method1, ko_to_kegg = TRUE) pa_results2 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method2, ko_to_kegg = TRUE) pa_results3 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method3, ko_to_kegg = TRUE)

pa_results <- rbind(pa_results1, pa_results2, pa_results3)

pa_results_padj <- subset(pa_results, p_adjust < 0.001)

Group <- meta_table$Day

daa_results_list <- pathway_errorbar(abundance = kegg_abundance, daa_results_df = pa_results, Group = Group, ko_to_kegg = TRUE, p_values_threshold = 0.05, order = "pathway_class", select = 'NULL', p_value_bar = TRUE, colors = NULL, x_lab = 'description')

_Output of pathwayerrorbar The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

Can you please help me with this error? Everything else worked fine! My pa_results table doesn't have a column called x_lab. Is this the problem?

cafferychen777 commented 1 year ago

Hi,

Thanks for using ggpicrust2 and for reaching out. Regarding the error message you received with the pathway_errorbar function, it seems that the error message "The feature with statistically significance is zero, pathway_errorbar can't do the visualization." means that there are no results in pa_results with p.adjust values less than 0.05. As a result, the function is unable to perform the visualization.

To resolve this issue, I suggest checking the pa_results dataframe to ensure that it contains the expected results with p.adjust values less than 0.05. If the results are not as expected, you may want to adjust the analysis parameters or thresholds to obtain the desired results.

Please let me know if you have any further questions or concerns.

Best, Chen YANG

andressamv @.***> 于2023年5月10日周三 06:37写道:

Hi (and thanks for creating this tool!)

I am running ggpicrust2 one function at a time but struggling with the pathway_errorbar function.

My code kegg_abundance <- ko2kegg_abundance(abundance_file)

daa_results <- pathway_daa(kegg_abundance, meta_table, 'Day', daa_method = 'ALDEx2', select = NULL, p.adjust = 'bonferroni', reference = NULL)

daa_results_method <- daa_results[daa_results$method == 'ALDEx2_Kruskal-Wallace test', ]

daa_results_method1 <- daa_results_method[1:100,] daa_results_method2 <- daa_results_method[101:200,] daa_results_method3 <- daa_results_method[201:253,]

pa_results1 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method1, ko_to_kegg = TRUE) pa_results2 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method2, ko_to_kegg = TRUE) pa_results3 <- pathway_annotation(pathway = 'KO', daa_results_df = daa_results_method3, ko_to_kegg = TRUE)

pa_results <- rbind(pa_results1, pa_results2, pa_results3)

pa_results_padj <- subset(pa_results, p_adjust < 0.001)

Group <- meta_table$Day

daa_results_list <- pathway_errorbar(abundance = kegg_abundance, daa_results_df = pa_results, Group = Group, ko_to_kegg = TRUE, p_values_threshold = 0.05, order = "pathway_class", select = 'NULL', p_value_bar = TRUE, colors = NULL, x_lab = 'description')

Output of pathway_errorbar The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

Can you please help me with this error? Everything else worked fine!

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTVV4IATRWDQXSERS4LXFLBLJANCNFSM6AAAAAAX35JFD4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

andressamv commented 1 year ago

Thank you for the quick response! In the pa_results, most p.adjust values are less than 0.05 (the lowest is 7.598821e-05).

I actually copy the error with one line missing. Please find the correct version below: There is no x_lab you defined in daa_results_df The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

cafferychen777 commented 1 year ago

I guess I know the problem. Please set x_lab = 'pathway_name' because there is no "description" column in pa_results.

andressamv @.***> 于2023年5月10日周三 12:59写道:

Thank you for the quick response! In the pa_results, most p.adjust values are less than 0.05 (the lowest is 7.598821e-05).

I actually copy the error with one line missing. Please find the correct version below: There is no x_lab you defined in daa_results_df The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25#issuecomment-1541353494, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTXJQHHT2KWFNN6QUJLXFMOCDANCNFSM6AAAAAAX35JFD4 . You are receiving this because you commented.Message ID: @.***>

andressamv commented 1 year ago

I changed and now I get the following error: The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

cafferychen777 commented 1 year ago

It's very strange if you do have the feature whose p value < 0.05. Would you like to share the data?

[image: Mailtrack] https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& Sender notified by Mailtrack https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& 23/05/10 下午01:46:42

andressamv @.***> 于2023年5月10日周三 13:41写道:

I changed and now I get the following error: The feature with statistically significance is zero, pathway_errorbar can't do the visualization. Error in h(simpleError(msg, call)) : error in evaluating the argument 'physeq' in selecting a method for function 'taxa_are_rows': invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25#issuecomment-1541381299, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTX34QYECUTYULO74YDXFMTATANCNFSM6AAAAAAX35JFD4 . You are receiving this because you commented.Message ID: @.***>

andressamv commented 1 year ago

Sure! kegg_abundance.txt pa_results.txt

cafferychen777 commented 1 year ago

Metadata is also needed. 😂

[image: Mailtrack] https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& Sender notified by Mailtrack https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& 23/05/10 下午01:56:02

andressamv @.***> 于2023年5月10日周三 13:55写道:

Sure! kegg_abundance.txt https://github.com/cafferychen777/ggpicrust2/files/11438523/kegg_abundance.txt pa_results.txt https://github.com/cafferychen777/ggpicrust2/files/11438525/pa_results.txt

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25#issuecomment-1541393226, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTSJOOPF7JGOFTB634TXFMUUNANCNFSM6AAAAAAX35JFD4 . You are receiving this because you commented.Message ID: @.***>

andressamv commented 1 year ago

sorry, it is quite late here :sweat_smile:

meta_table.txt

cafferychen777 commented 1 year ago
image
library(readr)
library(ggpicrust2)
library(tidyverse)
library(patchwork)
library(ggprism)

daa_results_df_annotated <-
  read.delim(
    "/Users/apple/Microbiome/ggpicrust2总/ggpicrust2测试/ggpicrust2_test/andressamv/pa_results.txt"
  ) 

kegg_abundance <-
  read.delim(
    "/Users/apple/Microbiome/ggpicrust2总/ggpicrust2测试/ggpicrust2_test/andressamv/kegg_abundance.txt"
  ) %>% column_to_rownames("X")

metadata <-
  read_delim(
    "/Users/apple/Microbiome/ggpicrust2总/ggpicrust2测试/ggpicrust2_test/andressamv/meta_table.txt",
    delim = "\t",
    escape_double = FALSE,
    trim_ws = TRUE
  )

daa_results_df_annotated$p_adjust <- round(daa_results_df_annotated$p_adjust,5)

p <- ggpicrust2::pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_results_df_annotated,
  Group = metadata$Day,
  p_values_threshold = 0.05,
  order = "pathway_class",
  select = c("ko05340", "ko00564", "ko00680", "ko00562", "ko03030", "ko00561", "ko00440", "ko00250", "ko00740", "ko04940", "ko00010", "ko00195", "ko00760", "ko00920", "ko00311", "ko00310", "ko04146", "ko00600", "ko04141", "ko04142", "ko00604", "ko04260", "ko00909", "ko04973", "ko00510", "ko04974"),
  ko_to_kegg = TRUE,
  p_value_bar = FALSE,
  colors = NULL,
  x_lab = "pathway_name"
)

It works well in my m1 macbook. I don't know what's going on. Maybe you can check the R session. The problem appearing in my mac is following. I use select parameter to solve it.

The feature with statistically significance are more than 30, the visualization will be terrible.
 Please use select to reduce the number.
 Now you have "ko05340", "ko00564", "ko00680", "ko00562", "ko03030", "ko00561", "ko00440", "ko00250", "ko00740", "ko04940", "ko00010", "ko00195", "ko00760", "ko00920", "ko00311", "ko00310", "ko04146", "ko00600", "ko04141", "ko04142", "ko00604", "ko04260", "ko00909", "ko04973", "ko00510", "ko04974", "ko00450", "ko01051", "ko00565", "ko00524", "ko00300", "ko00260", "ko00190", "ko03440", "ko00750", "ko00950", "ko00591", "ko00590", "ko00061", "ko03070", "ko00253", "ko03060", "ko04740", "ko00380", "ko00500", "ko05120", "ko05322", "ko04964", "ko04962", "ko00625", "ko00627", "ko00626", "ko00621", "ko00620", "ko00623", "ko00622", "ko00270", "ko00940", "ko01053", "ko01056", "ko00071", "ko00072", "ko05219", "ko05215", "ko05214", "ko05211", "ko01055", "ko04910", "ko00531"

Here is my R session info.

R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Ventura 13.3

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] ggprism_1.0.4    patchwork_1.1.2  lubridate_1.9.2  forcats_1.0.0    stringr_1.5.0    dplyr_1.1.1     
 [7] purrr_1.0.1      tidyr_1.3.0      tibble_3.2.1     ggplot2_3.4.2    tidyverse_2.0.0  ggpicrust2_1.6.3
[13] readr_2.1.4     

loaded via a namespace (and not attached):
  [1] nlme_3.1-162           bitops_1.0-7           fs_1.6.1               usethis_2.1.6         
  [5] phyloseq_1.42.0        bit64_4.0.5            devtools_2.4.5         RColorBrewer_1.1-3    
  [9] GenomeInfoDb_1.34.9    tools_4.2.2            profvis_0.3.7          utf8_1.2.3            
 [13] R6_2.5.1               vegan_2.6-4            BiocGenerics_0.44.0    mgcv_1.8-42           
 [17] colorspace_2.1-0       permute_0.9-7          rhdf5filters_1.10.1    ade4_1.7-22           
 [21] withr_2.5.0            urlchecker_1.0.1       tidyselect_1.2.0       prettyunits_1.1.1     
 [25] GGally_2.1.2           processx_3.8.0         bit_4.0.5              compiler_4.2.2        
 [29] cli_3.6.1              Biobase_2.58.0         labeling_0.4.2         scales_1.2.1          
 [33] callr_3.7.3            digest_0.6.31          XVector_0.38.0         pkgconfig_2.0.3       
 [37] htmltools_0.5.5        sessioninfo_1.2.2      fastmap_1.1.1          htmlwidgets_1.6.2     
 [41] rlang_1.1.0            rstudioapi_0.14        shiny_1.7.4            farver_2.1.1          
 [45] generics_0.1.3         jsonlite_1.8.4         vroom_1.6.1            RCurl_1.98-1.12       
 [49] magrittr_2.0.3         GenomeInfoDbData_1.2.9 biomformat_1.26.0      Matrix_1.5-4          
 [53] Rcpp_1.0.10            munsell_0.5.0          S4Vectors_0.36.2       Rhdf5lib_1.20.0       
 [57] fansi_1.0.4            ape_5.7-1              lifecycle_1.0.3        stringi_1.7.12        
 [61] MASS_7.3-59            zlibbioc_1.44.0        rhdf5_2.42.0           pkgbuild_1.4.0        
 [65] plyr_1.8.8             grid_4.2.2             parallel_4.2.2         promises_1.2.0.1      
 [69] crayon_1.5.2           miniUI_0.1.1.1         lattice_0.21-8         Biostrings_2.66.0     
 [73] splines_4.2.2          multtest_2.54.0        hms_1.1.3              ps_1.7.4              
 [77] pillar_1.9.0           igraph_1.4.1           reshape2_1.4.4         codetools_0.2-19      
 [81] stats4_4.2.2           pkgload_1.3.2          glue_1.6.2             data.table_1.14.8     
 [85] remotes_2.4.2          renv_0.16.0            tzdb_0.3.0             vctrs_0.6.1           
 [89] httpuv_1.6.9           foreach_1.5.2          gtable_0.3.3           reshape_0.8.9         
 [93] cachem_1.0.7           mime_0.12              xtable_1.8-4           later_1.3.0           
 [97] survival_3.5-5         iterators_1.0.14       memoise_2.0.1          IRanges_2.32.0        
[101] cluster_2.1.4          timechange_0.2.0       ellipsis_0.3.2       
andressamv commented 1 year ago

Now it worked! The last two packages were not loaded (patchwork and ggprism). Thank you so much!

I tried using the code Top 30 for selection but it is not working. How can I select the top 30 manually? Is this base on abundance or p values?

cafferychen777 commented 1 year ago
daa_results_df_annotated <- daa_results_df_annotated[!is.na(daa_results_df_annotated$pathway_name),]

daa_results_df_annotated$p_adjust <- round(daa_results_df_annotated$p_adjust,5)

low_p_feature <- daa_results_df_annotated[order(daa_results_df_annotated$p_adjust), ]$feature[1:20]

p <- ggpicrust2::pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_results_df_annotated,
  Group = metadata$Day,
  p_values_threshold = 0.05,
  order = "pathway_class",
  select = low_p_feature,
  ko_to_kegg = TRUE,
  p_value_bar = FALSE,
  colors = NULL,
  x_lab = "pathway_name"
)

I think it can solve the problem.

andressamv commented 1 year ago

That works perfectly! Thank you so much!

I am sorry for all the questions, but I have one more. I am trying to change the order of the factor "Day" (i.e., Day1, Day2, Day4, etc) for cosmetic purposes. I normally do this before plotting using ggplot. However, my way is not really working using the package. What do you recommend for it?

I am also using and getting the following error:

kegg_abundance_t <- as.data.frame(t(kegg_abundance)) pca_plot <- pathway_pca(t(kegg_abundance_t), meta_table, 'Day'); print(pca_plot)

Error in pca[, group] : subscript out of bounds

cafferychen777 commented 1 year ago
截屏2023-05-11 13 59 05 截屏2023-05-11 13 59 45

I just updated ggpicrust2 to version 1.6.4. You can use the following codes after update the version through Github.

Group <- factor(metadata$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36"))

p <- pathway_errorbar(
  abundance = kegg_abundance,
  daa_results_df = daa_results_df_annotated,
  Group = Group,
  p_values_threshold = 0.05,
  order = "pathway_class",
  select = low_p_feature,
  ko_to_kegg = TRUE,
  p_value_bar = FALSE,
  colors = NULL,
  x_lab = "pathway_name"
)

metadata$Day <- factor(metadata$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36"))
pathway_pca(kegg_abundance, metadata, 'Day')

Best regards,

andressamv commented 1 year ago

Thank you so much for your help! I can not describe how thankful I am. This tool is definitely a lifesaver!

I am still finding some errors after running this script. Do you think is related to the dependencies in my computer?

Group <- factor(meta_table$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36")) pathway_errorbar(abundance = kegg_abundance, daa_results_df = annotation_results, Group = Group, p_values_threshold = 0.05, order = 'pathway_class', ko_to_kegg = TRUE, select = low_p_feature, p_value_bar = FALSE, colors = NULL, x_lab = 'pathway_name')

Error in factor(error_bar_df$group, levels = levels(as.factor(Group))) : object 'error_bar_df' not found

meta_table$Day <- factor(meta_table$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36")) pathway_pca(kegg_abundance, meta_table, 'Day')

Error in pca[, group] : subscript out of bounds

cafferychen777 commented 1 year ago

I am sorry about the pathway_errorbar error. I upload the wrong version. Please re-install ggpicrust2 with Github. Sorry again 😭. About the pathway_pca. I have no idea about it. Please you can check the session info. Following is my sessioninfor.

R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Ventura 13.3

Matrix products: default
LAPACK:
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base

other attached packages:
 [1] ggprism_1.0.4    patchwork_1.1.2  lubridate_1.9.2  forcats_1.0.0
 stringr_1.5.0    dplyr_1.1.2
 [7] purrr_1.0.1      tidyr_1.3.0      tibble_3.2.1     ggplot2_3.4.2
 tidyverse_2.0.0  readr_2.1.4
[13] ggpicrust2_1.6.4

loaded via a namespace (and not attached):
  [1] colorspace_2.1-0       ellipsis_0.3.2         rprojroot_2.0.3
 XVector_0.38.0
  [5] fs_1.6.2               aplot_0.1.10           rstudioapi_0.14
 farver_2.1.1
  [9] remotes_2.4.2          bit64_4.0.5            fansi_1.0.4
 codetools_0.2-19
 [13] splines_4.2.2          cachem_1.0.7           pkgload_1.3.2
 ade4_1.7-22
 [17] jsonlite_1.8.4         phyloseq_1.42.0        cluster_2.1.4
 shiny_1.7.4
 [21] compiler_4.2.2         Matrix_1.5-4           fastmap_1.1.1
 cli_3.6.1
 [25] later_1.3.0            htmltools_0.5.5        prettyunits_1.1.1
 tools_4.2.2
 [29] igraph_1.4.1           gtable_0.3.3           glue_1.6.2
GenomeInfoDbData_1.2.9
 [33] reshape2_1.4.4         Rcpp_1.0.10            Biobase_2.58.0
vctrs_0.6.2
 [37] Biostrings_2.66.0      rhdf5filters_1.10.1    multtest_2.54.0
 ape_5.7-1
 [41] nlme_3.1-162           iterators_1.0.14       ps_1.7.5
timechange_0.2.0
 [45] mime_0.12              miniUI_0.1.1.1         lifecycle_1.0.3
 renv_0.16.0
 [49] devtools_2.4.5         zlibbioc_1.44.0        MASS_7.3-59
 scales_1.2.1
 [53] vroom_1.6.3            hms_1.1.3              promises_1.2.0.1
parallel_4.2.2
 [57] biomformat_1.26.0      rhdf5_2.42.0           RColorBrewer_1.1-3
curl_5.0.0
 [61] memoise_2.0.1          ggfun_0.0.9            yulab.utils_0.0.6
 reshape_0.8.9
 [65] stringi_1.7.12         S4Vectors_0.36.2       desc_1.4.2
foreach_1.5.2
 [69] permute_0.9-7          BiocGenerics_0.44.0    pkgbuild_1.4.0
GenomeInfoDb_1.34.9
 [73] rlang_1.1.1            pkgconfig_2.0.3        bitops_1.0-7
lattice_0.21-8
 [77] Rhdf5lib_1.20.0        labeling_0.4.2         htmlwidgets_1.6.2
 bit_4.0.5
 [81] processx_3.8.1         tidyselect_1.2.0       GGally_2.1.2
plyr_1.8.8
 [85] magrittr_2.0.3         R6_2.5.1               IRanges_2.32.0
generics_0.1.3
 [89] profvis_0.3.7          pillar_1.9.0           withr_2.5.0
 mgcv_1.8-42
 [93] survival_3.5-5         RCurl_1.98-1.12        crayon_1.5.2
utf8_1.2.3
 [97] tzdb_0.3.0             urlchecker_1.0.1       usethis_2.1.6
 grid_4.2.2
[101] data.table_1.14.8      callr_3.7.3            vegan_2.6-4
 digest_0.6.31
[105] xtable_1.8-4           httpuv_1.6.9           gridGraphics_0.5-1
stats4_4.2.2
[109] munsell_0.5.0          ggplotify_0.1.0        sessioninfo_1.2.2

[image: Mailtrack] https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& Sender notified by Mailtrack https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& 23/05/12 上午12:42:44

andressamv @.***> 于2023年5月12日周五 00:24写道:

Thank you so much for your help! This tool is definitely a lifesaver!

I am still finding some errors after running this script. Do you think is related to the dependencies in my computer?

Group <- factor(meta_table$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36")) pathway_errorbar(abundance = kegg_abundance, daa_results_df = annotation_results, Group = Group,

-

                 p_values_threshold = 0.05, order = 'pathway_class', ko_to_kegg = TRUE,

-

                 select = low_p_feature, p_value_bar = FALSE, colors = NULL, x_lab = 'pathway_name')

Error in factor(error_bar_df$group, levels = levels(as.factor(Group))) : object 'error_bar_df' not found

meta_table$Day <- factor(meta_table$Day, levels = c("Day0","Day1","Day2","Day4","Day8","Day16","Day24","Day36")) pathway_pca(kegg_abundance, meta_table, 'Day')

Error in pca[, group] : subscript out of bounds

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25#issuecomment-1544298890, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTQTZA4ZHYMOPOLHCW3XFUHD5ANCNFSM6AAAAAAX35JFD4 . You are receiving this because you commented.Message ID: @.***>

andressamv commented 1 year ago

Ohh, thanks! However, I am still getting the same errors. I am copying here the info I receive after installing ggpicrust2.

Downloading GitHub repo cafferychen777/ggpicrust2@HEAD Skipping 8 packages not available: SummarizedExperiment, metagenomeSeq, Maaslin2, limma, lefser, edgeR, DESeq2, ALDEx2 ── R CMD build ─────────────────────────────────────────────────────────────────────────────────────────── ✔ checking for file ‘/private/var/folders/bj/9nzmk8w529gbq758w4rm5b6w0000gn/T/RtmpSedrhY/remotes117815c018d32/cafferychen777-ggpicrust2-b676db5/DESCRIPTION’ (396ms) ─ preparing ‘ggpicrust2’: ✔ checking DESCRIPTION meta-information ... ─ checking for LF line-endings in source and make files and shell scripts ─ checking for empty or unneeded directories Removed empty directory ‘ggpicrust2/data’ Omitted ‘LazyData’ from DESCRIPTION NB: this package now depends on R (>= 3.5.0) WARNING: Added dependency on R >= 3.5.0 because serialized objects in serialize/load version 3 cannot be read in older versions of R. File(s) containing such objects: ‘ggpicrust2/inst/extdata/EC_reference.RData’ ‘ggpicrust2/inst/extdata/KO_reference.RData’ ‘ggpicrust2/inst/extdata/MetaCyc_reference.RData’ ‘ggpicrust2/inst/extdata/kegg_reference.RData’ ─ building ‘ggpicrust2_1.6.4.tar.gz’

cafferychen777 commented 1 year ago

It's very strange. Because I update the newest version to Github, but Github doesn't update the content which we downloaded. It may need some time to let github update the content.

[image: Mailtrack] https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& Sender notified by Mailtrack https://mailtrack.io?utm_source=gmail&utm_medium=signature&utm_campaign=signaturevirality11& 23/05/12 上午01:06:32

andressamv @.***> 于2023年5月12日周五 00:51写道:

Ohh, thanks! However, I am still getting the same error. I am copying here the info I receive after installing ggpicrust2.

Downloading GitHub repo @.*** Skipping 8 packages not available: SummarizedExperiment, metagenomeSeq, Maaslin2, limma, lefser, edgeR, DESeq2, ALDEx2 ── R CMD build ─────────────────────────────────────────────────────────────────────────────────────────── ✔ checking for file ‘/private/var/folders/bj/9nzmk8w529gbq758w4rm5b6w0000gn/T/RtmpSedrhY/remotes117815c018d32/cafferychen777-ggpicrust2-b676db5/DESCRIPTION’ (396ms) ─ preparing ‘ggpicrust2’: ✔ checking DESCRIPTION meta-information ... ─ checking for LF line-endings in source and make files and shell scripts ─ checking for empty or unneeded directories Removed empty directory ‘ggpicrust2/data’ Omitted ‘LazyData’ from DESCRIPTION NB: this package now depends on R (>= 3.5.0) WARNING: Added dependency on R >= 3.5.0 because serialized objects in serialize/load version 3 cannot be read in older versions of R. File(s) containing such objects: ‘ggpicrust2/inst/extdata/EC_reference.RData’ ‘ggpicrust2/inst/extdata/KO_reference.RData’ ‘ggpicrust2/inst/extdata/MetaCyc_reference.RData’ ‘ggpicrust2/inst/extdata/kegg_reference.RData’ ─ building ‘ggpicrust2_1.6.4.tar.gz’

  • installing source package ‘ggpicrust2’ ... using staged installation R inst byte-compile and prepare package for lazy loading ** help * installing help indices building package indices testing if installed package can be loaded from temporary location testing if installed package can be loaded from final location ** testing if installed package keeps a record of temporary installation path
  • DONE (ggpicrust2) Warning message: In readLines(old_path) : incomplete final line found on '/Users/andressa/.R/Makevars'

— Reply to this email directly, view it on GitHub https://github.com/cafferychen777/ggpicrust2/issues/25#issuecomment-1544336581, or unsubscribe https://github.com/notifications/unsubscribe-auth/ATZEQTSDXGD7XHVPLS2AXU3XFUKI7ANCNFSM6AAAAAAX35JFD4 . You are receiving this because you commented.Message ID: @.***>

cafferychen777 commented 1 year ago
image

I found the problem and I fixed it. Every time I can't wait to add new features and fix code, something always goes wrong. I am also troubled.😂 Maybe you can reinstall it after 2 hours or just the following function. Sorry again.

pathway_errorbar <-
  function(abundance,
           daa_results_df,
           Group,
           ko_to_kegg = FALSE,
           p_values_threshold = 0.05,
           order = "group",
           select = NULL,
           p_value_bar = TRUE,
           colors = NULL,
           x_lab = NULL) {
    if (is.null(x_lab)){
      if (ko_to_kegg == TRUE){
        x_lab <- "pathway_name"
      }else{
        x_lab <- "description"
      }
      if (is.null(daa_results_df$pathway_name)&is.null(daa_results_df$description)){
        message("Please use pathway_annotation to annotate the daa_results_df")
      }
    }
    if (!(x_lab %in% colnames(daa_results_df))){
      message("There is no x_lab you defined in daa_results_df")
    }
    if (nlevels(factor(daa_results_df$method)) != 1) {
      message("There are more than one method in daa_results_df$method, please filter it.")
    }
    if (is.null(colors)) {
      colors <- c("#d93c3e", "#3685bc", "#6faa3e", "#e8a825", "#c973e6", "#ee6b3d", "#2db0a7", "#f25292")[1:nlevels(as.factor(Group))]
    }
    errorbar_abundance_mat <- as.matrix(abundance)
    daa_results_filtered_df <-
      daa_results_df[daa_results_df$p_adjust < p_values_threshold,]
    if (!is.null(select)) {
      daa_results_filtered_sub_df <-
        daa_results_filtered_df[daa_results_filtered_df$feature %in% select, ]
    } else {
      daa_results_filtered_sub_df <- daa_results_filtered_df
    }
    if (nrow(daa_results_filtered_sub_df) > 30) {
      stop(
        paste0(
          "The feature with statistically significance are more than 30, the visualization will be terrible.\n Please use select to reduce the number.\n Now you have ",
          paste(paste0('"', daa_results_filtered_sub_df$feature, '"'), collapse = ", ")
        )
      )
    }
    if (nrow(daa_results_filtered_sub_df) == 0){
      message("The feature with statistically significance is zero, pathway_errorbar can't do the visualization.")
    }
    errorbar_sub_abundance_mat <-
      errorbar_abundance_mat[rownames(errorbar_abundance_mat) %in% daa_results_filtered_sub_df$feature,]
    errorbar_sub_relative_abundance_mat <-
      as.matrix(as.data.frame(phyloseq::transform_sample_counts(phyloseq::otu_table(errorbar_sub_abundance_mat,taxa_are_rows = TRUE), function(x) x/sum(x))))
    error_bar_matrix <-
      cbind(
        sample = colnames(errorbar_sub_relative_abundance_mat),
        group = Group,
        t(errorbar_sub_relative_abundance_mat)
      )
    error_bar_df <- as.data.frame(error_bar_matrix)
    error_bar_df$group <- factor(Group,levels = levels(as.factor(Group)))
      error_bar_pivot_longer_df <- tidyr::pivot_longer(error_bar_df,-c(sample, group))
    error_bar_pivot_longer_tibble <-
      mutate(error_bar_pivot_longer_df, group = as.factor(group))
    error_bar_pivot_longer_tibble$sample <-
      factor(error_bar_pivot_longer_tibble$sample)
    error_bar_pivot_longer_tibble$name <-
      factor(error_bar_pivot_longer_tibble$name)
    error_bar_pivot_longer_tibble$value <-
      as.numeric(error_bar_pivot_longer_tibble$value)
    error_bar_pivot_longer_tibble_summarised <-
      error_bar_pivot_longer_tibble %>% group_by(name, group) %>%
      summarise(mean = mean(value), sd = stats::sd(value))
    error_bar_pivot_longer_tibble_summarised <-
      error_bar_pivot_longer_tibble_summarised %>% mutate(group2 = "nonsense")
    switch(
      order,
      "p_values" = {
        order <- order(daa_results_filtered_sub_df$p_adjust)
      },
      "name" = {
        order <- order(daa_results_filtered_sub_df$feature)
      },
      "group" = {
        daa_results_filtered_sub_df$pro <- 1
        for (i in levels(error_bar_pivot_longer_tibble_summarised$name)) {
          error_bar_pivot_longer_tibble_summarised_sub <-
            error_bar_pivot_longer_tibble_summarised[error_bar_pivot_longer_tibble_summarised$name ==
                                                       i,]
          pro_group <-
            error_bar_pivot_longer_tibble_summarised_sub[error_bar_pivot_longer_tibble_summarised_sub$mean ==
                                                           max(error_bar_pivot_longer_tibble_summarised_sub$mean),]$group
          pro_group <- as.vector(pro_group)
          daa_results_filtered_sub_df[daa_results_filtered_sub_df$feature ==
                                        i,]$pro <- pro_group
        }
        order <-
          order(daa_results_filtered_sub_df$pro,
                daa_results_filtered_sub_df$p_adjust)
      },
      "pathway_class" = {
        if (!"pathway_class" %in% colnames(daa_results_filtered_sub_df)) {
          stop("Please use pathway_annotation function to annotate the pathway_daa results")
        }
        order <- order(
          daa_results_filtered_sub_df$pathway_class,
          daa_results_filtered_sub_df$p_adjust
        )
      },
      {
        order <- order
      }
    )
    daa_results_filtered_sub_df <-
      daa_results_filtered_sub_df[order,]
    error_bar_pivot_longer_tibble_summarised_ordered <-
      data.frame(
        name = NULL,
        group = NULL,
        mean = NULL,
        sd = NULL
      )
    for (i in daa_results_filtered_sub_df$feature) {
      error_bar_pivot_longer_tibble_summarised_ordered <-
        rbind(
          error_bar_pivot_longer_tibble_summarised_ordered,
          error_bar_pivot_longer_tibble_summarised[error_bar_pivot_longer_tibble_summarised$name ==
                                                     i,]
        )
    }
    if (ko_to_kegg == FALSE){
      error_bar_pivot_longer_tibble_summarised_ordered[, x_lab] <-
        rep(daa_results_filtered_sub_df[, x_lab], each = length(levels(
          factor(error_bar_pivot_longer_tibble_summarised_ordered$group)
        )))
    }

    # levels(error_bar_pivot_longer_tibble_summarised_ordered$name) <-
    #   rev(daa_results_filtered_sub_df$feature)
    #
    if (ko_to_kegg == TRUE) {
      error_bar_pivot_longer_tibble_summarised_ordered$pathway_class <-
        rep(daa_results_filtered_sub_df$pathway_class,
            each = length(levels(
              factor(error_bar_pivot_longer_tibble_summarised_ordered$group)
            )))
    }
    error_bar_pivot_longer_tibble_summarised_ordered$name <- factor(error_bar_pivot_longer_tibble_summarised_ordered$name, levels = rev(daa_results_filtered_sub_df$feature))

    #error_bar_pivot_longer_tibble_summarised_ordered$order <- rep(0:(nrow(daa_results_filtered_sub_df)-1),each=2)
    bar_errorbar <-
      ggplot2::ggplot(error_bar_pivot_longer_tibble_summarised_ordered, # nolint: object_usage_linter.
             ggplot2::aes(mean, name, fill = group)) + # nolint
      ggplot2::geom_errorbar(
        ggplot2::aes(xmax = mean + sd, xmin = 0),
        position = ggplot2::position_dodge(width = 0.8),
        width = 0.5,
        size = 0.5,
        color = "black"
      ) +
      ggplot2::geom_bar(stat = "identity",
               position = ggplot2::position_dodge(width = 0.8),
               width = 0.8) +
      GGally::geom_stripped_cols() +
      ggplot2::scale_fill_manual(values = colors) +
      ggplot2::scale_color_manual(values = colors) +
      ggprism::theme_prism() +
      ggplot2::scale_x_continuous(expand = c(0, 0),
                         guide = "prism_offset_minor",) +
      ggplot2::scale_y_discrete(labels = rev(daa_results_filtered_sub_df[, x_lab])) +
      ggplot2::labs(x = "Relative Abundance(%)", y = NULL) +
      ggplot2::theme(
        axis.ticks.y = ggplot2::element_blank(),
        axis.line.y = ggplot2::element_blank(),
        axis.line.x = ggplot2::element_line(size = 0.5),
        axis.ticks.x = ggplot2::element_line(size = 0.5),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        axis.text = ggplot2::element_text(size = 10, color = "black"), # nolint
        axis.text.x = ggplot2::element_text(margin = ggplot2::margin(r = 0)), # nolint
        axis.text.y = ggplot2::element_text(
          size = 10,
          color = "black",
          margin = ggplot2::margin(b = 6)
        ),
        axis.title.x =  ggplot2::element_text(
          size = 10,
          color = "black",
          hjust = 0.5
        ),
        legend.position = "top",
        legend.key.size = ggplot2::unit(0.1, "cm"),
        legend.direction = "vertical",
        legend.justification = "left",
        legend.text = ggplot2::element_text(size = 8, face = "bold"),
        legend.box.just = "right",
        plot.margin = ggplot2::margin(0, 0.5, 0.5, 0, unit = "cm")
      ) + ggplot2::coord_cartesian(clip = "off")

    if (ko_to_kegg == TRUE) {
      pathway_class_group_mat <-
        daa_results_filtered_sub_df$pathway_class %>%
        table() %>%
        data.frame() %>% column_to_rownames(".")
      pathway_class_group <- data.frame(.= unique(daa_results_filtered_sub_df$pathway_class),Freq = pathway_class_group_mat[unique(daa_results_filtered_sub_df$pathway_class),])
      start <-
        c(1, rev(pathway_class_group$Freq)[1:(length(pathway_class_group$Freq) - 1)]) %>%
        cumsum()
      end <- cumsum(rev(pathway_class_group$Freq))
      ymin <- start - 1 / 2
      ymax <- end + 1 / 2
      nPoints <- length(start)
      pCol <- c("#D51F26","#272E6A","#208A42","#89288F","#F47D2B",
                 "#FEE500","#8A9FD1","#C06CAB","#E6C2DC","#90D5E4",
                 "#89C75F","#F37B7D","#9983BD","#D24B27","#3BBCA8",
                 "#6E4B9E","#0C727C", "#7E1416","#D8A767","#3D3D3D")[1:nPoints]
      pFill <- pCol
      for (i in 1:nPoints)  {
        bar_errorbar <- bar_errorbar +
          ggplot2::annotation_custom(
            grob = grid::rectGrob(
              gp = grid::gpar(
                col = pCol[i],
                fill = pFill[i],
                lty = NULL,
                lwd = NULL,
                alpha = 0.2
              )
            ),
            xmin = ggplot2::unit(-2, 'native'),
            xmax = ggplot2::unit(0, 'native'),
            ymin = ggplot2::unit(ymin[i], 'native'),
            ymax = ggplot2::unit(ymax[i], 'native')
          )
      }
    }
    daa_results_filtered_sub_df <-
      cbind(
        daa_results_filtered_sub_df,
        negative_log10_p = -log10(daa_results_filtered_sub_df$p_adjust),
        group_nonsense = "nonsense",
        log_2_fold_change = NA
      )

    for (i in daa_results_filtered_sub_df$feature){
      mean <- error_bar_pivot_longer_tibble_summarised_ordered[error_bar_pivot_longer_tibble_summarised_ordered$name %in% i,]$mean
      daa_results_filtered_sub_df[daa_results_filtered_sub_df$feature==i,]$log_2_fold_change <- log2(mean[1]/mean[2])
    }
    daa_results_filtered_sub_df$feature <- factor(daa_results_filtered_sub_df$feature,levels = rev(daa_results_filtered_sub_df$feature))
    p_values_bar <- daa_results_filtered_sub_df %>%
      ggplot2::ggplot(ggplot2::aes(feature, log_2_fold_change, fill = group_nonsense)) +
      ggplot2::geom_bar(stat = "identity",
               position = ggplot2::position_dodge(width = 0.8),
               width = 0.8) +
      ggplot2::labs(y = "log2 fold change", x = NULL) +
      GGally::geom_stripped_cols() +
      ggplot2::scale_fill_manual(values = "#87ceeb") +
      ggplot2::scale_color_manual(values = "#87ceeb") +
      ggplot2::geom_hline(ggplot2::aes(yintercept = 0),
                 linetype = 'dashed',
                 color = 'black') +
      ggprism::theme_prism() +
      ggplot2::scale_y_continuous(expand = c(0, 0),
                         guide = "prism_offset_minor") +
      ggplot2::theme(
        axis.ticks.y = ggplot2::element_blank(),
        axis.line.y = ggplot2::element_blank(),
        axis.line.x = ggplot2::element_line(size = 0.5),
        axis.ticks.x = ggplot2::element_line(size = 0.5),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        axis.text = ggplot2::element_text(size = 10, color = "black"),
        axis.text.y = ggplot2::element_blank(),
        axis.text.x = ggplot2::element_text(
          size = 10,
          color = "black",
          margin = ggplot2::margin(b = 6)
        ),
        axis.title.x =  ggplot2::element_text(
          size = 11,
          color = "black",
          hjust = 0.5
        ),
        legend.position = "non"
      ) +
      ggplot2::coord_flip()

    if (ko_to_kegg == TRUE) {
      pathway_class_y <- (ymax + ymin) / 2 - 0.5
      pathway_class_plot_df <-
        as.data.frame(
          cbind(
            nonsense = "nonsense",
            pathway_class_y = pathway_class_y,
            pathway_class = rev(unique(
              daa_results_filtered_sub_df$pathway_class
            ))
          )
        )
      pathway_class_plot_df$pathway_class_y <-
        as.numeric(pathway_class_plot_df$pathway_class_y)
      pathway_class_annotation <-
        pathway_class_plot_df %>% ggplot2::ggplot(ggplot2::aes(nonsense, pathway_class_y)) + ggplot2::geom_text(
          ggplot2::aes(nonsense, pathway_class_y, label = pathway_class),
          size = 3.5,
          color = "black",
          fontface = "bold",
          family = "sans"
        ) +
        ggplot2::scale_y_discrete(position = "right") +
        ggprism::theme_prism() +
        ggplot2::theme(
          axis.ticks = ggplot2::element_blank(),
          axis.line = ggplot2::element_blank(),
          panel.grid.major.y = ggplot2::element_blank(),
          panel.grid.major.x = ggplot2::element_blank(),
          panel.background = ggplot2::element_blank(),
          axis.text = ggplot2::element_blank(),
          plot.margin = ggplot2::unit(c(0, 0.2, 0, 0), "cm"),
          axis.title.y =  ggplot2::element_blank(),
          axis.title.x = ggplot2::element_blank(),
          legend.position = "non"
        )
    }
    daa_results_filtered_sub_df$p_adjust <-
      as.character(daa_results_filtered_sub_df$p_adjust)
    daa_results_filtered_sub_df$unique <-
      nrow(daa_results_filtered_sub_df) - seq_len(nrow(daa_results_filtered_sub_df)) + 1
    daa_results_filtered_sub_df$p_adjust <-
      substr(daa_results_filtered_sub_df$p_adjust, 1, 5)
    p_annotation <- daa_results_filtered_sub_df %>%
      ggplot2::ggplot(ggplot2::aes(group_nonsense, p_adjust)) +
      ggplot2::geom_text(
        ggplot2::aes(group_nonsense, unique, label = p_adjust),
        size = 3.5,
        color = "black",
        fontface = "bold",
        family = "sans"
      ) +
      ggplot2::labs(y = "p-value (adjusted)") +
      ggplot2::scale_y_discrete(position = "right") +
      ggprism::theme_prism() +
      ggplot2::theme(
        axis.ticks = ggplot2::element_blank(),
        axis.line = ggplot2::element_blank(),
        panel.grid.major.y = ggplot2::element_blank(),
        panel.grid.major.x = ggplot2::element_blank(),
        panel.background = ggplot2::element_blank(),
        axis.text = ggplot2::element_blank(),
        plot.margin = ggplot2::unit(c(0, 0.2, 0, 0), "cm"),
        axis.title.y =  ggplot2::element_text(
          size = 11,
          color = "black",
          vjust = 0
        ),
        axis.title.x = ggplot2::element_blank(),
        legend.position = "non"
      )
    if (p_value_bar == TRUE) {
      if (ko_to_kegg == TRUE) {
        combination_bar_plot <-
          pathway_class_annotation + bar_errorbar + p_values_bar + p_annotation + patchwork::plot_layout(ncol = 4, widths =
                                                                                                c(1, 1.2, 0.5, 0.1))
      }
      else{
        combination_bar_plot <-
          bar_errorbar + p_values_bar + p_annotation + patchwork::plot_layout(ncol = 3, widths = c(2.3, 0.7, 0.3))
      }
    }else{
      if (ko_to_kegg == TRUE) {
        combination_bar_plot <-
          pathway_class_annotation + bar_errorbar + p_annotation + patchwork::plot_layout(ncol = 3, widths =
                                                                                                           c(1, 1.2, 0.1))
      }
      else{
        combination_bar_plot <-
          bar_errorbar + p_annotation + patchwork::plot_layout(ncol = 2, widths = c(2.5,  0.2))
      }
    }
    return(combination_bar_plot)
  }
andressamv commented 1 year ago

Sorry for the delay response, but using your function directly worked perfectly! Thank you!