<repo_name>santiago1234/corrcodon<file_sep>/R/compute_correlations.R
import::from(dplyr, "%>%")

#' Computes the correlation of tile and codon enrichment
#'
#' \code{compute_correlations} returns the correlations of tile (rna expression)
#' with codon enrichment for all the codons starting at \code{min_tile}
#'
#' This function computes the correlation (spearman or pearson)
#' between the tile number (rna expression) and the codon enrichment value (codon usage
#' transcriptome level) starting at \code{min_tile}. A positive correlation indicates that the
#' codon is depleted at genes with low expression and enriched at genes with high expression
#' and vice versa for a negative correlation, a correlation near zero may indicate
#' that there is not an enrichment pattern for that particular codon. The \code{min_tile} value
#' is \strong{important}, it indicates from where to start  computing the correlation,
#' this min tile value should be where we start having gene expression, (if we choose a tile
#' value of 60 all the genes in tiles less than 60 will not be used in the computation), for instance
#' if we have 100 tiles and we pick \code{min_tile} to be 1, the first tiles will be genes
#' that are potentially no expressed and it will add noise to the resulting value. We shoul
#' include a \code{min_tile} value where the genes start to get expressed (low expression but not abssent expression).
#' @param enrichment A tibble,  enrichment of each codon at each tile, this
#' tibble is the output of the function \code{codon_enrichmen}.
#' @param min_tile An integer, tile less than min_tile will not be used in the computation.
#' @param corr_method A character string, either spearman or pearson
#' @return A tibble with the correlation of each codon.
#' @examples
#' encode_rna_seq %>%
#'   subset_genes(hsapiens_codon_composition) %>%
#'   tile_expression(n = 100) %>%
#'   quantify_codons_tile(hsapiens_codon_composition) %>%
#'   codon_enrichment() -> enrichment
#'   compute_correlations(enrichment, 60)
#' @note Correlations for the stop codons are misleading
compute_correlations <- function(enrichment, min_tile, corr_method = "spearman") {
  if (!min_tile %in% enrichment$tile) {
    stop("wrong tile value, see ?compute_correlations")
  }
  correlation <- function(codon_) {
    enrichment %>%
      dplyr::filter(codon == codon_, tile >= min_tile) -> selected_cols
    cor(selected_cols$tile, selected_cols$enrichment, method = corr_method) -> correlation_result

    tibble::tibble(codon = codon_, correlation = correlation_result)
  }

  CODONS %>%
    purrr::map_df(correlation)
}
<file_sep>/R/data.R
#' Example of expected data input for the functions in this packages
#'
#' This RNA quantification data sets comes from the encode project
#'
#' @format A tibble: 61,471 × 2
#' \describe{
#'   \item{gene_id}{ensmble, gene ids}
#'   \item{gene_quantification, trancript quantifications TPMs}
#'   ...
#' }
#' @source \url{https://github.com/santiago1234/ccrnaec/tree/master/human-encode-analysis/data}
"encode_rna_seq"

#' H. Sapiens open reading frames codon composition
#'
#' This data set contains the codon composition of the human genes
#' the longest transcript isoform
#'
#' @format A tibble: 19,673 × 67
#' \describe {
#'   \item{gene_id}{ensmble gene id}
#'   \item{cds_length}{coding sequence length in nucleotides}
#'   \item{codon}{the frequency of each codon in the cds}
#'   ...
#'}
#' @source \url{https://github.com/santiago1234/ccrnaec/tree/master/data}
"hsapiens_codon_composition"

#' Genetic Code
#'
#' the DNA codons
#'
#' @format A charcter vector of length 64
"CODONS"


<file_sep>/tests/testthat/test_codon_enrichment.R
# test the function codon_enrichment
import::from(dplyr, "%>%")

encode_rna_seq %>%
  subset_genes(hsapiens_codon_composition) %>%
  tile_expression(n = 10) %>%
  quantify_codons_tile(codon_composition = hsapiens_codon_composition) %>%
  codon_enrichment() -> results_to_test

