FTICR Visualizations

# Clean workspace
rm(list=ls())

library("ftmsRanalysis")
library("vegan")
library("ggplot2")
library("ape")
library("ggVennDiagram")
library("reshape2")
library("plyr")

Load Peak Object

load("ficusPeakObj.Robj")
str(peakObj)

List of 3
 $ e_data:'data.frame': 3704 obs. of  33 variables:
  ..$ Mass: chr [1:3704] "200.0564425" "200.092833" "201.004077" "201.0404701" ...
  ..$ H1  : num [1:3704] 21.7 20.2 20.2 22.3 NA ...
  ..$ H2  : num [1:3704] NA NA NA 20.3 NA ...
  ..$ H3  : num [1:3704] 21.3 NA NA 22.2 NA ...
  ..$ H4  : num [1:3704] NA NA NA 20.9 NA ...
  ..$ H5  : num [1:3704] 21.1 NA NA 21.8 20.3 ...
  ..$ H6  : num [1:3704] 20.3 NA 20.2 21.2 NA ...
  ..$ H7  : num [1:3704] 21.1 NA NA 21.5 NA ...
  ..$ H8  : num [1:3704] NA NA NA 20.7 NA ...
  ..$ H9  : num [1:3704] 21.2 20.6 20.3 22.7 NA ...
  ..$ H10 : num [1:3704] NA NA NA NA NA ...
  ..$ H11 : num [1:3704] 21.7 20.7 NA 22.2 20.6 ...
  ..$ H12 : num [1:3704] 20.4 NA NA 21.2 NA ...
  ..$ H13 : num [1:3704] 21.5 NA NA 21.8 NA ...
  ..$ H14 : num [1:3704] 20.8 NA 20.2 21.9 NA ...
  ..$ H15 : num [1:3704] 21.8 NA 20.7 22.6 21 ...
  ..$ H16 : num [1:3704] 21 NA NA 21.8 NA ...
  ..$ C1  : num [1:3704] 21.7 20.2 20.4 22.8 NA ...
  ..$ C2  : num [1:3704] NA NA NA 21 NA ...
  ..$ C3  : num [1:3704] 21.1 20.3 20.2 22.4 20.8 ...
  ..$ C4  : num [1:3704] 20.2 NA NA 21.7 NA ...
  ..$ C5  : num [1:3704] 21.2 NA 20.2 22.7 20.3 ...
  ..$ C6  : num [1:3704] 20.3 NA 20.4 21.7 NA ...
  ..$ C7  : num [1:3704] 21.9 NA 20.8 23 20.3 ...
  ..$ C8  : num [1:3704] 20.3 20.2 NA 21.5 NA ...
  ..$ C9  : num [1:3704] 21.5 20.9 20.9 22.8 NA ...
  ..$ C10 : num [1:3704] 20.7 NA NA 21.5 NA ...
  ..$ C11 : num [1:3704] 21.5 20.2 NA 22.7 20.6 ...
  ..$ C12 : num [1:3704] 20.3 NA NA 20.8 NA ...
  ..$ CS13: num [1:3704] 21.8 20.9 20.2 22.6 20.5 ...
  ..$ C14 : num [1:3704] NA NA NA 20.6 NA ...
  ..$ C15 : num [1:3704] 22 20.8 20.9 22.9 20.7 ...
  ..$ C16 : num [1:3704] NA NA NA 20.9 NA ...
 $ f_data:'data.frame': 32 obs. of  7 variables:
  ..$ Sample_ID  : chr [1:32] "H1" "H2" "H3" "H4" ...
  ..$ Sample.name: chr [1:32] "H2O" "H2M" "H4O" "H4M" ...
  ..$ Horizon    : chr [1:32] "Organic" "Mineral" "Organic" "Mineral" ...
  ..$ Treatment  : chr [1:32] "Heated" "Heated" "Heated" "Heated" ...
  ..$ T_H        : chr [1:32] "Organic_Heated" "Mineral_Heated" "Organic_Heated" "Mineral_Heated" ...
  ..$ Sample.ID  : chr [1:32] "BW.H.2.O" "BW.H.2.M" "BW.H.4.O" "BW.H.4.M" ...
  ..$ wMT        : chr [1:32] "Y" "Y" "Y" "Y" ...
 $ e_meta:'data.frame': 3704 obs. of  27 variables:
  ..$ Mass         : chr [1:3704] "200.0564425" "200.092833" "201.004077" "201.0404701" ...
  ..$ C            : int [1:3704] 8 9 7 8 7 9 8 10 9 11 ...
  ..$ H            : int [1:3704] 11 15 6 10 10 14 14 18 18 22 ...
  ..$ O            : int [1:3704] 5 4 7 6 5 5 4 4 3 3 ...
  ..$ N            : int [1:3704] 1 1 0 0 2 0 2 0 2 0 ...
  ..$ C13          : int [1:3704] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ S            : int [1:3704] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ P            : int [1:3704] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ Error_ppm    : num [1:3704] 0.01752 -0.00742 -0.00455 -0.04212 0.02847 ...
  ..$ NeutralMass  : num [1:3704] 201 201 202 202 202 ...
  ..$ MolForm      : chr [1:3704] "C8H11O5N" "C9H15O4N" "C7H6O7" "C8H10O6" ...
  ..$ AI           : num [1:3704] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ AI_Mod       : num [1:3704] 0.111 0 0.429 0.2 0 ...
  ..$ DBE          : num [1:3704] 4 3 5 4 4 3 3 2 2 1 ...
  ..$ DBE_O        : num [1:3704] -1 -1 -2 -2 -1 -2 -1 -2 -1 -2 ...
  ..$ DBE_AI       : num [1:3704] -2 -2 -2 -2 -3 -2 -3 -2 -3 -2 ...
  ..$ GFE          : num [1:3704] 53.2 73 27.7 53.2 35.9 ...
  ..$ kmass.CH2    : num [1:3704] 200 200 201 201 201 ...
  ..$ kdefect.CH2  : num [1:3704] 0.167 0.131 0.22 0.184 0.173 ...
  ..$ NOSC         : num [1:3704] 0.25 -0.444 1.143 0.25 0.857 ...
  ..$ OtoC_ratio   : num [1:3704] 0.625 0.444 1 0.75 0.714 ...
  ..$ HtoC_ratio   : num [1:3704] 1.375 1.667 0.857 1.25 1.429 ...
  ..$ NtoC_ratio   : num [1:3704] 0.125 0.111 0 0 0.286 ...
  ..$ PtoC_ratio   : num [1:3704] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ NtoP_ratio   : num [1:3704] NA NA NA NA NA NA NA NA NA NA ...
  ..$ ElComposition: chr [1:3704] "CHON" "CHON" "CHO" "CHO" ...
  ..$ bs1_class    : chr [1:3704] "Lignin;Amino Sugar" "Protein" "Tannin" "Tannin" ...
 - attr(*, "cnames")=List of 29
  ..$ edata_cname     : chr "Mass"
  ..$ fdata_cname     : chr "Sample_ID"
  ..$ mass_cname      : chr "Mass"
  ..$ extraction_cname: NULL
  ..$ mf_cname        : chr "MolForm"
  ..$ c_cname         : chr "C"
  ..$ h_cname         : chr "H"
  ..$ o_cname         : chr "O"
  ..$ n_cname         : chr "N"
  ..$ s_cname         : chr "S"
  ..$ p_cname         : chr "P"
  ..$ isotopic_cname  : chr "C13"
  ..$ o2c_cname       : chr "OtoC_ratio"
  ..$ h2c_cname       : chr "HtoC_ratio"
  ..$ kmass_cname     : Named chr "kmass.CH2"
  .. ..- attr(*, "names")= chr "CH2"
  ..$ kdefect_cname   : Named chr "kdefect.CH2"
  .. ..- attr(*, "names")= chr "CH2"
  ..$ nosc_cname      : chr "NOSC"
  ..$ gfe_cname       : chr "GFE"
  ..$ mfname_cname    : NULL
  ..$ aroma_cname     : chr "AI"
  ..$ modaroma_cname  : chr "AI_Mod"
  ..$ dbe_cname       : chr "DBE"
  ..$ dbeo_cname      : chr "DBE_O"
  ..$ dbeai_cname     : chr "DBE_AI"
  ..$ elcomp_cname    : chr "ElComposition"
  ..$ n2c_cname       : chr "NtoC_ratio"
  ..$ p2c_cname       : chr "PtoC_ratio"
  ..$ n2p_cname       : chr "NtoP_ratio"
  ..$ bs1class_cname  : chr "bs1_class"
 - attr(*, "data_info")=List of 2
  ..$ data_scale     : chr "log2"
  ..$ instrument_type: chr "12T"
 - attr(*, "class")= chr [1:2] "peakData" "ftmsData"
 - attr(*, "filters")=List of 3
  ..$ massFilt    :List of 3
  .. ..$ report_text: chr "A mass filter was applied to the data, removing Masses that had a mass less than 200 or a mass greater than 900"| __truncated__
  .. ..$ threshold  : num [1:2] 200 900
  .. ..$ filtered   : chr [1:6869] "98.656748" "98.833389" "98.833638" "98.85848" ...
  ..$ moleculeFilt:List of 3
  .. ..$ report_text: chr "A molecule filter was applied to the data, removing Masses that were present in fewer than 2 samples. A total o"| __truncated__
  .. ..$ threshold  : num 2
  .. ..$ filtered   : chr [1:4258] "200.020042" "200.032592" "200.040329" "200.046097" ...
  ..$ formulaFilt :List of 3
  .. ..$ report_text: chr "A formula filter was applied to the data, removing Masses that had NoFormula assigned. A total of 612 Masses we"| __truncated__
  .. ..$ threshold  : chr "NoFormula"
  .. ..$ filtered   : chr [1:612] "200.9892554" "201.053514" "202.027165" "202.0376923" ...

