Examples

Example 1: Basic Usage

In the first example we will correct a measured response vector of the M-57 alanine fragment. First the chemical formula of the fragment and the response vector have to be defined, after which the uncorrected mass isotopomer distribution (MID) can be calculated.

using IsotopeCorrection

ala_m57 = "C11H26N1O2Si2LabC3"
response_vec = [2500000, 2000000, 2500000, 800000]
MID = response_vec ./ sum(response_vec)

4-element Vector{Float64}:
 0.32051282051282054
 0.2564102564102564
 0.32051282051282054
 0.10256410256410256

The function isotope_correction()` can be directly used on a single measurement vector. Here we are assuming that the tracer used for this experiment had an 80% purity.

corr_resp, corr_MID, mean_enrich, residuum = isotope_correction(response_vec,
                                                                ala_m57;
                                                                tracer_purity = 0.8)

([3.0266026323459675e6, 819243.5134320183, 3.666924789018221e6, 587930.8175325774], [0.37362227679544946, 0.10113241278096713, 0.4526675467300171, 0.07257776369356636], 0.5560501993304251, [-2.022641591536693e-13, -2.999336291582156e-13, 3.788572473403735e-13, 2.28532231771029e-13])

corr_resp` is the corrected measurement vector,

corr_resp

4-element Vector{Float64}:
      3.0266026323459675e6
 819243.5134320183
      3.666924789018221e6
 587930.8175325774

corr_MID is the MID of the corrected measurement vector,

corr_MID

4-element Vector{Float64}:
 0.37362227679544946
 0.10113241278096713
 0.4526675467300171
 0.07257776369356636

mean_enrich is the mean enrichment of the corrected measurement vector,

mean_enrich

0.5560501993304251

and residuum is the residals of the optimization.

residuum

4-element Vector{Float64}:
 -2.022641591536693e-13
 -2.999336291582156e-13
  3.788572473403735e-13
  2.28532231771029e-13

The optimization can also be turned off, using the optimization = false keyword argument, though this isn't recommended as the optimization prevents negative values in the corrected response.

To visualize the effects of isotope_correction(), the MID can be plotted before and after correction.

using ColorSchemes, CairoMakie

colors = colorschemes[:bamako][1:64:end];
group_color = [PolyElement(color = color, strokecolor = :transparent)
    for color in colors]
labels = ["m+0", "m+1", "m+2", "m+3"]
Mplus = [0, 1, 2, 3]

fig = Figure();
ax = Axis(fig[1,1],
    xticks = (1:2, ["uncorrected", "corrected"]),
    title = "Alanine (M-57) MID")
barplot!(ax, repeat([1, 2], inner = length(corr_MID)),
    cat(MID, corr_MID, dims = 1),
    stack = repeat(Mplus, 2),
    colormap = colors,
    color = repeat(Mplus .+ 1, 2))

Legend(fig[1,2], group_color, labels, framevisible = false)
fig

Example 2: Usage with a DataFrame

For this example we will make a very simple DataFrame with 2 samples for each of which 2 methionine fragments were measured, specifically M-159 and M-57. We need the chemical formulas of both fragments, as shown in the code below, as well as measured response vectors. For this example the measured response for sample 1 is all possible m+ being recorded with equal frequency while the measured response for sample 2 is simply a random vector.

using IsotopeCorrection
using DataFrames, DataFramesMeta

df = DataFrame(
    Sample = repeat([1, 2], inner = 11),
    AminoAcidFragment = repeat(cat(repeat(["Met_(M-159)"], 5), repeat(["Met_(M-57)"], 6), dims = 1), 2),
    Molecule = repeat(cat(repeat(["C10H24N1S1Si1LabC4"], 5), repeat(["C13H30N1O2S1Si2LabC5"], 6), dims = 1), 2),
    MZ0 = repeat(cat(repeat([218], 5), repeat([320], 6), dims = 1), 2),
    Mplus = repeat(cat(repeat(0:4, 1), repeat(0:5, 1), dims = 1), 2),
    Response = cat(repeat([3e6/5], 5), repeat([3e6/6], 6), rand(0:3e6, 11), dims = 1))

# Of course normally the first step would be to read in (and prepare) real data,
# using CSV.jl or XLSX.jl, for example.

# using CSV
# df = CSV.read("DATA.csv", DataFrame)

Using DataFramesMeta the DataFrame can be grouped based on the sample ID and amino acid fragment and each group corrected using the main function of this package: isotope_correction()

Since this function outputs the corrected response, MID, mean enrichment and residuum, we index into element 1 to recieve only the corrected response. Additionally we are assuming that the tracer used for this experiment had a 99% purity.

gdfs = groupby(df, [:Sample, :AminoAcidFragment])
df_corrected = @transform(gdfs,
    :CorrectedResponse = isotope_correction(Vector(:Response), :Molecule[1]; tracer_purity = 0.99)[1])

Now, again using DataFramesMeta, both the corrected and uncorrected MIDs can be calculated.

gdfs_corrected = groupby(df_corrected, [:Sample, :AminoAcidFragment])
df_corrected_MIDs = @transform(gdfs_corrected,
    :MID = :Response ./ sum(:Response),
    :CorrectedMID = :CorrectedResponse ./ sum(:CorrectedResponse))
select(df_corrected_MIDs, :Sample, :AminoAcidFragment, :Response, :CorrectedResponse, :MID, :CorrectedMID)
# And of course the finished DataFrame can then be saved.
# CSV.write("Corrected_MIDs.csv", df_corrected_MIDs)

To visualize what isotope_correction() does we can plot the MIDs before and after correction.

using ColorSchemes, CairoMakie

gdfs = @groupby(df_corrected_MIDs, :Sample, :AminoAcidFragment)

colors = colorschemes[:bamako][1:43:256];
group_color = [PolyElement(color = color, strokecolor = :transparent)
    for color in colors]
labels = ["m+0", "m+1", "m+2", "m+3", "m+4", "m+5"]

fig = Figure();

function myplot(group, fig, colors)
    name = group.AminoAcidFragment[1]
    sample = group.Sample[1]
    ax = Axis(fig,
        xticks = (1:2, ["uncorrected", "corrected"]),
        title = "Sample $sample $name MID")
    barplot!(ax, repeat([1, 2], inner = length(group.MID)),
        cat(group.MID, group.CorrectedMID, dims = 1),
        stack = repeat(group.Mplus, 2),
        colormap = colors,
        color = repeat(group.Mplus .+ 1, 2))
end
myplot(gdfs[1], fig[1,1], colors[1:5])
myplot(gdfs[2], fig[1,2], colors)
myplot(gdfs[3], fig[2,1], colors[1:5])
myplot(gdfs[4], fig[2,2], colors)

Legend(fig[1:2,3], group_color, labels, framevisible = false)

fig
# save("Methionine_MIDs.png", fig)

This page was generated using Literate.jl.