version 0.5.0

cran · Jan 10, 2024 · 9db9a80 · 9db9a80
1 parent 58ca55e
commit 9db9a80
Show file tree

Hide file tree

Showing 18 changed files with 697 additions and 484 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: cubar
 Title: Codon Usage Bias Analysis
-Version: 0.4.2
+Version: 0.5.0
 Authors@R: 
     person(given = "Hong",
            family = "Zhang",
@@ -31,8 +31,8 @@ VignetteBuilder: knitr
 RoxygenNote: 7.2.3
 Config/testthat/edition: 3
 NeedsCompilation: no
-Packaged: 2023-11-18 03:00:19 UTC; mt1022
+Packaged: 2024-01-09 03:24:52 UTC; mt1022
 Author: Hong Zhang [aut, cre, cph] (<https://orcid.org/0000-0002-4064-9432>)
 Maintainer: Hong Zhang <mt1022.dev@gmail.com>
 Repository: CRAN
-Date/Publication: 2023-11-18 15:30:10 UTC
+Date/Publication: 2024-01-09 04:00:02 UTC
diff --git a/MD5 b/MD5
@@ -1,26 +1,26 @@
-21c39d529f2c6c13f2fe3b55704b75df *DESCRIPTION
+6d80040fede76749a44d80c660046a3d *DESCRIPTION
 9bc3d376c8fc8af894cfa9eb6177f596 *LICENSE
 1553d0571fdfb109c110a1deac4ad63a *NAMESPACE
-c452a4382e883791ae5bc418c7f62629 *NEWS.md
+0b64b07cf5c06e70511386a7874e17c5 *NEWS.md
 4a7e366db6aeb2799b1950b7de2cd210 *R/codon_info.R
-4e0d4b3b102817b820ea4e7cdff64d5c *R/codon_weight.R
+d3b971a265db410d20f1ca53e88aac70 *R/codon_weight.R
 8628f2b40a9a7b80faed596d0fcb168e *R/data.R
 05b19aa8f69e4a0e392334bda0086dba *R/gene_index.R
 4bec2b4108e903c4f8365baa29c2125a *R/sequences.R
-ca21156c41a0c212d9381a7e7c1a8cb2 *README.md
-c6ad5532b0349fa09a335683bc47288f *build/vignette.rds
+d96f5e16e526627bcadcf1393131f2b9 *README.md
+78a5b0c52b9c89097f9ca40a1b95aba7 *build/vignette.rds
 1f349d33a5d568d0391852b7967ab029 *data/aa2codon.rda
 20c59defd31164fa0b1f8f0d7b3a453c *data/human_mt.rda
 a1ad4e7cd932e86fe53f82ad3bfae345 *data/yeast_cds.rda
 1640256abff1c6ecedbf894bc835f67c *data/yeast_exp.rda
 3b1a684c94051e25d2c9a717ae1cb952 *data/yeast_half_life.rda
 ac63a1404bc7666aace3dc425340bdfe *data/yeast_trna_gcn.rda
-751b3d128f133b4afe3ec7d23ebf33e8 *inst/doc/Introduction.R
-bc1bc1b40d0c0dbe52c68fb2615d55c5 *inst/doc/Introduction.Rmd
-730a1072dabcc152d66469de4c3b9bd7 *inst/doc/Introduction.html
-897842b2d76822440e54c66f5095e0f0 *inst/doc/Non_standard_genetic_code.R
-d6cc1b5daf923a028a0ae3b9f8420ff7 *inst/doc/Non_standard_genetic_code.Rmd
-27ecf2904ca3b1159e660991a7f689f8 *inst/doc/Non_standard_genetic_code.html
+808a812ea48a05e564c45542f0106418 *inst/doc/cubar.R
+ef7efced21805b4a933266bb00ad4e79 *inst/doc/cubar.Rmd
+7bf0aa085d64fcc987fbb8c6716e92a1 *inst/doc/cubar.html
+47ddbe56b7d1390fdd66ade78282545b *inst/doc/non_standard_genetic_code.R
+1ec2413b38d204e8b1a54a42515a098e *inst/doc/non_standard_genetic_code.Rmd
+723470d96cb00d7f7118b519456569e1 *inst/doc/non_standard_genetic_code.html
 06d9239a330458908f7946eab2e144e5 *man/aa2codon.Rd
 a6561adc65c4bca7d9ec7a5e5c420364 *man/check_cds.Rd
 3e072ee96300049a088a86901f152be7 *man/count_codons.Rd
@@ -49,5 +49,5 @@ e468aeaa30b867f6269d5307e9e22f1d *man/yeast_half_life.Rd
 4bfa782e433bcbb8b08b04daff42cb22 *man/yeast_trna_gcn.Rd
 ce5944b0e391747a5b0141d648c0e557 *tests/testthat.R
 6c7f21ff880577f2f90546df430e3920 *tests/testthat/test-count_codons.R
-bc1bc1b40d0c0dbe52c68fb2615d55c5 *vignettes/Introduction.Rmd
-d6cc1b5daf923a028a0ae3b9f8420ff7 *vignettes/Non_standard_genetic_code.Rmd
+ef7efced21805b4a933266bb00ad4e79 *vignettes/cubar.Rmd
+1ec2413b38d204e8b1a54a42515a098e *vignettes/non_standard_genetic_code.Rmd
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,12 @@
+# cubar 0.5.0
+
+* fixed a bug in `est_trna_weight`. Now zero w values were replaced with
+  geometric mean (rather than the arithmetic mean) of non-zero w values.
+* fixed unexpected warnings in `est_optimal_codons`.
+* fixed bugs that update input codon table due to `data.table` reference
+  semantics.
+* added a new vignette explaining the mathematical details of implementation.
+
 # cubar 0.4.2
 
 * adjust formatting

diff --git a/R/codon_weight.R b/R/codon_weight.R
@@ -56,6 +56,7 @@ plot_ca_pairing <- function(codon_table = get_codon_table(), plot = TRUE){
     anticodon <- codon <- codon_b1 <- codon_b2 <- codon_b3 <- NULL # due to NSE notes in R CMD check
     . <- aa_code <- base_codon <- base_anti <- type <- NULL
     anticodon_aa <- codon_aa <- i.aa_code <- NULL
+    codon_table <- data.table::copy(codon_table)
     codon_table[, anticodon := as.character(rev_comp(codon_table$codon))]
     codon_table[, c('codon_b1', 'codon_b2', 'codon_b3') := data.table::tstrsplit(codon, '')]
     bases <- c('T', 'C', 'A', 'G')
@@ -161,7 +162,11 @@ est_trna_weight <- function(trna_level, codon_table = get_codon_table(),
     dtt_W <- ca_pairs[, .(W = sum(ac_level * (1 - penality))), by = .(codon)]
     codon_table[dtt_W, W := i.W, on = .(codon)]
     codon_table[, w := W/max(W, na.rm = TRUE)]
-    mean_w <- mean(codon_table$w, na.rm = TRUE)
+
+    # using geometric mean of w for rare cases that a codon has no matching tRNA
+    # probably due to incomplete tRNA annotation
+    w0 <- codon_table[!is.na(w), w]
+    mean_w <- exp(mean(log(w0)))
     codon_table[is.na(w), w := mean_w]
     return(codon_table)
 }
@@ -214,7 +219,7 @@ est_optimal_codons <- function(seqs, codon_table = get_codon_table()){
     bingreg <- data.table::rbindlist(binreg, idcol = 'subfam')
     bingreg[, qvalue := stats::p.adjust(pvalue, method = 'BH')]
     bingreg <- codon_table[bingreg, on = .(codon, subfam)]
-    bingreg[, c('codon_b1', 'codon_b2', 'codon_b3', 'se', 'zvalue') := NULL]
+    bingreg[, c('se', 'zvalue') := NULL]
     return(bingreg)
 }
 

diff --git a/README.md b/README.md
@@ -4,7 +4,8 @@
 <!-- badges: start -->
 [![CRAN status](https://www.r-pkg.org/badges/version/cubar)](https://CRAN.R-project.org/package=cubar)
 [![](https://cranlogs.r-pkg.org/badges/cubar)](https://cran.r-project.org/package=cubar)
-[![Lifecycle: stable](https://img.shields.io/badge/lifecycle-stable-brightgreen.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable)
+[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10155990.svg)](https://doi.org/10.5281/zenodo.10155990)
 <!-- badges: end -->
 
 ### Overview
@@ -52,6 +53,13 @@ The latest developmental version of `cubar` can be installed with:
 devtools::install_github("mt1022/cubar", dependencies = TRUE)
 ```
 
+### Usage
+Documentation can be found within R (by typing `?function_name`). The following tutorials are available from our [website](https://mt1022.github.io/cubar/):
+
+- [Get Started](https://mt1022.github.io/cubar/articles/cubar.html)
+- [Non-standard Genetic Code](https://mt1022.github.io/cubar/articles/non_standard_genetic_code.html)
+- [Theories behind cubar](https://mt1022.github.io/cubar/articles/theory.html)
+
 ### Suggests
 - [Biostrings](https://bioconductor.org/packages/release/bioc/html/Biostrings.html) for sequence input/output and manipulation.
 - [Peptides](https://github.com/dosorio/Peptides) for peptide- or protein-related indices.

diff --git a/build/vignette.rds b/build/vignette.rds
diff --git a/inst/doc/Introduction.Rmd b/inst/doc/Introduction.Rmd
diff --git a/inst/doc/Non_standard_genetic_code.Rmd b/inst/doc/Non_standard_genetic_code.Rmd
diff --git a/inst/doc/Introduction.R → inst/doc/cubar.R b/inst/doc/Introduction.R → inst/doc/cubar.R
@@ -12,14 +12,19 @@ library(ggplot2)
 ## -----------------------------------------------------------------------------
 # example data
 yeast_cds
+
+# qc
 yeast_cds_qc <- check_cds(yeast_cds)
+yeast_cds_qc
 
+## -----------------------------------------------------------------------------
 # convert a CDS to codon sequence
 seq_to_codons(yeast_cds_qc[['YDR320W-B']])
 
 # convert a CDS to amino acid sequence
 Biostrings::translate(yeast_cds_qc[['YDR320W-B']])
 
+## -----------------------------------------------------------------------------
 # get codon frequency
 yeast_cf <- count_codons(yeast_cds_qc)
 
@@ -45,14 +50,14 @@ plot_dist <- function(x, xlab = 'values'){
 plot_dist(enc, 'ENC')
 
 ## -----------------------------------------------------------------------------
-# estimate optimal codons
-optimal_codons <- est_optimal_codons(yeast_cds_qc, codon_table = ctab)
-head(optimal_codons)
-
 # get fop
 fop <- get_fop(yeast_cds)
 plot_dist(fop, 'Fop')
 
+## -----------------------------------------------------------------------------
+optimal_codons <- est_optimal_codons(yeast_cds_qc, codon_table = ctab)
+head(optimal_codons[optimal_codons$coef < 0 & optimal_codons$qvalue < 0.01, ])
+
 ## -----------------------------------------------------------------------------
 # estimate RSCU of highly expressed genes
 yeast_heg <- head(yeast_exp[order(-yeast_exp$fpkm), ], n = 500)