peakObj

peakData object
# Peaks: 3704
# Samples: 32
Meta data columns: [Mass, C, H, O, N, C13, S, P, Error_ppm, NeutralMass, MolForm, AI, AI_Mod, DBE, DBE_O, DBE_AI, GFE, kmass.CH2, kdefect.CH2, NOSC, OtoC_ratio, HtoC_ratio, NtoC_ratio, PtoC_ratio, NtoP_ratio, ElComposition, bs1_class]
A mass filter was applied to the data, removing Masses that had a mass less than 200 or a mass greater than 900. A total of 6869 Masses were filtered out of the dataset by this filter.
A molecule filter was applied to the data, removing Masses that were present in fewer than 2 samples. A total of 4258 Masses were filtered out of the dataset by this filter.
A formula filter was applied to the data, removing Masses that had NoFormula assigned. A total of 612 Masses were filtered out of the dataset by this filter.

Binary Matrix

# Set rownames
edata <- peakObj$e_data
rownames(edata) <- as.character(edata$Mass)
edata$Mass <- NULL
# Binary matrix
edata.binary <- apply(edata,2,function(y)as.numeric(y>0))
# Na's as zero
edata.binary[is.na(edata.binary)] <- 0
# As dataframe edata.binary
edata.binary <- as.data.frame(edata.binary)
# Add Mass to binary matrix
edata.binary$Mass <- rownames(edata)
# Add metadata to binary matrix
tmp <- peakObj$e_meta[,c("Mass","bs1_class")] # save to tmp; direct merge isn't working
str(tmp)

