curveFitting

[UriaMorP]

Hey, I wish to get the output of metalondaAll in some data-structure (especially the curve fits) so I can later plot it myself using other frameworks. I saw that metalondaAll only saves jpg files, so I tried to call curveFittig myself. I got strange errors and when tried to follow the exact snippet as in the vignette I get the same:

Error in gssanova(Count ~ Time, data = group.0, family = "nbinomial",  : 
  gss error in gssanova: unpenalized terms are linearly dependent
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf
5: In min(x, y) : no non-missing arguments to min; returning Inf
6: In max(x, y) : no non-missing arguments to max; returning -Inf

Obviously, the curveFitting function works when called from within metalondaAll, and again, my only wish is to have the data points for the fitted curve and the normalized data in a nice dataframe so I can plot these myself. Can I get some help with that?

My sessionInfo:

R version 3.4.2 (2017-09-28)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS: /apps/RH7U2/gnu/R/3.4.2/lib64/R/lib/libRblas.so
LAPACK: /apps/RH7U2/gnu/R/3.4.2/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] dplyr_0.7.8          readr_1.3.1          MetaLonDA_1.1.0      BiocInstaller_1.28.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0           lattice_0.20-38      prettyunits_1.0.2    gtools_3.8.1         utf8_1.1.4          
 [6] assertthat_0.2.0     glmnet_2.0-16        digest_0.6.18        foreach_1.4.4        R6_2.3.0            
[11] plyr_1.8.4           stats4_3.4.2         RSQLite_2.1.1        httr_1.4.0           ggplot2_3.1.0       
[16] pillar_1.3.1         gplots_3.0.1         rlang_0.3.1          progress_1.2.0       lazyeval_0.2.1      
[21] gdata_2.18.0         rstudioapi_0.9.0     blob_1.1.1           S4Vectors_0.16.0     Matrix_1.2-15       
[26] labeling_0.3         stringr_1.3.1        RCurl_1.95-4.11      bit_1.1-14           biomaRt_2.34.2      
[31] munsell_0.5.0        compiler_3.4.2       pkgconfig_2.0.2      BiocGenerics_0.24.0  tidyselect_0.2.5    
[36] tibble_2.0.1         gss_2.1-9            IRanges_2.12.0       codetools_0.2-16     matrixStats_0.54.0  
[41] XML_3.98-1.16        fansi_0.4.0          crayon_1.3.4         bitops_1.0-6         grid_3.4.2          
[46] gtable_0.2.0         DBI_1.0.0            magrittr_1.5         scales_1.0.0         KernSmooth_2.23-15  
[51] cli_1.0.1            stringi_1.2.4        doParallel_1.0.14    bindrcpp_0.2.2       metagenomeSeq_1.20.1
[56] limma_3.34.9         RColorBrewer_1.1-2   iterators_1.0.10     tools_3.4.2          bit64_0.9-7         
[61] Biobase_2.38.0       glue_1.3.0           purrr_0.2.5          hms_0.4.2            parallel_3.4.2      
[66] yaml_2.2.0           AnnotationDbi_1.40.0 colorspace_1.3-2     caTools_1.17.1.1     memoise_1.1.0       
[71] bindr_0.1.1  

Thanks!

[aametwally]

Hi,

I've updated the package on Github to include the fitted model in the output of metalonda and metalondaAll.

First, download/install the latest development code of MetaLonDA (v1.1.5) from GitHub:

library(devtools)
install_github("aametwally/MetaLonDA", ref = "master")

Follow instructions on README file to run MetaLonDA. For example, if you run it on one feature as in:

output.metalonda.f5 = metalonda(Count = metalonda_test_data[5,], Time = Time, Group = Group,
                                ID = ID, n.perm = 100, fit.method = "nbinomial", points = points,
                                text = rownames(metalonda_test_data)[5], parall = FALSE, pvalue.threshold = 0.05,
                                adjust.method = "BH", time.unit = "days", ylabel = "Normalized Count",
                                col = c("black", "green"), prefix = "Test_F5")

All details about the fitted model can be extracted from "output.metalonda.f5$model". The estimated points of each group are saved in data frames dd.0 and dd.1. So, you can access them through

output.metalonda.f5$model$dd.0
output.metalonda.f5$model$dd.1

Similarly for metalondaAll example:

output.metalonda.all = metalondaAll(Count = metalonda_test_data, Time = Time, Group = Group,
                                    ID = ID, n.perm = 100, fit.method = "nbinomial", num.intervals = 100, 
                                    parall = FALSE, pvalue.threshold = 0.05, adjust.method = "BH", 
time.unit ="hours", 
                                    norm.method = "none", prefix = "Test_metalondaALL", ylabel = "Read Counts",
                                    col = c("black", "green"))

You can access the estimated points from first feature's model using:

output.metalonda.all$output.model[[1]]$dd.0
output.metalonda.all$output.model[[1]]$dd.1

Moreover, in this updated MetaLoDA version, all internal calculations are exported in .RData format so they can be accessed later on.

please let me know if I can be of any further help.

[UriaMorP]

Thanks for the quick response! I'll follow your instructions and give an update!

[UriaMorP]

Hey @aametwally , Sorry for the delayed answer! It works perfectly. Another (unrelated) issue I noticed is that now, in the CurveFit_nbinomial.jpg files, the curve for actual data and the ones for fitted values have the same transparency value so they cannot be distinguished.

Thanks, Uria

[aametwally]

Wonderful!

Regarding the fitted spline figure, I'll be adding this enhancement to the next version, which should be released in less than a month.

Thanks,