inst/extdata/171003_ATF_fluxscript.R

In adsteen/metafluxr: Process Isotope-Labeled Metabolic Flux Data

library(devtools)
#install_github("adsteen/metablize", ref="dev")
#library(metabalize)

library(ggplot2)
library(reshape2)
library(plyr)
# library(knitr)

# Load package by hand
source("R/calc_12C.R")
source("R/calc_ks_and_plot.R")
source("R/exp_coef.R")
source("R/fit_raw.R")
source("R/generate_exp_guess.R")
source("R/generate_facet_formula.R")
source("R/get_preds.R")
source("R/key_fn.R")
source("R/metabalize.R")
source("R/plot_single_metabolite.R")
source("R/plot_timecourse.R")
source("R/read_long.R")
source("R/read_MAVEN.R")
source("R/safe_NLS_pred.R")
source("R/safe_NLS.R")

#######
# DATA SETUP
#######

# Set data and key filenames #*****# Change filenames
data.fn <- "inst/extdata/abigail_25_fluxdata.csv"
key.fn <- "inst/extdata/salmonella_samplekey.csv"

# Read the raw data #*****# CHange NH4 to glucose or glucosamine or something
raw_list_glucose <- read_MAVEN(data.fn=data.fn, key.fn=key.fn)

# Set the replicates to a factor
raw_glucose <- raw_list_glucose$raw_data
raw_glucose$replicate <- as.factor(raw_glucose$replicate)

# Pull out the N14 peaks
system.time({ #****# DON'T change calc_12C, but DO change raw_NH4 to raw_glu or raw_GLN or whatever
  raw_glucose <- calc_12C(raw_glucose) # Note that calc_12C identifies 12C OR 14N peaks (really, just whetever the lowest mass isotopomer is)
}) # about 6 seconds on my macbook air

# Filter 14N peaks from 15N-containing peaks
C12_only <- subset(raw_glucose, is.12C==TRUE) # Again, is.12C actually tells you whether youve got the 14N isotopomer

##########
# CREATE NLS FITS AND PREDICTIONS
##########

# Make nls fits for each compound
system.time({ # This throws an error - I don't know why
  nls_fits_12C <- dlply(C12_only,
                    c("compound", "treatment", "sample.type"), safe_NLS)
}) # ~1.5 seconds on my laptop, not so bad

#### Make 'predicted' fits through each data set (i.e., the fit line)

# Define the 'domain' for predicted fits
dom <- seq(from=0, to=max(C12_only$time), length.out = 100)
# Make the predicted fits
pred_fits_12C <- ldply(nls_fits_12C, safe_NLS_pred, domain=dom)


# Find the number of cases where the comuputer couldn't fit
bad_fits <- is.na(nls_fits_12C)
frac_bad_fits <- sum(bad_fits)/length(nls_fits_12C)
print(frac_bad_fits)
############
# DATA OUTPUT: Plotting and writing table of values
############

# Make one, giant graph of the 14N-vs-time data
p_raw_huge <- ggplot(C12_only, aes(x=time, y=relative.ion.count, colour=replicate, shape=treatment)) +
  geom_point() +
  geom_line() +
  facet_grid(treatment~compound, scales="free_y") +
  theme(text=element_text(size=6))
print(p_raw_huge) #****# Change filename for plot
ggsave("plots/salmonella_glucose_all.png", width=50, height=3, units="in", dpi=150,limitsize=FALSE)

# Make individual plots of each metabolite
all_plot_list <- dlply(C12_only, c("compound", "treatment"), plot_single_metabolite,
                       print.plots=FALSE, save.plots=TRUE, base.fn=("plots/try 1/"), #****# Change NH4_singleplots to glu_singleplots etc
                       predictions=pred_fits_12C,
                       height=2.5, width=4, units="in", dpi=150)

# Get NLS fit parameters
coefs_df <- ldply(nls_fits_12C, exp_coef)
names(coefs_df)[4:8] <- c("A.est", "k.est", "A.std.err", "k.std.err", "n")
# Drop the ones where it couldn't calculate a fit
coefs_df_good <- na.omit(coefs_df)

# Make a .csv file of the fit results
write.csv(coefs_df_good, "data/salmonella_glucose_rates.csv", row.names = FALSE) #***# Change filename


##defined negative in the script