Closed mt1022 closed 11 months ago
est_optimal_codons <- function(seqs, codon_table = get_codon_table()) { aa_code <- qvalue <- pvalue <- . <- codon <- subfam <- NULL # due to NSE notes in R CMD check cf_all <- count_codons(seqs) enc <- get_enc(cf_all, codon_table = codon_table)
# exclude stop codons
codon_table <- codon_table[aa_code != "*"]
cf_all <- cf_all[, colnames(cf_all) %in% codon_table$codon]
# regression analysis for each codon sub-family
binreg <- lapply(split(codon_table$codon, f = codon_table$subfam), function(x) {
cf <- cf_all[, x, drop = FALSE]
if (ncol(cf) == 1) {
data.table(codon = colnames(cf), coef = 0, se = 0, zvalue = 0, pvalue = 0)
} else {
total <- rowSums(cf)
res <- apply(cf, 2, function(x) {
x
fit <- stats::glm(cbind(x, total - x) ~ enc, family = "binomial")
summary(fit)$coefficients[-1, ]
})
res <- data.table::as.data.table(as.data.frame(t(res)), keep.rownames = "codon")
data.table::setnames(res, c("codon", "coef", "se", "zvalue", "pvalue"))
}
})
bingreg <- data.table::rbindlist(binreg, idcol = "subfam")
bingreg[, `:=`(qvalue = stats::p.adjust(pvalue, method = "BH"))]
# Check if columns exist before removing them
cols_to_remove <- c("codon_b1", "codon_b2", "codon_b3", "se", "zvalue")
cols_to_remove <- cols_to_remove[cols_to_remove %in% colnames(bingreg)]
if (length(cols_to_remove) > 0) {
bingreg[, (cols_to_remove) := NULL]
}
return(bingreg)
}
删除est_optimal_codons中的:bingreg[, c('codon_b1', 'codon_b2', 'codon_b3', 'se', 'zvalue') := NULL]
These warnings do not influence the results but should be solved.