search

microbiome / mia

Microbiome analysis
https://microbiome.github.io/mia/
Artistic License 2.0
45 stars 25 forks source link

Error using mia::transformCounts #565

Open CecileBlanchon opened 4 weeks ago

CecileBlanchon commented 4 weeks ago

Hello,

I am new on mia and I want to do Abundance plot. I work with a phyloseq object so I have first transform it on TreeSE using TreeSE_PRJNA556449_1.2 <- mia::makeTreeSEFromPhyloseq(physeq_PRJNA556449_1.2) TreeSE_PRJNA556449_1.2

class: TreeSummarizedExperiment dim: 758 2 metadata(0): assays(1): counts rownames(758): f7a27f7f3420f0b11f58eb518867d01b 2a488db8a698f43f871b55b3f3744c73 ... 9fb7271f930f020eefd76224a9a33cd3 27499fff67a4e6954307eb4868faae7b rowData names(7): Kingdom Phylum ... Genus Species colnames(2): SRR9826960 SRR9826961 colData names(26): Experiment.ID Biosample.ID ... Name Tax.ID reducedDimNames(0): mainExpName: NULL altExpNames(0): rowLinks: a LinkDataFrame (758 rows) rowTree: 1 phylo tree(s) (758 leaves) colLinks: NULL colTree: NULL

rowData(TreeSE_PRJNA556449_1.2)

DataFrame with 758 rows and 7 columns Kingdom Phylum Class Order

f7a27f7f3420f0b11f58eb518867d01b d__Archaea Aenigmarchaeota Deep_Sea_Euryarchaeo.. Deep_Sea_Euryarchaeo.. 2a488db8a698f43f871b55b3f3744c73 d__Bacteria Bacteroidota Bacteroidia Chitinophagales 8fecc5747436fa0ec37ba97f73f60725 d__Bacteria Bacteroidota Bacteroidia Bacteroidales 527d95d66a874bc23f450cb1e2009b56 d__Bacteria Bacteroidota Bacteroidia Bacteroidales 8906d4b2dbc2e4f5fb79c0a1bca3e68a d__Bacteria Bacteroidota Bacteroidia Bacteroidales ... ... ... ... ... Family Genus Species f7a27f7f3420f0b11f58eb518867d01b Deep_Sea_Euryarchaeo.. Deep_Sea_Euryarchaeo.. uncultured_archaeon 2a488db8a698f43f871b55b3f3744c73 Chitinophagaceae Sediminibacterium Sediminibacterium_sp. 8fecc5747436fa0ec37ba97f73f60725 Rikenellaceae Rikenellaceae_RC9_gu.. Unknown_Species_Rike.. 527d95d66a874bc23f450cb1e2009b56 Rikenellaceae Rikenellaceae_RC9_gu.. uncultured_organism 8906d4b2dbc2e4f5fb79c0a1bca3e68a Rikenellaceae Rikenellaceae_RC9_gu.. uncultured_bacterium ... ... ... ...

assay(TreeSE_PRJNA556449_1.2)[1:5,1:2]

SRR9826960 SRR9826961 f7a27f7f3420f0b11f58eb518867d01b 58 53 2a488db8a698f43f871b55b3f3744c73 0 15 8fecc5747436fa0ec37ba97f73f60725 35 0 527d95d66a874bc23f450cb1e2009b56 28 0 8906d4b2dbc2e4f5fb79c0a1bca3e68a 13 0

Then I tried to create transformed the counts into relative abundance but I always get an error.... Can you help me? TreeSE_PRJNA556449_1.2_relativeabundance<-mia::transformCounts(TreeSE_PRJNA556449_1.2, assay_name = "counts", method = "relabundance", name = "relabundance")

Erreur : [matrixStats (>= 1.2.0)] useNames = NA is defunct. Instead, specify either useNames = TRUE or useNames = FALSE. See also ?matrixStats::matrixStats.options

Thanks !

TuomasBorman commented 4 weeks ago

Hello! Thanks, hmmm, can you give here session info? Also, is the data something that you could share?

-Tuomas

CecileBlanchon commented 4 weeks ago

Yes, thanks,

Here are the session info :