test_that("enrichment values add 1", {
  results_to_test %>%
    dplyr::group_by(tile) %>%
    dplyr::summarise(add_one = sum(enrichment)) %>%
    .[["add_one"]] -> ones
  tile <- results_to_test$tile %>% max
  expect_equal(sum(ones), tile)
})
<file_sep>/R/corrcodon_pipeline.R
import::from(dplyr, "%>%")

#' Pipeline function to compute the codon enrichment correlation with gene expression
#' @inheritParams subset_genes
#' @return A tibble with the correlation values
#' @examples
#' subset_genes(encode_rna_seq, hsapiens_codon_composition) %>% corrcodon_pipeline(hsapiens_codon_composition)
corrcodon_pipeline <- function(rna_experiment, codon_composition, tiles = 100, min_tile = 59) {
  message("running pipeline ...")

  if (any(!rna_experiment$gene_id %in% codon_composition$gene_id)) {
    warning("running function subset_genes ..., not all genes ids of rna experiment in codon composition file" )
  }
  rna_experiment %>%
    subset_genes(codon_composition) %>%
    tile_expression(n = tiles) %>%
    quantify_codons_tile(codon_composition = codon_composition) %>%
    codon_enrichment() %>%
    compute_correlations(min_tile = min_tile, corr_method = "spearman")

}
<file_sep>/tests/testthat/test_tile_expression.R
import::from(dplyr, "%>%")

## test the function subset_genes
experiment_subset <- subset_genes(encode_rna_seq, hsapiens_codon_composition)
test_that("gene ids functio output subset of codon_composition", {
  subset_ids <- experiment_subset$gene_id
  contained <- any(!subset_ids %in% hsapiens_codon_composition$gene_id)
  expect_false(contained)
})


## test the tile expression function
tiles <- 20
f_out <- tile_expression(encode_rna_seq, n = tiles)

test_that("added a new column to the data frame", {
  expect_equal(ncol(f_out), ncol(encode_rna_seq) + 1)
  expect_true("tile" %in% colnames(f_out))
})

test_that('ranges tiles', {
  expect_equal(min(f_out$tile), 1)
  expect_equal(max(f_out$tile), tiles)
})

## test the function expected behaviour
## sum gene_quantifications in tile i should always be less
## than sum ene_quantifications in tile i + 1

test_that(
  "expected behaviour",
  {
     for (i in 1:(tiles - 1)) {
       quanitfications_i <- f_out$gene_quantification[f_out$tile == i]
       quanitfications_i_plus_1 <- f_out$gene_quantification[f_out$tile == (i + 1)]
       expect_true(sum(quanitfications_i) <= sum(quanitfications_i_plus_1))
     }
  }

)
<file_sep>/R/tile_expression.R
#' Gets the genes that are present in both files: codon composition & rna experiment
#'
#' \code{subset_genes} returns a subset of rna_experiment where all the ids are present in
#' the codon composition file
#'
#' This is the first function to be run, if the returned file is empty we are not using the
#' same gene ids in both tables, gene_ids from codon composition should match
#' gene ids in rna_seq experiment
#'
#' @param rna_experiment A tibble with the gene quantifications (TPMs or RPKMS)
#' in a collumn called gene_quantification
#' @param codon_composition A tibble with with the codon composition, 66 columns,
#' densities of each codon, cds_length & gene_id
#' @return A tibble where all the gene ids in rna_exeriment and contained in gene
#' ids codon composition
#' @examples
#' subset_genes(encode_rna_seq, hsapiens_codon_composition)
subset_genes <- function(rna_experiment, codon_composition) {
  rna_experiment[rna_experiment$gene_id %in% codon_composition$gene_id, ]
}