'data.frame':   3704 obs. of  2 variables:
 $ Mass     : chr  "200.0564425" "200.092833" "201.004077" "201.0404701" ...
 $ bs1_class: chr  "Lignin;Amino Sugar" "Protein" "Tannin" "Tannin" ...

# Merge 
edata.binary <- merge(edata.binary,tmp,by="Mass",all = T)
# Remove Mass column from df
edata.binary$Mass <- NULL

# Carbon compounds #
carbon.cpds <- aggregate(.~bs1_class,data=edata.binary,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using bs1_class as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

# Color palette #
pal12 <- c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99","#E31A1C", "#FDBF6F", "#FF7F00", "#CAB2D6", "#6A3D9A","#FFFF99", "#B15928")
# New color palette #
pal4 <- c("#1F78B4", "#FF7F00", "#33A02C", "#E31A1C")
# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=bs1_class,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=3,position=position_dodge(width=0.75), color="red") +
      ylab("Percentage") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(bs1_class,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Ordination Plot:Treatment

rm(list = ls())
load("ficusPeakObj.Robj")
peakObj <- group_designation(peakObj, main_effects=c("Treatment"))
pcoaMat <- getPrincipalCoordinates(ftmsObj = peakObj,dist_metric = "bray")

Data has already been transformed. Data is on the log2 scale. Changing to pres scale.

plotPrincipalCoordinates(pcoaMat,x=1,y=2,ftmsObj = peakObj )

minimal value for n is 3, returning requested palette with 3 different levels
minimal value for n is 3, returning requested palette with 3 different levels
minimal value for n is 3, returning requested palette with 3 different levels
minimal value for n is 3, returning requested palette with 3 different levels

Ordination Plot:Treatment x Horizon

peakObj <- group_designation(peakObj, main_effects=c("T_H"))
#getGroupDF(peakObj)
pcoaMat <- getPrincipalCoordinates(ftmsObj = peakObj,dist_metric = "bray")

Data has already been transformed. Data is on the log2 scale. Changing to pres scale.

plotPrincipalCoordinates(pcoaMat,x=1,y=2,ftmsObj = peakObj )

Venny Diagram

rm(list = ls())
# load organic
load("organicUnique.Robj")
# get values
organic.unique.df <- organic_unique$e_data
# load mineral
load("mineralUnique.Robj")
# get values
mineral.unique.df <- mineral_unique$e_data
# Organic data
All.organic <- organic.unique.df[organic.unique.df$uniqueness_nsamps %in% "Observed in Both","Mass"]
# Mineral data
All.mineral <- mineral.unique.df[mineral.unique.df$uniqueness_nsamps %in% "Observed in Both","Mass"]
# Create List
venn.data <- list(
  Control.Mineral = c(mineral.unique.df[mineral.unique.df$uniqueness_nsamps %in% "Unique to Control","Mass"],All.mineral),
  Control.Organic = c(organic.unique.df[organic.unique.df$uniqueness_nsamps %in% "Unique to Control","Mass"],All.organic),
  Heated.Mineral = c(mineral.unique.df[mineral.unique.df$uniqueness_nsamps %in% "Unique to Heated","Mass"],All.mineral),
  Heated.Organic = c(organic.unique.df[organic.unique.df$uniqueness_nsamps %in% "Unique to Heated","Mass"],All.organic)
  )

# Venn plot
ggVennDiagram(venn.data,label_alpha = 0)

Elements Plot

Load Peak Object

load("ficusPeakObj.Robj")

Binary Matrix

# Set rownames
edata <- peakObj$e_data
rownames(edata) <- as.character(edata$Mass)
edata$Mass <- NULL
# Binary matrix
edata.binary <- apply(edata,2,function(y)as.numeric(y>0))
# Na's as zero
edata.binary[is.na(edata.binary)] <- 0
# As dataframe edata.binary
edata.binary <- as.data.frame(edata.binary)
# Add Mass to binary matrix
edata.binary$Mass <- rownames(edata)
# Add metadata to binary matrix
tmp <- peakObj$e_meta[,c("Mass","ElComposition")] # save to tmp; direct merge isn't working
str(tmp)

'data.frame':   3704 obs. of  2 variables:
 $ Mass         : chr  "200.0564425" "200.092833" "201.004077" "201.0404701" ...
 $ ElComposition: chr  "CHON" "CHON" "CHO" "CHO" ...

# Merge 
edata.binary <- merge(edata.binary,tmp,by="Mass",all = T)
# Remove Mass column from df
edata.binary$Mass <- NULL

# Carbon compounds #
carbon.cpds <- aggregate(.~ElComposition,data=edata.binary,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

Elements Plot

# New color palette #
pal4 <- c("#1F78B4", "#FF7F00", "#33A02C", "#E31A1C")

# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=ElComposition,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=5,position=position_dodge(width=0.75), color="red") +  
      ylab("Relative Frequency (%)") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

Elements Stats

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(ElComposition,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Nitrogen Compounds

# Filter by N-containing elements
tmp <- edata.binary[grep(pattern = "N",edata.binary$ElComposition),]
# Carbon compounds #
carbon.cpds <- aggregate(.~ElComposition,data=tmp,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

Nitrogen Elements: Plot

# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=ElComposition,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=5,position=position_dodge(width=0.75), color="red") +  
      ylab("Relative Frequency (%)") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

Nitrogen Stats

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(ElComposition,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Phosphorus Compounds

# Filter by N-containing elements
tmp <- edata.binary[grep(pattern = "P",edata.binary$ElComposition),]
# Carbon compounds #
carbon.cpds <- aggregate(.~ElComposition,data=tmp,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

Phosphorus Elements: Plot

# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=ElComposition,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=5,position=position_dodge(width=0.75), color="red") +
      ylab("Relative Frequency (%)") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

Phosphorus Stats

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(ElComposition,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Sulfur Compounds

# Filter by N-containing elements
tmp <- edata.binary[grep(pattern = "S",edata.binary$ElComposition),]
# Carbon compounds #
carbon.cpds <- aggregate(.~ElComposition,data=tmp,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

Sulfur Elements: Plot

# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=ElComposition,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=5,position=position_dodge(width=0.75), color="red") +
      ylab("Relative Frequency (%)") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

Sulfur Stats

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(ElComposition,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Oxygen Compounds

# Filter by N-containing elements
tmp <- edata.binary[grep(pattern = "O",edata.binary$ElComposition),]
# Carbon compounds #
carbon.cpds <- aggregate(.~ElComposition,data=tmp,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,-1] <- decostand(carbon.cpds[,-1],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,-1])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

Oxygen Elements: Plot

# Order levels
carbon.cpds$Horizon <- factor(carbon.cpds$Horizon,levels = c("Organic","Mineral"))

# Plotting
p <- ggplot(carbon.cpds,aes(x=ElComposition,y=value,fill=Treatment))+ geom_boxplot()+
      stat_summary(fun.y=mean, geom="point", shape=20, size=5,position=position_dodge(width=0.75), color="red") +
      ylab("Relative Frequency (%)") + xlab("") + scale_fill_manual(values=pal4) + 
      theme(axis.text.x=element_text(size=12,face="bold",angle = 90,vjust = 0.7),
            axis.text.y=element_text(size=18,face="bold"),
            strip.text.x = element_text(size = 16),
            legend.text=element_text(size=18))

`fun.y` is deprecated. Use `fun` instead.

p + facet_grid(~Horizon)

Oxygen Stats

#stats on relative abundances #
results.temp<-ddply(carbon.cpds,.(ElComposition,Horizon),function(x) summarize(x,Pvalue=t.test(value~Treatment,data=x,na.rm=TRUE,paired=F)$p.value))
results.temp

Map Elements to Compounds

# Set rownames
edata <- peakObj$e_data
rownames(edata) <- as.character(edata$Mass)
edata$Mass <- NULL
# Binary matrix
edata.binary <- apply(edata,2,function(y)as.numeric(y>0))
# Na's as zero
edata.binary[is.na(edata.binary)] <- 0
# As dataframe edata.binary
edata.binary <- as.data.frame(edata.binary)
# Add Mass to binary matrix
edata.binary$Mass <- rownames(edata)
# Add metadata to binary matrix
tmp <- peakObj$e_meta[,c("Mass","bs1_class","ElComposition")] # save to tmp; direct merge isn't working
str(tmp)

'data.frame':   3704 obs. of  3 variables:
 $ Mass         : chr  "200.0564425" "200.092833" "201.004077" "201.0404701" ...
 $ bs1_class    : chr  "Lignin;Amino Sugar" "Protein" "Tannin" "Tannin" ...
 $ ElComposition: chr  "CHON" "CHON" "CHO" "CHO" ...

# Merge 
edata.binary <- merge(edata.binary,tmp,by="Mass",all = T)
# Remove Mass column from df
edata.binary$Mass <- NULL

# Carbon compounds #
carbon.cpds <- aggregate(.~bs1_class*ElComposition,data=edata.binary,FUN=sum)
# Transfor count to relative abundance #
carbon.cpds[,c(-1,-2)] <- decostand(carbon.cpds[,c(-1,-2)],method = "total",MARGIN = 2)*100
# Sum check
colSums(carbon.cpds[,c(-1,-2)])

  H1   H2   H3   H4   H5   H6   H7   H8   H9  H10  H11  H12  H13  H14  H15  H16   C1   C2   C3   C4   C5   C6   C7   C8   C9  C10 
 100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
 C11  C12 CS13  C14  C15  C16 
 100  100  100  100  100  100 

# melt df #
carbon.cpds <- melt(carbon.cpds,variable.name = "Sample_ID")

Using bs1_class, ElComposition as id variables

# Add sample info
carbon.cpds <- merge(carbon.cpds,peakObj$f_data,by="Sample_ID")
head(carbon.cpds)

CHON Subset

tmp <- carbon.cpds[carbon.cpds$ElComposition %in% "CHON",]
tmp

CHOSP Subset

tmp <- carbon.cpds[carbon.cpds$ElComposition %in% "CHOSP",]
tmp