R version 4.2.3 (2023-03-15) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Big Sur ... 10.16

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

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

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

other attached packages: [1] miaViz_1.6.0 ggraph_2.2.1 mia_1.6.0
[4] MultiAssayExperiment_1.24.0 TreeSummarizedExperiment_2.6.0 Biostrings_2.66.0
[7] XVector_0.38.0 SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0
[10] Biobase_2.58.0 GenomicRanges_1.50.2 GenomeInfoDb_1.34.9
[13] IRanges_2.32.0 S4Vectors_0.36.2 BiocGenerics_0.44.0
[16] MatrixGenerics_1.10.0 matrixStats_1.1.0 mixOmics_6.22.0
[19] lattice_0.22-6 MASS_7.3-58.2 stringr_1.5.1
[22] microbiome_1.20.0 phyloseq_1.42.0 ggpubr_0.6.0
[25] ggplot2_3.5.1

loaded via a namespace (and not attached): [1] utf8_1.2.4 tidyselect_1.2.1 RSQLite_2.3.7 htmlwidgets_1.6.4
[5] grid_4.2.3 BiocParallel_1.32.6 Rtsne_0.17 munsell_0.5.1
[9] ScaledMatrix_1.6.0 codetools_0.2-20 DT_0.33 withr_3.0.0
[13] colorspace_2.1-0 knitr_1.47 rstudioapi_0.16.0 ggsignif_0.6.4
[17] labeling_0.4.3 GenomeInfoDbData_1.2.9 polyclip_1.10-6 bit64_4.0.5
[21] farver_2.1.2 rhdf5_2.42.1 vctrs_0.6.5 treeio_1.22.0
[25] generics_0.1.3 xfun_0.44 R6_2.5.1 ggbeeswarm_0.7.2
[29] graphlayouts_1.1.1 rsvd_1.0.5 bitops_1.0-7 rhdf5filters_1.10.1
[33] cachem_1.1.0 gridGraphics_0.5-1 DelayedArray_0.24.0 scales_1.3.0
[37] nnet_7.3-19 beeswarm_0.4.0 gtable_0.3.5 beachmat_2.14.2
[41] processx_3.8.4 tidygraph_1.3.0 rlang_1.1.4 splines_4.2.3
[45] rstatix_0.7.2 lazyeval_0.2.2 broom_1.0.6 checkmate_2.3.1
[49] BiocManager_1.30.23 yaml_2.3.8 reshape2_1.4.4 abind_1.4-5
[53] backports_1.5.0 Hmisc_5.1-1 tools_4.2.3 zCompositions_1.5.0-3
[57] ggplotify_0.1.2 decontam_1.18.0 biomformat_1.26.0 RColorBrewer_1.1-3
[61] Rcpp_1.0.12 plyr_1.8.9 base64enc_0.1-3 sparseMatrixStats_1.10.0
[65] zlibbioc_1.44.0 purrr_1.0.2 RCurl_1.98-1.14 ps_1.7.6
[69] rpart_4.1.23 viridis_0.6.5 ggrepel_0.9.5 cluster_2.1.6
[73] fs_1.6.4 DECIPHER_2.26.0 magrittr_2.0.3 data.table_1.15.4
[77] RSpectra_0.16-1 truncnorm_1.0-9 ggnewscale_0.4.10 patchwork_1.2.0
[81] evaluate_0.23 gridExtra_2.3 compiler_4.2.3 scater_1.26.1
[85] ellipse_0.5.0 tibble_3.2.1 crayon_1.5.2 htmltools_0.5.8.1
[89] ggfun_0.1.5 mgcv_1.9-1 corpcor_1.6.10 Formula_1.2-5
[93] tidyr_1.3.1 aplot_0.2.2 DBI_1.2.3 tweenr_2.0.3
[97] Matrix_1.5-3 ade4_1.7-22 car_3.1-2 permute_0.9-7
[101] cli_3.6.2 parallel_4.2.3 igraph_1.5.1 pkgconfig_2.0.3
[105] foreign_0.8-86 scuttle_1.8.4 foreach_1.5.2 ggtree_3.6.2
[109] rARPACK_0.11-0 vipor_0.4.7 DirichletMultinomial_1.40.0 multtest_2.54.0
[113] NADA_1.6-1.1 yulab.utils_0.1.4 callr_3.7.6 digest_0.6.35
[117] vegan_2.6-4 rmarkdown_2.27 tidytree_0.4.6 htmlTable_2.4.2
[121] DelayedMatrixStats_1.20.0 curl_5.2.1 lifecycle_1.0.4 nlme_3.1-164
[125] jsonlite_1.8.8 Rhdf5lib_1.20.0 carData_3.0-5 BiocNeighbors_1.16.0
[129] qiime2R_0.99.6 desc_1.4.3 viridisLite_0.4.2 fansi_1.0.6
[133] pillar_1.9.0 pkgbuild_1.4.4 fastmap_1.2.0 survival_3.6-4
[137] glue_1.7.0 remotes_2.5.0 iterators_1.0.14 bit_4.0.5
[141] ggforce_0.4.2 stringi_1.8.4 blob_1.2.4 BiocSingular_1.14.0
[145] memoise_2.0.1 dplyr_1.1.4 irlba_2.3.5.1 ape_5.7-1