#' Groups gene expression in tiles
#'
#' @param rna_experiment A tibble with the gene quantifications (TPMs or RPKMS) in a collumn called gene_quantification
#' @param n An interger, number of tiles
#' @return A tibble with a new column: tile
tile_expression <- function(rna_experiment, n) {

  if (!"gene_quantification" %in% colnames(rna_experiment)) {
    stop(
      "not column called gene_quantification, provided a collumn with the genes quantification (TPMs or RPKMs)"
      )
  }

    dplyr::mutate(
      rna_experiment,
      tile = dplyr::ntile(gene_quantification, n = n)
    )
}
<file_sep>/tests/testthat.R
library(testthat)
library(corrcodon)

test_check("corrcodon")
<file_sep>/tests/testthat/test_count_codons_in_tiles.R
import::from(dplyr, "%>%")

## test the functions in count_codons_in_tiles.R
til <- 50
subset_genes(encode_rna_seq, hsapiens_codon_composition) %>%
  tile_expression(n = til) -> experiment


# test tha function uses the apposite genes -------------------------------

test_that(
  "output the same genes that are in the tile",
  {
    message("testing function get_codon_expression_tile ...")
    test_f <- function(tile_number) {
      n_genes <- dplyr::filter(experiment, tile == tile_number) %>% .$gene_id
      get_codon_expression_tile(
        experiment, codon = "AAA", tile_value = tile_number, hsapiens_codon_composition
      ) %>%
        .$gene_id -> out_genes
      expect_false(any(!n_genes %in% out_genes))
    }
    purrr::map(
      1:til,
      test_f
    )
})


# test function quantify codons tile --------------------------------------

codon_quantification <- quantify_codons_tile(experiment,codon_composition =  hsapiens_codon_composition)

test_that(
  "outputs the correct number of codons",
  {
    # sum all codons
    test_f2 <- function(tile_number) {
      get_codon_expression_tile(
        experiment, codon = "AAA", tile_value = tile_number, hsapiens_codon_composition
      ) -> out_f
      (out_f$number_of_codons * out_f[["AAA"]]) %>% sum -> tc
      expect_true(dplyr::filter(codon_quantification, codon == "AAA", tile == tile_number)[["total_codons"]] == tc)
    }
    purrr::map(
      1:til,
      test_f2
    )
  }
)
<file_sep>/R/codon_enrichment.R
import::from(dplyr, "%>%")

#' Compute the enrichment of each codon at each tile
#'
#' \code{codon_enrichment} returns a table with the enrichment of each codon at
#' each tile
#'
#' This function computes the enrichment of each codon in each tile, the enrichment
#' is computed: \deqn{E_{codon} = \frac{codon_in_tile}{total_codons_in_tile}}
#' @param codon_quantifications a tibble, the output of function \code{quantify_codons_tile}
#' @return A tibble, with enrichment column added indicating the value of the enrichment
#' of each particular codon in tile
#' @examples
#' encode_rna_seq %>%
#'   subset_genes(hsapiens_codon_composition) %>%
#'   tile_expression(n = 100) %>%
#'   quantify_codons_tile(hsapiens_codon_composition) %>%
#'   codon_enrichment()
codon_enrichment <- function(codon_quantifications) {
  # compute the total number of codons in each tile
  codon_quantifications %>%
    dplyr::group_by(tile) %>%
    dplyr::summarise(total_in_tile = sum(total_codons)) -> total_tiles

  enrichment_tile <- function(tile_value) {
    # compute the total codons in the tile
    total_tiles %>%
      dplyr::filter(tile == tile_value) %>%
      .[["total_in_tile"]] -> tile_total

    # compute the codon enrichment in tile
    dplyr::filter(codon_quantifications, tile == tile_value) %>%
      dplyr::mutate(enrichment = total_codons / tile_total)
  }

  # apply the above function to all the tiles
  1:max(codon_quantifications$tile) %>%
    purrr::map_df(.f = enrichment_tile)

}
<file_sep>/R/count_codons_in_tiles.R
## FUNCTIONS TO COUNT THE NUMBER OF CODONS IN EACH TILE

## import the pipe with the import package
import::from(dplyr, "%>%")

