# Hello, world!
#
# This is an example function named 'hello'
# which prints 'Hello, world!'.
#
# You can learn more about package authoring with RStudio at:
#
# http://r-pkgs.had.co.nz/
#
# Some useful keyboard shortcuts for package authoring:
#
# Build and Reload Package: 'Cmd + Shift + B'
# Check Package: 'Cmd + Shift + E'
# Test Package: 'Cmd + Shift + T'
read_peaks <- function(fname, format = "xlsx", ...) {
if (format == "xlsx"){
raw <- readxl::read_excel(fname, ...)
}
else if (format == "csv"){
raw <- readr::read_csv(fname, ...)
}
else{
stop("unsupported file format", fname)
}
return(raw)
}
#' renames kathrins dataset to useful names and reshapres a little. Specific to one file.
rename_kathrin <-function(df){
df %>% dplyr::rename(
minutes = `Time (min)`,
start_time = `flg22 added`,
cell_number = `cell#`,
tech_rep = TechRepID,
bio_rep = BioRepID
) %>%
tidyr::gather('YFP', 'CFP', key = "reporter", value = "intensity") %>%
dplyr::mutate(has_started = minutes >= start_time) %>%
dplyr::mutate(cell_number = as.character(cell_number),
tech_rep = as.character(tech_rep),
bio_rep = as.character(bio_rep)
)
}
#' adds smoothed intensity for each curve. Time before `has_started` counts as seperate curve
#' @export
smooth <- function(df, col_to_smooth = ratio, ...){
col_to_smooth <- dplyr::enquo(col_to_smooth)
df %>%
dplyr::filter(! is.na(!! col_to_smooth )) %>%
dplyr::mutate(smoothed = smooth.spline(!!col_to_smooth)$y )
}
#' @export
ratio_curve <- function(df, value = corrected_intensity, key = reporter, test_reporter = YFP, base_reporter = CFP ){
value <- dplyr::enquo(value)
key <- dplyr::enquo(key)
test_reporter <- dplyr::enquo(test_reporter)
base_reporter <- dplyr::enquo(base_reporter)
df %>%
tidyr::spread(!!key, !!value) %>%
dplyr::mutate(ratio = !!test_reporter / !!base_reporter)
}
get_proportions <- function(df){
return( df %>% dplyr::group_by(cell_number, tech_rep, bio_rep, line) %>%
tidyr::spread(reporter, intensity) %>%
dplyr::mutate(total = YFP + CFP,
prop_yfp = YFP / total,
prop_cfp = CFP / total
)
)
}
#' converts a grouped_dataframe to a list, each member of which is an original group
grouped_df_to_list <- function(gdf){
indices <- attributes(gdf)$indices
df <- dplyr::ungroup(gdf) %>%
dplyr::mutate(idx = rownames(gdf) )
groups <- lapply(indices, function(x,df){df[x + 1,]}, df = df )
return(groups)
}
find_loess_points <- function(group_list,
#grouped_df,
x_col = "minutes",
y_col = "intensity",
reporter_col = "reporter",
test_reporter_name = "YFP",
base_reporter_name = "CFP",
test_reporter_quantile = 0.75,
base_reporter_quantile = 0.25 ){
#group_list <- grouped_df_to_list(grouped_df)
result <- lapply(group_list, find_loess_points_ ,
x_col = x_col,
y_col = y_col,
reporter_col = reporter_col,
test_reporter_name = test_reporter_name,
base_reporter_name = base_reporter_name,
test_reporter_quantile = test_reporter_quantile,
base_reporter_quantile = base_reporter_quantile)
return(result)
}
find_loess_points_ <- function(df, x_col = "minutes",
y_col = "intensity",
reporter_col = "reporter",
test_reporter_name = "YFP",
base_reporter_name = "CFP",
test_reporter_quantile = 0.75,
base_reporter_quantile = 0.25
){
mod <- lm(y ~ x, data = data.frame(x = df[[x_col]], y = df[[y_col]] ))
if ( df[[reporter_col]][[1]] == base_reporter_name ){
limit <- quantile( resid(mod), c(base_reporter_quantile) )
peaks <- 1:length(df[[y_col]]) %in% quantmod::findPeaks(df[[y_col]])
df$is_loess_point <- (resid(mod) > limit) & peaks
return( df )
}
if ( df[[reporter_col]][[1]] == test_reporter_name ){
limit <- quantile( resid(mod), c(test_reporter_quantile) )
valleys <- 1:length(df[[y_col]]) %in% quantmod::findValleys(df[[y_col]])
df$is_loess_point <- (resid(mod) > limit) & valleys
return( df )
}
}
get_corrected <- function(df,
x_col = "minutes",
y_col = "intensity"){
points <- dplyr::filter(df, is_loess_point == TRUE )
mod_d <- tibble::tibble(x = points[[x_col]], y = points[[y_col]])
baseline_mod <- loess(y ~ x, data = mod_d )
values <- predict( baseline_mod, newdata = data.frame( x = df[[x_col]] ))
values <- abs(values - max(values, na.rm = TRUE))
df$corrected_intensity <- df[[y_col]] + values
return(df)
}
correct_intensity <- function(grouped_df,
x_col = "minutes",
y_col = "intensity",
reporter_col = "reporter",
test_reporter_name = "YFP",
base_reporter_name = "CFP",
test_reporter_quantile = 0.75,
base_reporter_quantile = 0.25
){
#points_col_quo <- dplyr::enquo(points_column)
#points <- dplyr::filter(df, !! points_col_quo == TRUE )
group_list <- grouped_df_to_list(grouped_df)
loess_point_list <- lapply(group_list, find_loess_points_ ,
x_col = x_col,
y_col = y_col,
reporter_col = reporter_col,
test_reporter_name = test_reporter_name,
base_reporter_name = base_reporter_name,
test_reporter_quantile = test_reporter_quantile,
base_reporter_quantile = base_reporter_quantile)
corrected_list <- lapply(loess_point_list, get_corrected, x_col = x_col, y_col = y_col)
return(dplyr::bind_rows(corrected_list))
}
#' @export
find_peak_floor <- function(df){
df %>% dplyr::group_by(reporter, cell_number, tech_rep, bio_rep) %>%
dplyr::summarize(peak_floor = mean(smoothed) )
}
single_curve <- function(df){
df %>% dplyr::filter(cell_number == "1", tech_rep == "1", bio_rep == "1", line == "WT")
}
#' df. Applies smmoth spline to a column called corrected_intensity
smooth_corrected <- function(df, col_to_smooth = corrected_intensity){
smooth_col_quo <- dplyr::enquo(col_to_smooth)
dplyr::filter(df, ! is.na( !!smooth_col_quo) ) %>%
dplyr::mutate(smoothed = smooth.spline(!! smooth_col_quo )$y )
}
find_peaks <- function(df){
quantmod::findPeaks(df[['smoothed']])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.