And for the data,

I first create the phyloseq object with the data (https://filesender.renater.fr/?s=download&token=398beb64-fad9-44dd-8e37-63b171312ec5)

physeq_PRJNA556449<-qiime2R::qza_to_phyloseq(features="table_PseudoPooling.qza",
                        tree="rooted-tree-PseudoPooling.qza",
                        taxonomy="taxonomy_Silva_PseudoPooling.qza")
physeq_PRJNA556449
sampledata <-read.csv("sra-metadata.tsv", sep="\t", row.names = 1)
sample_data(physeq_PRJNA556449) <- sampledata

Then I perform some cleaning.

new_rank <- c("Kingdom", "Phylum",  "Class",   "Order",   "Family",  "Genus",   "Species")
colnames(tax_table(physeq_PRJNA556449)) <- new_rank
nochlo = subset_taxa(physeq_PRJNA556449, !Class=="Chloroplast") # delete chloroplastes
nochlo
physeq_PRJNA556449_1.1 = prune_taxa(taxa_sums(nochlo)!=0, nochlo) # delete otu without counts
physeq_PRJNA556449_1.2 <- physeq_PRJNA556449_1.1
tax <- data.frame(tax_table(physeq_PRJNA556449_1.1))
tax.clean <- data.frame(row.names = row.names(tax),
                          Kingdom = str_replace(tax[,1], "D_0__",""),
                          Phylum = str_replace(tax[,2], "D_1__",""),
                          Class = str_replace(tax[,3], "D_2__",""),
                          Order = str_replace(tax[,4], "D_3__",""),
                          Family = str_replace(tax[,5], "D_4__",""),
                          Genus = str_replace(tax[,6], "D_5__",""),
                          Species = str_replace(tax[,7], "D_6__",""),
                          stringsAsFactors = FALSE)
tax.clean[is.na(tax.clean)] <- "Unknown"
for (j in 1:7){tax.clean[,j] <- as.character(tax.clean[,j])}
tax.clean[is.na(tax.clean)] <- ""
for (j in 1:nrow(tax.clean)) {
    if (tax.clean[j,6]=="uncultured"){
      genus <- paste("Uncultured_Genus_", tax.clean[j,5], sep = "")
      tax.clean[j,6:7] <- genus
    }
    if (tax.clean[j,6]=="Unknown"){
      genus <- paste("Unknown_Genus_", tax.clean[j,5], sep = "")
      tax.clean[j,6:7] <- genus
    }
        if (tax.clean[j,7]=="Unknown"){
      genus <- paste("Unknown_Species_", tax.clean[j,6], sep = "")
      tax.clean[j,7] <- genus
    }
  }
  tax_table(physeq_PRJNA556449_1.2) <- as.matrix(tax.clean)

Thanks

TuomasBorman commented 4 weeks ago

Thanks! Hmm, you have very old mia. The current release version is 1.12.0 https://bioconductor.org/packages/release/bioc/html/mia.html

I advice you to first update the mia since the package has improved a lot from your version.