#' Joins the codon composition with the rna experiment
#'
#' \code{get_codon_expression_tile} generates the tibble with columns gene_id,
#' number_of_codons, codon, gene_quantification for the tile
#'
#' This function combines the codon_composition file and the rna_experiment by
#' selecting the genes in tile: tile_value and then reports the columns:
#' gene_id, number_of_codons, codon, gene_quantification, Also the function
#' calculates the number of codons per gene.
#' @param rna_experiment_tiled The output of function \code{tile_expression}
#' @param codon A string one of the codons
#' @param tile_value The tile to subset the data, integer
#' @param codon_composition A tibble with with the codon composition, 66 columns
#' @examples
#' tile_data <- encode_rna_seq %>% subset_genes(hsapiens_codon_composition) %>% tile_expression(n = 10)
#' get_codon_expression_tile(tile_data, codon = "AAA", tile_value = 8, codon_composition = hsapiens_codon_composition)
get_codon_expression_tile <- function(rna_experiment_tiled, codon, tile_value, codon_composition) {
  tile_x <- dplyr::filter(rna_experiment_tiled, tile == tile_value)
  codon_composition %>%
    dplyr::mutate(
      number_of_codons = cds_length / 3
    ) %>%
    dplyr::select(
      gene_id,
      number_of_codons,
      get(codon)
    ) %>%
    dplyr::inner_join(
      dplyr::select(tile_x, gene_id, gene_quantification),
      by = "gene_id"
    )

}


# codons_in_tile ----------------------------------------------------------

#' Computes the total number of codons in tile
#'
#' This function runs \code{get_codon_expression_tile} and computes the total number
#' of codons in the tile: the sum of the codons of each gene in the tile
#'
#' @inheritParams get_codon_expression_tile
#' @return A tibble: 1 × 3, codon, total_codons tile
#' @examples
#' tile_data <- encode_rna_seq %>% subset_genes(hsapiens_codon_composition) %>% tile_expression(n = 10)
#' codons_in_tile(tile_data, codon = "AAA", tile_value = 8, codon_composition = hsapiens_codon_composition)
#' @author <NAME>
codons_in_tile <- function(rna_experiment_tiled, codon, tile_value, codon_composition) {
   get_codon_expression_tile(rna_experiment_tiled, codon, tile_value, codon_composition) %>%
    dplyr::mutate(num_codon = .$number_of_codons * .[[codon]]) -> codons_in_tile_per_gene

  dplyr::tibble(
    codon = codon,
    total_codons = sum(codons_in_tile_per_gene$num_codon),
    tile = tile_value
  )

 }


# quantify_codons_tile ----------------------------------------------------

#' Counts the codons in each tile
#'
#' This function runs \code{codons_in_tile} for all the codons and all tiles
#' values.
#' @param rna_experiment_tiled rna_experiment_tiled The output of function \code{tile_expression}
#' @param codon_composition A tibble with with the codon composition, 66 columns
#' @return A tibble with the number of codons in each tile per codon
#' @examples
#' subset_genes(encode_rna_seq, hsapiens_codon_composition) %>%
#'   tile_expression(n = 100) %>%
#'   quantify_codons_tile(hsapiens_codon_composition)

quantify_codons_tile <- function(rna_experiment_tiled, codon_composition) {
  ## send warning if not all genes id of experiment are in codon_composition gene ids
  if(any(!rna_experiment_tiled$gene_id %in% codon_composition$gene_id)) {
    warning("not all ids of experiment in codon composition ids, you should run: subset_genes before this :/")
  }

  max_tile <- max(rna_experiment_tiled$tile)
  dplyr::tibble(
  codones = rep(CODONS, max_tile),
  tile = rep(1:max_tile, 64) %>% sort
  ) -> to_map


  ## parallel version
  parallel::mclapply(
    1:nrow(to_map),
    FUN = function(index) codons_in_tile(
      rna_experiment_tiled,
      codon = to_map$codones[index],
      tile_value = to_map$tile[index],
      codon_composition
    ),
    mc.cores = 4
  ) %>%
    dplyr::bind_rows()

}


