diff --git a/omics/kelly/workshop.qmd b/omics/kelly/workshop.qmd index 697ad89..3810f06 100644 --- a/omics/kelly/workshop.qmd +++ b/omics/kelly/workshop.qmd @@ -136,6 +136,9 @@ vfa_cummul <- vfa_cummul |> sep = "-", remove = FALSE) ``` +📢 This code depends on the `sample_replicate` column being in the form treatment-replicate. In the sample data CN10 and NC are the treatments. The replicate is a number from 1 to 3. The value does include a encoding for time. You might want to edit your file to match this format. + + The provided data is cumulative/absolute. We need to calculate the change in VFA with time. There is a function, `lag()` that will help us @@ -147,7 +150,7 @@ sure the data is in the right order so we will arrange by ## 1. Calculate *Change* in VFA g/l with time -🎬 Create dataframe for the change in VFA +🎬 Create dataframe for the change in VFA 📢 and the change in time ```{r} vfa_delta <- vfa_cummul |> @@ -160,11 +163,15 @@ vfa_delta <- vfa_cummul |> isopentanoate = isopentanoate - lag(isopentanoate), pentanoate = pentanoate - lag(pentanoate), isohexanoate = isohexanoate - lag(isohexanoate), - hexanoate = hexanoate - lag(hexanoate)) + hexanoate = hexanoate - lag(hexanoate), + delta_time = time_day - lag(time_day)) ``` Now we have two dataframes, one for the cumulative data and one for the -change in VFA. +change in VFA and time. Note that the VFA values have been replaced by the change in VFA but the change in time is in a separate column. I have done this because we later want to plot flux (not yet added) against time + +📢 This code also depends on the `sample_replicate` column being in the form treatment-replicate. lag is calculating the difference between a value at one time point and the next for a treatment-replicate combination. + ## 2. Recalculate the data into grams per litre @@ -209,14 +216,15 @@ View `vfa_cummul` to check you understand what you have done. Repeat for the delta data. -🎬 Pivot the change data, `delta_vfa` to long format: +🎬 Pivot the change data, `delta_vfa` to long format (📢 delta_time is added to the list of columns that do not need to be pivoted but repeated): ```{r} vfa_delta <- vfa_delta |> pivot_longer(cols = -c(sample_replicate, treatment, replicate, - time_day), + time_day, + delta_time), values_to = "conc_mM", names_to = "vfa") ``` @@ -515,13 +523,39 @@ hexanoate is especially high in the NC replicate 2 at time = 22 Calculate the flux(change in VFA concentration over a period of time, divided by weight or volume of material) of each VFA, by mM and by -weight. +weight. +Emma's note: I think the terms flux and reaction rate are used +interchangeably I've requested clarification: *for the flux measurements, do they need graphs of the rate of change wrt time? And is the sludge volume going to be a constant for all samples or something they measure and varies by vial?* +Answer: **The sludge volume is constant, at 30 mls within a 120ml vial. +Some students will want to graph reaction rate with time, others will +want to compare the measured GC-FID concentrations against the model +output.** + +📢 Kelly asked for ".. a simple flux measurement, which is the change in VFA concentration over a period of time, divided by weight or volume of material. In this case it might be equal to == Delta(Acetate at 3 days - Acetate at 1 day)/Delta (3days - 1day)/50 mls sludge. This would provide a final flux with the units of mg acetate per ml sludge per day." + +Note: Kelly says mg/ml where earlier he used g/L. These are the same (but I called +my column `conc_g_l`) + +We need to use the `vfa_delta` data frame. It contains the change in VFA concentration and the change in time. We will add a column for the flux of each VFA in g/L/day. (mg/ml/day) + +```{r} + +sludge_volume <- 30 # ml +vfa_delta <- vfa_delta |> + mutate(flux = conc_g_l / delta_time / sludge_volume) + +``` + +NAs at time 1 are expected because there's no time before that to calculate a changes + + + ## 5. Graph and extract the reaction rate - pending Graph and extract the reaction rate assuming a first order @@ -532,6 +566,75 @@ fitting, I assume x is time but I'm not clear what the Y variable is - concentration? or change in concentration? or rate of change of concentration?* +Answer: **The non-linear equation describes concentration change with +time. Effectively the change in concentration is dependent upon the +available concentration, in this example \[Hex\] represents the +concentration of Hexanoic acid, while the T0 and T1 represent time +steps.** + +**\[Hex\]T1 = \[Hex\]T0 - \[Hex\]T0 \* k** + +**Or. the amount of Hexanoic acid remaining at T1 (let's say one hour +from the last data point) is equal to the starting concentration +(\[Hex\]T0) minus the concentration dependent metabolism (\[Hex\]To \* +k).** + +📢 We can now plot the observed fluxes (reaction rates) over time + +I've summarised the data to add error bars and means +```{r} +vfa_delta_summary <- vfa_delta |> + group_by(treatment, time_day, vfa) |> + summarise(mean_flux = mean(flux), + se_flux = sd(flux)/sqrt(length(flux))) |> + ungroup() +``` + + +```{r} +ggplot(data = vfa_delta, aes(x = time_day, colour = vfa)) + + geom_point(aes(y = flux), alpha = 0.6) + + geom_errorbar(data = vfa_delta_summary, + aes(ymin = mean_flux - se_flux, + ymax = mean_flux + se_flux), + width = 1) + + geom_errorbar(data = vfa_delta_summary, + aes(ymin = mean_flux, + ymax = mean_flux), + width = 0.8) + + scale_color_viridis_d(name = NULL) + + scale_x_continuous(name = "Time (days)") + + scale_y_continuous(name = "VFA Flux mg/ml/day") + + theme_bw() + + facet_wrap(~treatment) + + theme(strip.background = element_blank()) +``` + +Or maybe this is easier to read: + +```{r} +ggplot(data = vfa_delta, aes(x = time_day, colour = treatment)) + + geom_point(aes(y = flux), alpha = 0.6) + + geom_errorbar(data = vfa_delta_summary, + aes(ymin = mean_flux - se_flux, + ymax = mean_flux + se_flux), + width = 1) + + geom_errorbar(data = vfa_delta_summary, + aes(ymin = mean_flux, + ymax = mean_flux), + width = 0.8) + + scale_color_viridis_d(name = NULL, begin = 0.2, end = 0.7) + + scale_x_continuous(name = "Time (days)") + + scale_y_continuous(name = "VFA Flux mg/ml/day") + + theme_bw() + + facet_wrap(~ vfa, nrow = 2) + + theme(strip.background = element_blank(), + legend.position = "top") +``` + + +I have not yet worked out the best way to plot the modelled reaction rate + Pages made with R [@R-core], Quarto [@allaire2022], `knitr` [@knitr], `kableExtra` [@kableExtra] diff --git a/renv.lock b/renv.lock index a098a78..3fa2bb5 100644 --- a/renv.lock +++ b/renv.lock @@ -1,6 +1,6 @@ { "R": { - "Version": "4.3.3", + "Version": "4.3.1", "Repositories": [ { "Name": "BioCsoft", @@ -196,14 +196,14 @@ }, "DBI": { "Package": "DBI", - "Version": "1.2.2", + "Version": "1.1.3", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R", "methods" ], - "Hash": "164809cd72e1d5160b4cb3aa57f510fe" + "Hash": "b2866e62bab9378c3cc9476a1954226b" }, "DESeq2": { "Package": "DESeq2", @@ -336,7 +336,7 @@ }, "MASS": { "Package": "MASS", - "Version": "7.3-60.0.1", + "Version": "7.3-60", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -347,16 +347,15 @@ "stats", "utils" ], - "Hash": "b765b28387acc8ec9e9c1530713cb19c" + "Hash": "a56a6365b3fa73293ea8d084be0d9bb0" }, "Matrix": { "Package": "Matrix", - "Version": "1.6-5", + "Version": "1.5-4.1", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R", - "grDevices", "graphics", "grid", "lattice", @@ -364,7 +363,7 @@ "stats", "utils" ], - "Hash": "8c7115cd3a0e048bda2a7cd110549f7a" + "Hash": "38082d362d317745fb932e13956dccbb" }, "MatrixGenerics": { "Package": "MatrixGenerics", @@ -842,7 +841,7 @@ }, "cluster": { "Package": "cluster", - "Version": "2.1.6", + "Version": "2.1.4", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -852,7 +851,7 @@ "stats", "utils" ], - "Hash": "0aaa05204035dc43ea0004b9c76611dd" + "Hash": "5edbbabab6ce0bf7900a74fd4358628e" }, "codetools": { "Package": "codetools", @@ -1330,14 +1329,14 @@ }, "glue": { "Package": "glue", - "Version": "1.7.0", + "Version": "1.6.2", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R", "methods" ], - "Hash": "e0b3a53876554bd45879e596cdb10a52" + "Hash": "4f2596dfb05dac67b9dc558e5c6fba2e" }, "googledrive": { "Package": "googledrive", @@ -1703,7 +1702,7 @@ }, "lattice": { "Package": "lattice", - "Version": "0.22-5", + "Version": "0.21-8", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1714,7 +1713,7 @@ "stats", "utils" ], - "Hash": "7c5e89f04e72d6611c77451f6331a091" + "Hash": "0b8a6d63c8770f02a8b5635f3c431e6b" }, "lazyeval": { "Package": "lazyeval", @@ -1728,7 +1727,7 @@ }, "lifecycle": { "Package": "lifecycle", - "Version": "1.0.4", + "Version": "1.0.3", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1737,7 +1736,7 @@ "glue", "rlang" ], - "Hash": "b8552d117e1b808b09a832f589b79035" + "Hash": "001cecbeac1cff9301bdc3775ee46a86" }, "limma": { "Package": "limma", @@ -1819,7 +1818,7 @@ }, "mgcv": { "Package": "mgcv", - "Version": "1.9-1", + "Version": "1.8-42", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1832,7 +1831,7 @@ "stats", "utils" ], - "Hash": "110ee9d83b496279960e162ac97764ce" + "Hash": "3460beba7ccc8946249ba35327ba902a" }, "mime": { "Package": "mime", @@ -1875,7 +1874,7 @@ }, "nlme": { "Package": "nlme", - "Version": "3.1-164", + "Version": "3.1-162", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1885,7 +1884,7 @@ "stats", "utils" ], - "Hash": "a623a2239e642806158bc4dc3f51565d" + "Hash": "0984ce8da8da9ead8643c5cbbb60f83e" }, "openssl": { "Package": "openssl", @@ -2020,17 +2019,16 @@ }, "progress": { "Package": "progress", - "Version": "1.2.3", + "Version": "1.2.2", "Source": "Repository", "Repository": "CRAN", "Requirements": [ - "R", "R6", "crayon", "hms", "prettyunits" ], - "Hash": "f4625e061cb2865f111b47ff163a5ca6" + "Hash": "14dc9f7a3c91ebb14ec5bb9208a07061" }, "promises": { "Package": "promises", @@ -2104,7 +2102,7 @@ }, "readr": { "Package": "readr", - "Version": "2.1.5", + "Version": "2.1.4", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2123,7 +2121,7 @@ "utils", "vroom" ], - "Hash": "9de96463d2117f6ac49980577939dfb3" + "Hash": "b5047343b3825f37ad9d3b5d89aa1078" }, "readxl": { "Package": "readxl", @@ -2180,7 +2178,7 @@ }, "reprex": { "Package": "reprex", - "Version": "2.1.0", + "Version": "2.0.2", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2198,7 +2196,7 @@ "utils", "withr" ], - "Hash": "1425f91b4d5d9a8f25352c44a3d914ed" + "Hash": "d66fe009d4c20b7ab1927eb405db9ee2" }, "reshape2": { "Package": "reshape2", @@ -2267,7 +2265,7 @@ }, "rvest": { "Package": "rvest", - "Version": "1.0.4", + "Version": "1.0.3", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2280,9 +2278,10 @@ "rlang", "selectr", "tibble", + "withr", "xml2" ], - "Hash": "0bcf0c6f274e90ea314b812a6d19a519" + "Hash": "a4a5ac819a467808c60e36e92ddf195e" }, "sass": { "Package": "sass", @@ -2300,23 +2299,21 @@ }, "scales": { "Package": "scales", - "Version": "1.3.0", + "Version": "1.2.1", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R", "R6", "RColorBrewer", - "cli", "farver", - "glue", "labeling", "lifecycle", "munsell", "rlang", "viridisLite" ], - "Hash": "c19df082ba346b0ffa6f833e92de34d1" + "Hash": "906cb23d2f1c5680b8ce439b44c6fa63" }, "scran": { "Package": "scran", @@ -2468,7 +2465,7 @@ }, "stringi": { "Package": "stringi", - "Version": "1.8.3", + "Version": "1.7.12", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2477,11 +2474,11 @@ "tools", "utils" ], - "Hash": "058aebddea264f4c99401515182e656a" + "Hash": "ca8bd84263c77310739d2cf64d84d7c9" }, "stringr": { "Package": "stringr", - "Version": "1.5.1", + "Version": "1.5.0", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2494,7 +2491,7 @@ "stringi", "vctrs" ], - "Hash": "960e2ae9e09656611e0b8214ad543207" + "Hash": "671a4d384ae9d32fc47a14e98bfa3dc8" }, "sys": { "Package": "sys", @@ -2547,7 +2544,7 @@ }, "tidyr": { "Package": "tidyr", - "Version": "1.3.1", + "Version": "1.3.0", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2566,7 +2563,7 @@ "utils", "vctrs" ], - "Hash": "915fb7ce036c22a6a33b5a8adb712eb1" + "Hash": "e47debdc7ce599b070c8e78e8ac0cfcf" }, "tidyselect": { "Package": "tidyselect", @@ -2626,14 +2623,14 @@ }, "timechange": { "Package": "timechange", - "Version": "0.3.0", + "Version": "0.2.0", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R", "cpp11" ], - "Hash": "c5f3c201b931cd6474d17d8700ccb1c8" + "Hash": "8548b44f79a35ba1791308b61e6012d7" }, "tinytex": { "Package": "tinytex", @@ -2678,13 +2675,13 @@ }, "utf8": { "Package": "utf8", - "Version": "1.2.4", + "Version": "1.2.3", "Source": "Repository", "Repository": "CRAN", "Requirements": [ "R" ], - "Hash": "62b65c52671e6665f803ff02954446e9" + "Hash": "1fe17157424bb09c48a8b3b550c753bc" }, "uuid": { "Package": "uuid", @@ -2824,10 +2821,10 @@ }, "yaml": { "Package": "yaml", - "Version": "2.3.8", + "Version": "2.3.7", "Source": "Repository", "Repository": "CRAN", - "Hash": "29240487a071f535f5e5d5a323b7afbd" + "Hash": "0d0056cc5383fbc240ccd0cb584bf436" }, "zlibbioc": { "Package": "zlibbioc",