In haiyueliu/Eskrate: Estimate time-dependent RNA kinetics

knitr::opts_chunk$set(cache=TRUE,cache.lazy=FALSE,message=FALSE,warning=FALSE, echo=FALSE, include=TRUE)
knitr::opts_knit$set(root.dir = "~/Documents/Eskrate/data/")

library(tibble)
library(dplyr)
library(magrittr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(Rmpfr)
library(gam)
library(mgcv)
library(foreach)
library(doParallel)
library(Eskrate)

Input data

1. cell info data containing cell_id and cc_time.

2. raw counts data containing cell_id, gene, capture_rate, p_u, p_total, m_u, m_total.

### parameters
# molecules.per.cell <- 10^6L
# cell.cycle.duration <- 19.33
# minutes.per.hour <- 60L
## cells with cell cycle time sorted
# cells <- readRDS(file.path(data.path, "cells_info.rds"))
# cells %<>%
#   arrange(cc_time) %>%
#   mutate(cc_time=rep(seq(0,cell.cycle.duration*minutes.per.hour,length.out=nrow(cells)))) %>%
# as.data.frame
# saveRDS(cells, "~/Documents/Eskrate/test_data/cells_info.rds")
# ## raw counts for different RNA types
# raw.counts <- readRDS(file.path(data.path,"raw_counts.rds"))
# 
# ## add cell time (cc_time in minutes)
# raw.counts %<>%
#   filter(cell_id %in% cells$cell_id) %>%
#   droplevels %>%
#   as.data.frame 
# 
# saveRDS(raw.counts, "test_data/raw_counts.rds")
#
# observed.data %<>%
#   left_join(cells) %>%
#   arrange(gene, cc_time) %>%
#   select(-cell_id) %>%
#   as.data.frame
#
#
# number.of.genes <- length((unique(observed.data$gene)))
# cat("Total number of genes in the observation is: ", number.of.genes)
# ## only keep 15 min cells
# number.of.cells <- length(unique(observed.data$cell_id))
# cat("Total number of cells in the observation is: ", number.of.cells)
# 
# saveRDS(observed.data, "observed.data.rds")

Get input data

data.path <- "~/Documents/Eskrate/data/"
raw_counts_labeled_mature <- readRDS(file.path(data.path, "raw_counts_labeled_mature.rds"))
raw_counts_unlabeled_mature <- readRDS(file.path(data.path, "raw_counts_unlabeled_mature.rds"))
raw_counts_labeled_precursor <- readRDS(file.path(data.path, "raw_counts_labeled_precursor.rds"))
raw_counts_unlabeled_precursor <- readRDS(file.path(data.path, "raw_counts_unlabeled_precursor.rds"))
cells_info <- readRDS(file.path(data.path, "cells_info.rds"))

observed.data <- 
  getInput(rawCountsLabeledMature = raw_counts_labeled_mature,
           rawCountsUnlabeledMature = raw_counts_unlabeled_mature,
           rawCountsLabeledPrecursor = raw_counts_labeled_precursor,
           rawCountsUnlabeledPrecursor = raw_counts_unlabeled_precursor,
           cells = cells_info,
           labelingTime = 15)

### select two genes (ZNF711, KIF14) for demo
observed.data %<>% 
  filter(gene %in% c("ZNF711", "KIF14")) %>%
  droplevels

Normalization

norm.data <- 
  normalizeCountsToConcentration(observed.data)

Smooth normalized gene expression in sorted cells

smooth.data <- smoothGeneProfileByPenalizedSpline(norm.data,
                                                  numberOfAnchorPoints = 20,
                                                  gamma=1.4)

Estimate RNA kinetics

rates.data <- 
  estimateKinetics(smoothedObservedData = smooth.data,
                   threadN=2)

Calculate Predictions

predicted.data <-
  getPredictions(dataRates = rates.data,
                 threadN=2)

Plot predicted gene expression and estimated rates

For example, plot PCNA gene expression and rates.

p1 <-
  plotMtotalAlphaGamma(df.prediction = predicted.data,   
                       gene.to.plot = "ZNF711",
                       color.alpha = "#D95F02",
                       color.gamma = "#1B9E77",
                       color.expression = "#666666",
                       show.phase.boundaries = FALSE,
                       line.size = 1,
                       font.size = 16)

p2 <-
  plotMtotalAlphaGamma(df.prediction = predicted.data,   
                       gene.to.plot = "KIF14",
                       color.alpha = "#D95F02",
                       color.gamma = "#1B9E77",
                       color.expression = "#666666",
                       show.phase.boundaries = FALSE,
                       line.size = 1,
                       font.size = 16)

plot_grid(p1, p2, ncol=2)

haiyueliu/Eskrate documentation built on Sept. 3, 2023, 3:33 p.m.