R code for paper Black soldier fly reared on pig manure: bioconversion efficiencies, nutrient reductions and greenhouse gasses

Journal of Waste Management
Alejandro Parodi ^a , Walter J.J. Gerrits ^b, Joop J.A. Van Loon ^c, Imke J.M. De Boer ^a, André J.A. Aarnink ^d, Hannah H.E. Van Zanten ^a

^a Animal Production Systems group, Wageningen University & Research, P.O. Box 338, 6700 AH Wageningen, the Netherlands
^b Animal Nutrition group, Wageningen University & Research, P.O. Box 338, 6700 AH Wageningen, the Netherlands
^c Laboratory of Entomology, Wageningen University & Research, P.O. Box 16, 6700 AA Wageningen, the Netherlands
^d Department of Livestock and Environment, Wageningen University & Research, P.O. Box 338, 6700 AH Wageningen, the Netherlands

Code elaborated by: Alejandro Parodi. alejandro.parodiparodi@wur.nl
https://orcid.org/0000-0003-1351-138X

This rnotebook contains all the steps made for data management, statistical analysis, tables and visualizations.

#Set working directory
setwd("write path here")

Graphical summary of this script

Activate packages

The following packages are be needed to run the code:

Part 1 - Nutrient balance

Step 1 - Data manipulation of datasets fresh_biomass.csv and chemical_analysis.csv

#Import datasets fresh_biomass and chemical_analysis
#Fresh biomass - contains initial amounts (in grams) of inputs and outputs (except condensed water and acid)
fresh_biomass <- read.csv("fresh_biomass.csv") 

#Chemical analysis - contains the chemical results of every input and output
chemical_analysis <- read.csv("chemical_analysis.csv")

#Explanation of columns
#Fresh biomass
#[1] rep - repetition number (R1, R2, R3 and R4)
#[2] treatment - treatment name (With_BSFL and Without)
#[3] chamber - chamber number (7 and 8)
#[4] part - name of the part (Fresh_manure, Starters (7-day old larvae), Larvae (16-day old larvae), Residues Without, Residues With_BSFL)
#[5] value - amount (in grams) of each part


# Chemical analysis
#[1] rep - repetition number (R1, R2, R3 and R4)
#[2] treatment name (With_BSFL and Without)
#[3] chamber - chamber number (7 and 8)
#[4] part - name of the part (Fresh_manure, Starters (7-day old larvae), Larvae (16-day old larvae), Residues Without, Residues With_BSFL)
#[5] DM - Dry matter content of every part (g/kg)
#[6] NH3 - Ammonia (g/kg)
#[7] C - Total carbon (g/kg)
#[8] N - Total N (g/kg)
#[9]P - Phosphorus (g/kg)
#[10]K - Potassium (g/kg)
#[11]E - Energy (kJ/kg)

#Create common column in both datasets to merge them
fresh_biomass$concatenate <- with(fresh_biomass, paste0(trial, "_", treatment, "_", part))
chemical_analysis$concatenate <- with(chemical_analysis, paste0(trial, "_", treatment, "_", part))

#merge both datasets
nutrient_balance <- merge(chemical_analysis, fresh_biomass,  by="concatenate", all=T)

#Remove unnecessary columns
nutrient_balance <- nutrient_balance [, -(13:16), drop=FALSE]

#reorder columnns. This dataset contains the fresh yields (g) of each part (columns 5:6) and next to it the chemical composition in dry weight of each part (columns 8:14)
nutrient_balance <- nutrient_balance[,c(1, 2, 3, 4, 5, 13, 6, 7, 8, 9, 10, 11, 12)]

#calculate the dry yields (g) in each part, and estimate how much of each nutrients is present in each part
#create new column for dryweight values
nutrient_balance$dryweight_grams <- nutrient_balance$value*((nutrient_balance$DM)/1000) 

#create new columns for the overall content of each nutrient (in grams)
nutrient_balance$N_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$N)/1000
nutrient_balance$C_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$C)/1000
nutrient_balance$E_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$ENERGY)/1000
nutrient_balance$K_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$K)/1000
nutrient_balance$P_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$P)/1000
nutrient_balance$NH3_overall <- (nutrient_balance$dryweight_grams*nutrient_balance$NH3)/1000

#select and re-name columns
#Select columns
nutrient_balance <- nutrient_balance %>% select(2:5,14:20)
names(nutrient_balance)[1] <- "trial"
names(nutrient_balance)[2] <- "treatment"
names(nutrient_balance)[3] <- "chamber"
names(nutrient_balance)[4] <- "part"

nutrient_balance

Step2 - Data manipulation for nitrogen emissions

#Import dataset Acid (AC) and condensed water (CW) - contains parameters related to AC and CW to estimate how much nitrogen was lost as gaseous losses. 
acid_condensed_water <- read.csv("acid_condensed_water.csv")

#Columns Acid_condensed water:
#[1] rep - repetition number (R1, R2, R3 and R4)
#[2] treatment name (With_BSFL and Without)
#[3] chamber - chamber number (7 and 8)
#[4] acid_bottleempty_grams - Weight of the empty container that contained the acid (grams)
#[5] acid_bottle_afterperiod_grams - Weight of the container with acid sample at the end of the experimental period (grams)
#[6] acid_gasmeter_valueupfront_m3 - Value of the gas meter before the experimental period (in m3). The flowmeter measured the amount of air (m3) that isflushed into the container with acid during the experimental period
#[7] acid_gasmeter_valueafterperiod_m3 - Value of the gas meter at the end of the experimental period (in m3).
#[8] net_air_flowm3 - net amount of air that flushed (column 7 - column 6)
#[9] overall_air_flow - Overall air flow to the climate respiration chambers (in m3)
#[10] N_acid_units - grams of N (g/kg) in acid
#[11] cw_net_grams -  grams of condensed water in each repetition (grams)
#[12] N_condensedwater_units - grams of N (g/kg) in condensed water

#separate the dataset acid_condensed_water in two sections: one for acid and the other for condensed water
acid <- acid_condensed_water %>% select(1:10)
condensedwater <- acid_condensed_water %>% select(1, 2, 3, 11, 12)

# Acid
#calculate net weight values of acid (grams) and total gas flow
acid$netbottle_grams <- acid$acid_bottle_afterperiod_grams - acid$acid_bottleempty_grams

#calculate the total amount of N in each acid sample
acid$N_sample_grams <- (acid$N_acid_units*acid$netbottle_grams)/1000

#calculate the total grams of N in the whole experimental period period
acid$N_overall <- acid$N_acid_units/1000 * (acid$overall_air_flow*acid$netbottle_grams/acid$net_air_flowm3)
acid$C_overall <- 0
acid$E_overall <- 0
acid$K_overall <- 0
acid$P_overall <- 0
acid$NH3_overall <- 0

#Reorder dataset to bind it to nutritional_overall
acid <- acid %>%  select(1:3, 13:18)
acid$part <- "Acid"
acid$dryweight_grams <- NA
acid <- acid[, c(1, 2, 3, 10, 11, 4, 5, 6, 7, 8, 9)]


#Condensed water
condensedwater$N_overall <- (condensedwater$cw_net_grams * condensedwater$N_condensedwater_units)/1000
condensedwater$C_overall <- 0
condensedwater$E_overall <- 0
condensedwater$K_overall <- 0
condensedwater$P_overall <- 0
condensedwater$NH3_overall <- 0


#format dataset as acid2
condensedwater <- condensedwater %>% select (1, 2, 3, 6:11)
condensedwater$part <- "Condensed_Water"
condensedwater$dryweight_grams <- NA
condensedwater <- condensedwater[, c(1, 2, 3, 10, 11, 4:9)]

#newdataset for total N emissions (sum acid and condensed water)
total_n <- rbind((condensedwater %>%select(1, 2, 4, 6)), (acid %>%select(1, 2, 4, 6)))
total_n <-  ddply(total_n, .(trial, treatment), summarize, 
                  N_g= sum(N_overall, na.rm=TRUE))

total_nitrogen <- rbind(condensedwater, acid)
total_nitrogen

Step3 - Data manipulation for CO2 and CH4 emissions, and heat loss

#Import dataset with gas and heat data obtained from the respiration chamber (CO2, CH4, O2, ENERGY)
chambers <- read.csv("chamberdata_manure.csv")

#Columns chamber data
# Chambers
#[1] chamber - chamber number (7 and 8)  
#[3] date - date of the measurement in format day/month/year
#[4] time - time of the measurement in format hour/minute/second
#[5] Heat - Heat production in KJ/day
#[6] O2 - Oxigen consumption (L/day)
#[7] CO2 - Carbon dioxide production (L/day)
#[8] RQ - Respiratory quotient (obtained by dividing CO2/O2)
#[9] Barometric_pressure_mmHg - Barometric_pressure inside the chamber
#[10] Temperature - Temperature (Celsius) inside the chamber
#[11] Humidity - Humidity (%) inside the chamber
#[12] STPD - Air ventilation (L/min) at standard pressure and temperature
#[13] CH4 - methane (L/day)
#[14] treatment - treatment name (With_BSFL and Without)
#[15] hour - hour of the day in 24h format
#[16] minute - minute at which the measurement was taken (from 0 to 60)
#[17] second - second at which the measurement was taken (from 0 to 60)
#[18] daynumber - day number (from day1 and day10) in which the measurement took place

#Calculate the production of each gas and heat per hour
chambers_summary2 <- ddply(chambers, .(treatment, trial, daynumber, hour, chamber), summarize, 
                           mean_CO2= mean(CO2, na.rm=TRUE), 
                           se_CO2= std.error(CO2, na.rm=TRUE),
                           mean_CH4= mean(CH4, na.rm=TRUE), 
                           se_CH4= std.error(CH4, na.rm=TRUE),
                           mean_O2=mean(O2, na.rm=TRUE),
                           se_O2=std.error(O2, na.rm=TRUE),
                           mean_heat= mean(Heat, na.rm=TRUE), 
                           se_heat= std.error(Heat, na.rm=TRUE),
                           mean_temp= mean(Temperature, na.rm=TRUE), 
                           se_temp= std.error(Temperature, na.rm=TRUE),
                           mean_RQ= mean(RQ, na.rm=TRUE), 
                           se_RQ= std.error(RQ, na.rm=TRUE),
                           mean_STDP= mean(STPD, na.rm=TRUE))

#Format column hour as numeric and set the correct order of the days and replicates

chambers_summary2$hour <- as.numeric(as.character(chambers_summary2$hour))
chambers_summary2$daynumber= factor(chambers_summary2$daynumber, levels =c("day1",
                                                                           "day2",
                                                                           "day3",
                                                                           "day4",
                                                                           "day5",
                                                                           "day6",
                                                                           "day7",
                                                                           "day8",
                                                                           "day9",
                                                                           "day10"))

chambers_summary2$trial= factor(chambers_summary2$trial, levels =c("T1", 
                                                               "T2", 
                                                               "T3", 
                                                               "T4"))
#Edit chamber number
chambers_summary2$chamber <- ifelse(chambers_summary2$chamber=="8", "C8", "C7")


#Format column to bind it to the ones from steps 1 and 2
chambers_summary2$part <- "Emissions/loss"
chambers_summary2$dryweight_grams <- NA
chambers_summary2 <- chambers_summary2[, c(2, 1, 5, 19, 20, 3, 4, 6, 8, 10, 12, 14, 16)]

#All gaseous losses are expressed in l/day. To get the overall production of each gas in grams (and heat in KJ), divide every record by 24 and multiply this value by a factor  to convert it from liters to grams. This factor is the molar mass/22.4
#CO2 -> 44.01/22.4 = 1.964732
#O2 -> 32/22.4 = 1.428571
#CH4 -> 16.06/22.4 = 0.7169643

chambers_summary3 <- ddply(chambers_summary2, .(trial, treatment, chamber, part, dryweight_grams), summarize, 
                           CO2= sum(mean_CO2/24, na.rm=TRUE) * 1.964735,
                           O2= sum(mean_O2/24, na.rm=TRUE) * 1.428571,
                           CH4= sum(mean_CH4/24, na.rm=TRUE)* 0.7169643,
                           HEAT= sum(mean_heat/24, na.rm=TRUE))
total_CO2_CH4_HEAT_O2 <- chambers_summary3
                           
#Calculate the content of C (grams) in CO2 and CH4.
chambers_summary3$N_overall <- 0
chambers_summary3$C_overall <- (chambers_summary3$CO2*0.27272727) + (chambers_summary3$CH4*0.747198) #these factors are used to get the grams of carbon.
chambers_summary3$E_overall <- chambers_summary3$HEAT
chambers_summary3$K_overall <- 0
chambers_summary3$P_overall <- 0
chambers_summary3$NH3_overall <- 0
chambers_summary3 <- chambers_summary3 %>% select (1, 2, 3, 4, 5, 10:15)

chambers_summary3

Step 4 - Bind datasets from step1, step2 and step3

#rbind all bases (fresh, acid, condensed water, chambers_summary)
NB <- rbind(nutrient_balance, acid, condensedwater, chambers_summary3)

#Create new columns (category (Input & Output) and part2 that will be used for the figures
NB$category <- ifelse((NB$part=="Starters" | NB$part=="Fresh_manure"), "Input", "Output")
NB$part2 <- ifelse((NB$part=="Acid" | NB$part=="Condensed_Water" | NB$part=="Aerial_losses"), "Emissions/loss", as.character(NB$part))
#reorder columns
NB <- NB[, c(1, 2, 3, 4, 12, 13, 5:11)]

###Step5 - Prepare dataset to calculate recoveries (sum of inputs and outputs)

#Sum the nutrients in starters and fresh manure to get the total inputs per nutrient (in grams) and format dataset to merge
NB_inputs <- subset(NB, (category=="Input"))
NB_inputs <- gather(NB_inputs, "variable", "value", 7:13)
NB_inputs <- ddply(NB_inputs, .(trial, treatment, chamber, variable), summarize, 
                    total_input= sum(value, na.rm=TRUE))

NB_inputs <- spread(NB_inputs, variable, total_input)
NB_inputs$concatenate <- paste(NB_inputs$trial, "_", NB_inputs$treatment)

# Get the total outputs per nutrient (in grams) and format dataset to merge
NB_outputs <- subset(NB, (category=="Output"))
NB_outputs <- gather(NB_outputs, "variable_output", "value_output", 7:13)
NB_outputs <- ddply(NB_outputs, .(trial, treatment, chamber, part2, variable_output), summarize,
                    sum_parts = sum(value_output, na.rm=T))

NB_outputs <- spread(NB_outputs, variable_output, sum_parts)
NB_outputs$concatenate <- paste(NB_outputs$trial, "_", NB_outputs$treatment)

#Add dry material emissions (calculated as equation 4)
DM_inputs <- NB_inputs %>% select(1, 2, 3, 5)
DM_inputs$concatenate <- paste0(DM_inputs$trial, DM_inputs$treatment, DM_inputs$chamber)

DM_outputs <- ddply(NB_outputs, .(trial, treatment, chamber), summarize,
                    sum_parts = sum(dryweight_grams, na.rm=T))
DM_outputs$concatenate <- paste0(DM_outputs$trial, DM_outputs$treatment, DM_outputs$chamber)

DM_emissions <- merge(DM_inputs, DM_outputs, by="concatenate")
DM_emissions$dm_emissions <- DM_emissions$dryweight_grams- DM_emissions$sum_parts
DM_emissions <- DM_emissions %>% select(2:10)
DM_emissions$concatenate <- paste(DM_emissions$trial.x, "_", DM_emissions$treatment.x)

NB_outputs <- merge(NB_outputs, DM_emissions, by="concatenate")
NB_outputs <- NB_outputs %>% select(1:12, 21)
NB_outputs$dryweight_grams.x <- ifelse(NB_outputs$dryweight_grams.x==0, NB_outputs$dm_emissions, NB_outputs$dryweight_grams.x)
NB_outputs <- NB_outputs %>% select(1:12)


#merge outputs with inputs to calculate the recovery in each part
NB_outputs2 <- merge(NB_outputs, NB_inputs, by="concatenate") #Columns 6:12 contain the outputs and columns 16:22 contain the inputs 
NB_outputs2

Step6 - Calculation of recoveries

#Calculate % of each output
NB_outputs2$DM_per <- (NB_outputs2$dryweight_grams.x/NB_outputs2$dryweight_grams)*100
NB_outputs2$N_per <- (NB_outputs2$N_overall.x/NB_outputs2$N_overall.y)*100
NB_outputs2$C_per <- (NB_outputs2$C_overall.x/NB_outputs2$C_overall.y)*100
NB_outputs2$E_per <- (NB_outputs2$E_overall.x/NB_outputs2$E_overall.y)*100
NB_outputs2$K_per <- (NB_outputs2$K_overall.x/NB_outputs2$K_overall.y)*100
NB_outputs2$P_per <- (NB_outputs2$P_overall.x/NB_outputs2$P_overall.y)*100
NB_outputs2$NH3_per <- (NB_outputs2$NH3_overall.x/NB_outputs2$NH3_overall.y)*100

NB_outputs2 <- NB_outputs2[, c(2:5, 23:28)]
NB_outputs2 <- gather(NB_outputs2, "variable", "value", 5:10)
NB_outputs2

Step7 - Figure S1 - recoveries per trial

#Prepare for Figure S1
NB_outputs2$variable2<- ifelse(NB_outputs2$variable=="DM_per", "Dry matter",
                          ifelse(NB_outputs2$variable=="N_per", "Nitrogen", 
                                 ifelse(NB_outputs2$variable=="C_per", "Carbon",
                                        ifelse(NB_outputs2$variable=="E_per", "Energy",
                                               ifelse(NB_outputs2$variable=="K_per", "Potassium", "Phosphorus")))))

NB_outputs2$variable2 = factor(NB_outputs2$variable2, levels = c("Dry matter",
                                                                 "Nitrogen",
                                                                 "Carbon",
                                                                 "Energy",
                                                                 "Phosphorus",
                                                                 "Potassium"))


NB_outputs2$trial.x=factor(NB_outputs2$trial.x, levels=c("T4",
                                                                 "T3",
                                                                 "T2",
                                                                 "T1"))

NB_outputs2$treatment.x=factor(NB_outputs2$treatment.x, levels=c("With_BSFL",
                                                                 "Without"))

#remove air emissions for P and K
NB_outputs2 <- subset(NB_outputs2, (value>0))

#Supplementary Figure 1

ggplot(NB_outputs2, aes(y= value, x=trial.x, fill=part2,), stat="identity")+
  geom_col() + 
  geom_text(aes(label = (paste0(round(value, digits = 1), "%")), group = part2), position = position_stack(vjust = 0.5), size=2.5)+
  scale_fill_manual(values=c("#FFA500", "#8FBC8F", "#808000"))+
  facet_grid(variable2~treatment.x)+
  labs(y=expression("Recovery (%)"))+ 
  xlab("") +
  theme(strip.background = element_rect(colour="white", fill="white"))+
  theme(strip.text = element_text(size=10, face="bold", hjust=0))+ 
  theme(panel.grid.major  = element_line(colour = "white", size = 0.2))+ 
  theme(panel.grid.minor  = element_line(colour = "white", size = 0.5))+ 
  ggtitle("") + 
  theme(plot.title=element_text(size=13, face="bold", hjust = 1, vjust=0)) + 
  theme(plot.title = element_text(hjust = 0)) + 
  theme(axis.text.x = element_text(size = 10, face = "plain",  angle = 0, hjust = 0.5, vjust=1, colour="black", lineheight = 0)) + 
  theme(strip.text.y = element_text(angle = 0)) + 
  theme(strip.background=element_rect(fill="white")) + 
  guides(fill=guide_legend(title="Part")) +
  theme(axis.title.y = element_text(size=12, face="bold")) + 
  theme(axis.title.x  = element_text(size=14, face="bold"))+ 
  coord_cartesian(ylim=c())+
  theme(legend.position="right")+
  scale_y_continuous(expand = c(0, 0),limits=c(0, 110),breaks = seq(0, 100, by = 25))+
  geom_hline(yintercept=100, linetype="dashed", 
             color = "black", size=0.5)+
  coord_flip()

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#ggsave(path="outputs", "FigureS2.png", width = 22, height = 15, units = "cm", dpi=300)

Step8 - Statistics - Do recovery rates differ between treatments?

statistics2 <- subset(NB_outputs2, (part2!="Larvae"))
statistics2 <- statistics2 %>%select(1, 2, 4, 7, 6)

#Anova
prueba <- statistics2 %>% 
  group_by(part2, variable2) %>%
  group_modify(~ broom::tidy(aov(value~treatment.x + trial.x, data = .x)))

prueba$significance <- ifelse(prueba$p.value<0.05 & prueba$p.value>0.01, "*",
                              ifelse(prueba$p.value<0.01 & prueba$p.value>0.001, "**",
                                     ifelse(prueba$p.value<0.001, "***", "")))

prueba

#write.csv(prueba, file="Anova_recoveries.csv", row.names = F)

Step9 - Figure 1 - Averaged recoveries per nutrient

Figure1 <- NB_outputs2
Figure1 <- ddply(Figure1, .(treatment.x, part2, variable2), summarize, 
                    mean= mean(value, na.rm=TRUE),
                    std_error= std.error(value, na.rm = TRUE))




ggplot(Figure1, aes(y= mean, x=treatment.x, fill=part2,), stat="identity")+
  geom_col(width = 0.85) + 
  geom_text(aes(label = (paste0(round(mean, digits = 0), " ± ", round(std_error, digits = 1))), group = part2), position = position_stack(vjust = 0.5), size=2.5)+
  scale_fill_manual(values=c("#FFA500", "#8FBC8F", "#808000"))+
  labs(y=expression("Recovery (%)"))+ 
  xlab("") +
  facet_grid(.~variable2)+
  theme(strip.background = element_rect(colour="white", fill="white"))+
  theme(strip.text = element_text(size=10, face="bold", hjust=0.5))+ 
  theme(panel.grid.major  = element_line(colour = "white", size = 0.2))+ 
  theme(panel.grid.minor  = element_line(colour = "white", size = 0.5))+ 
  ggtitle("") + 
  theme(plot.title=element_text(size=13, face="bold", hjust = 1, vjust=0)) + 
  theme(plot.title = element_text(hjust = 0)) + 
  theme(axis.text.x = element_text(size = 10, face = "plain",  angle = 90, hjust = 0.5, vjust=1, colour="black", lineheight = 0)) + 
  theme(strip.text.y = element_text(angle = 0)) + 
  theme(strip.background=element_rect(fill="white")) + 
  guides(fill=guide_legend(title="Legend")) +
  theme(axis.title.y = element_text(size=12, face="bold")) + 
  theme(axis.title.x  = element_text(size=14, face="bold"))+ 
  coord_cartesian(ylim=c())+
  theme(legend.position="bottom")+
  scale_y_continuous(expand = c(0, 0),limits=c(0, 110),breaks = seq(0, 100, by = 25))+
  geom_hline(yintercept=100, linetype="dashed", 
             color = "black", size=0.5)

#Save figure. This figure then was editted in Inkscape to add the bars and P.values showing the statistical differences bewtween treatments
#To save:
#ggsave(path="outputs", "Figure1.pdf", width = 22, height = 15, units = "cm", dpi=300)

Part 2 - Ammonia-N balance and nutrient ratios

Step10 - Ammonia balance

#get dataset
ammonia_balance <- NB
ammonia_balance <- subset(ammonia_balance, ((part2=="Fresh_manure" & treatment=="With_BSFL") | part2=="Residues" | part2=="Larvae"))

#select relevant columns
ammonia_balance <- ammonia_balance[, c(1, 2, 6, 8, 13)]

#calculate non nh3 N
ammonia_balance$Non_ammonia_N <- ifelse(is.na(ammonia_balance$NH3_overall), ammonia_balance$N_overall, (ammonia_balance$N_overall-ammonia_balance$NH3_overall))
colnames(ammonia_balance)[4] <- "Total_N"
colnames(ammonia_balance)[5] <- "Ammonia_N"

#save overall N in fresh manure (will be used in a later stage)
manure_n <- ammonia_balance[, c(1:4)]
manure_n <- subset(manure_n, (part2=="Fresh_manure"))

#Re-arrange dataset
ammonia_balance <- ammonia_balance[, c(1, 2, 3, 5, 6)]
ammonia_balance <- gather(ammonia_balance, "type_n", "value", 4:5)

#add  nh3 emissions
emissions_n <- total_n 
emissions_n$part2 <- "Emissions"
emissions_n$type_n <- "Ammonia_N"
colnames(emissions_n)[3] <- "value"
emissions_n <- emissions_n[, c(1, 2, 4, 5, 3)]

#bind two datasets
ammonia_balance <- rbind(ammonia_balance, emissions_n)

#add total n input to calculate percentages
manure_n <- manure_n[, c(1, 4)]
colnames(manure_n) [2] <- "N_input_manure"
ammonia_balance <- merge(ammonia_balance, manure_n, by="trial")

#calculate recovery
ammonia_balance$percentage <- ammonia_balance$value/ammonia_balance$N_input_manure*100
#remove NA
ammonia_balance <- na.omit(ammonia_balance)

ammonia_balance <- ammonia_balance[, c(1:4, 7)]

Step11 - Do ammonia and non-ammonia fractions differ between treatments?

#Anova
statistics_ammonia <- subset(ammonia_balance, (part2!="Fresh_manure" & part2!="Larvae"))

statistics_ammonia2 <- statistics_ammonia %>% 
  group_by(part2, type_n) %>%
  group_modify(~ broom::tidy(aov(percentage~treatment + trial, data = .x)))

statistics_ammonia2

#write.csv(statistics_ammonia2, file="Anova_Nammonia.csv", row.names = F)

Step12 - Figure 2 - Nitrogen balance

figure3 <- ammonia_balance

#New categories for the different fractions
figure3$categories <- ifelse(((ammonia_balance$part2=="Fresh_manure" |  ammonia_balance$part2=="Residues") & (ammonia_balance$type_n=="Ammonia_N")), "Ammonia-N in manure or residues",
                        ifelse(((ammonia_balance$part2=="Fresh_manure" |  ammonia_balance$part2=="Residues") & (ammonia_balance$type_n=="Non_ammonia_N")), "Non ammonia-N in manure or residues",
                               ifelse(ammonia_balance$part2=="Larvae", "N in larval biomass", "Ammonia-N emissions")))

figure3$name2_treatments <- ifelse(figure3$treatment=="With_BSFL" & figure3$part2=="Fresh_manure", "Fresh_manure", 
                   ifelse(figure3$treatment=="With_BSFL" & figure3$part2=="Residues", "With_BSFL",
                          ifelse(figure3$treatment=="With_BSFL" & figure3$part2=="Larvae", "With_BSFL", 
                                  ifelse(figure3$treatment=="With_BSFL" & figure3$part2=="Emissions", "With_BSFL", 
                                          ifelse(figure3$treatment=="Without" & figure3$part2=="Residues", "Without", "Without")))))
                        
figure3_s <- ddply(figure3, .(name2_treatments, categories),
                       summarize,
                       mean=mean(percentage),
                       std_error= std.error(percentage))

#calculate confidence interal for fresh manure
ci_ammonia_manure <- subset(figure3, (part2=="Fresh_manure"))
ci_ammonia_manure <- ci_ammonia_manure[, c(1, 7, 4, 5)]
ci_ammonia_manure <- ci_ammonia_manure %>%
                     group_by(type_n) %>%
                      do(tidy(CI(.$percentage)))

'tidy.numeric' is deprecated.
See help("Deprecated")`data_frame()` was deprecated in tibble 1.1.0.
Please use `tibble()` instead.'tidy.numeric' is deprecated.
See help("Deprecated")

ci_ammonia_manure <- spread(ci_ammonia_manure, names, x)
ci_ammonia_manure$meanlow<- ci_ammonia_manure$mean - ci_ammonia_manure$lower
ci_ammonia_manure$meanup <- ci_ammonia_manure$upper - ci_ammonia_manure$mean

#add column with confidence interval for fresh manure, the rest keep the standard error

figure3_s$variation <- ifelse(figure3_s$name2_treatments=="Fresh_manure" , as.numeric(14), figure3_s$std_error)

figure3_s$categories = factor(figure3_s$categories, 
                                         levels = c("N in larval biomass",
                                                    "Ammonia-N in manure or residues",
                                                    "Ammonia-N emissions",
                                                    "Non ammonia-N in manure or residues"))



ggplot(figure3_s, aes(fill=categories, y=mean, x=name2_treatments)) + 
  geom_bar(position="stack", stat="identity")+
  geom_text(aes(label = (paste0(round(mean, digits = 0), " ± ", round(variation, digits = 0))), group = categories), position = position_stack(vjust = 0.5), size=2.5)+
  theme_classic()+
  ylab("Percentage (%)")+
  xlab("")+
  scale_fill_manual(values = c("#D8E0BB", '#B6CEC7', "#86A3C3",'#7268A6'))+
  theme(legend.position = "right")+
  theme(axis.title.y = element_text(size=12)) + 
  theme(axis.title.x  = element_text(size=12))+
  theme(axis.text.x = element_text(size = 12))+
  theme(axis.text.y = element_text(size = 12))+
  scale_y_continuous(expand = c(0, 0),limits=c(0, 105),breaks = seq(0, 100, by = 25))+
  theme(legend.title=element_blank())

#ggsave(path="outputs", "Figure3.pdf", width = 22, height = 15, units = "cm", dpi=300)

Step13 - Data manipulation to calculate nutrient ratios

#Create new dataset from NB (part 1)
ratio <- NB

#Subset to obtain only the relevant parts
ratio <- subset(ratio, (part=="Fresh_manure" & treatment=="With_BSFL") | part2=="Residues")
#Modify the names of the categories to have "Fresh_manure", "Residues_larvae" and "Residues_manure"
ratio$part2 <- ifelse(ratio$part2=="Residues" & ratio$treatment=="With_BSFL", "Residues_larvae",
                          ifelse(ratio$part2=="Residues" & ratio$treatment=="Without", "Residues_manure", "Fresh_manure"))

#Select relevant columns
ratio <- ratio[, c(1, 6:13)]
#re-name columns
colnames(ratio) <- c("trial", 
                              "Item", 
                              "Dry matter (g)", 
                              "Total nitrogen (g)",
                              "Carbon (g)",
                              "Energy (g)",
                              "Potassium (g)",
                              "Phosphorus (g)",
                              "Ammonia (g)") 

#Remove energy & NH3
ratio <- ratio [, c(1:5,7, 8)]

#Calculate ratios
ratio$C_N <- ratio$`Carbon (g)`/ratio$`Total nitrogen (g)`
ratio$N_P <- ratio$`Total nitrogen (g)`/ratio$`Phosphorus (g)`
ratio$N_K <- ratio$`Total nitrogen (g)`/ratio$`Potassium (g)`
ratio$P_K <- ratio$`Phosphorus (g)`/ratio$`Potassium (g)`
ratio <- ratio[, c(1, 2, 8:11)]

ratio <- gather(ratio, ratio, value, 3:6)
ratio

NA

Step14 - Statistics - Do nutrient ratios differ between treatments?

#Anova
ratio_statistics <- subset(ratio, (Item!="Fresh_manure"))
ratio_statistics <- ratio_statistics%>% 
  group_by(ratio) %>%
  group_modify(~ broom::tidy(aov(value~Item + trial, data = .x)))

ratio_statistics$significance <- ifelse(ratio_statistics$p.value<0.05 & ratio_statistics$p.value>0.01, "*",
                              ifelse(ratio_statistics$p.value<0.01 & ratio_statistics$p.value>0.001, "**",
                                     ifelse(ratio_statistics$p.value<0.001, "***", "")))

ratio_statistics

#write.csv(ratio_statistics, file="Anova_ratios.csv", row.names = F)

Step15- Figure 3 - Averaged nutrient ratios with confidence interval for fresh manure

#Calculate averaged nutrient ratios
ratio_average <- subset(ratio, (Item!="Fresh_manure"))
ratio_average <- ddply(ratio_average, .(Item, ratio),
                       summarize,
                       mean=mean(value, na.rm=T),
                       std_error = std.error(value, na.rm=T))

#confidence interval fresh pig manure
ci_freshmanure <- subset(ratio, (Item=="Fresh_manure"))

ci_freshmanure <- ci_freshmanure %>%
                  group_by(ratio) %>%
                  do(tidy(CI(.$value)))

'tidy.numeric' is deprecated.
See help("Deprecated")'tidy.numeric' is deprecated.
See help("Deprecated")'tidy.numeric' is deprecated.
See help("Deprecated")'tidy.numeric' is deprecated.
See help("Deprecated")

ci_freshmanure <- spread(ci_freshmanure, names, x)

#merge confidence interval in the ratio dataset
ratio_average <- merge(ratio_average, ci_freshmanure, by="ratio")

#Change names of ratios
ratio_average$ratio <- ifelse(ratio_average$ratio=="C_N", "C:N",
                              ifelse(ratio_average$ratio=="N_K", "N:K",
                                     ifelse(ratio_average$ratio=="N_P", "N:P", "P:K")))


#Figure

ggplot(ratio_average, aes(x=Item, y=mean.x)) + 
  geom_point()+ 
  geom_errorbar(aes(min=mean.x-std_error, ymax=mean.x+std_error), width=.3) +
  geom_ribbon(aes(ymin = lower, ymax = upper, group=ratio), fill = "#20639B", alpha=0.6) +
  facet_wrap(~ratio, scales="free_y", ncol = 5) +
  labs(y=expression("Nutrient ratio"))+ 
  xlab("") +
  theme(strip.text = element_text(size=10, face="bold", hjust=0.5))+ 
  theme(panel.grid.major  = element_line(colour = "white", size = 0.2))+ 
  theme(panel.grid.minor  = element_line(colour = "white", size = 0.5))+ 
  ggtitle("") + 
  theme(plot.title=element_text(size=13, face="bold", hjust = 1, vjust=0)) + 
  theme(plot.title = element_text(hjust = 0)) + 
  theme(axis.text.x = element_text(size = 10, face = "plain",  angle = 65, hjust = 0.5, vjust=0.5, colour="black", lineheight = 0)) + 
  theme(strip.text.y = element_text(angle = 0)) + 
  theme(axis.line = element_line(size = 0.5, colour = "black"))+
  theme(strip.background = element_rect(colour="white", fill="white"))+
  theme(panel.background = element_rect(fill = "white", colour = "white"))+
  theme(axis.title.y = element_text(size=12, face="bold")) + 
  theme(axis.title.x  = element_text(size=14, face="bold"))+ 
  coord_cartesian(ylim=c())+
  theme(legend.position="none")

#Save figure. This figure was then editted in Inkscape.
#ggsave(path="outputs", "Figure2.pdf", width = 22, height = 15, units = "cm", dpi=300)

Part3 - Gas emissions and heat loss

Step16 - Gaseous emissions and heat losses (CO2, CH4, heat, oxygen)

#Get CO2, CH4, O2, HEAT and RQ
Gas_average <- chambers_summary2
Gas_average <- Gas_average[, c(1, 2, 3, 6, 7, 8, 9, 10, 11,13)]
Gas_average <- gather(Gas_average, "variable", "value", 6:10)


##import initial dry matter per trial (in grams) to estimate the emissions of gas/hour/kg initial manure
drymanure_g <- subset(NB, ((part2=="Fresh_manure" & treatment=="With_BSFL")))
drymanure_g <- drymanure_g[, c(1, 6, 7)]

#add dry matter manure per trial to calculate emissions per kg of DM initial manure
Gas_average <- merge(Gas_average, drymanure_g, by="trial")

#Calculate emissions (g/hour) per kg of initial DM manure input (except for RQ)
Gas_average$newvalue <- ifelse(Gas_average$variable!="mean_RQ", (((Gas_average$value*1000)/Gas_average$dryweight_grams)/24), Gas_average$value)
Gas_average <- Gas_average[, c(1:6, 10)]

#function to make first letter capital
firstup <- function(x) {
  x <- tolower(x)
  substr(x, 1, 1) <- toupper(substr(x, 1, 1))
  x
}

Gas_average$daynumber <- firstup(Gas_average$daynumber)

Step17 - Ammonia emissions

#Ammonia emissions on time
ammonia <- read.csv("ammonia.csv")

##name of the columns
# Ammonia
#[1] chamber - chamber number (7 and 8)  
#[2] date - date of the measurement in format day/month/year
#[3] time - time of the measurement in format hour/minute/second
#[4] verv_stpd - Air ventilation (L/min) corrected for temperature and pressure
#[5] vol_stp - Volume of the chamber (L) corrected for pressure and temperature
#[6] NH3_prod_lday - Ammonia in L/day
#[7] trial - trial number (T1, T2, T3 and T4)
#[8] treatment - treatment name (With_BSFL and Without)
#[9] hour - hour of the day in 24h format
#[10] minute - minute at which the measurement was taken (from 0 to 60)
#[11] second - second at which the measurement was taken (from 0 to 60)
#[12] daynumber - day number (from day1 and day10) in which the measurement took place

#Create dataset for time series plots (average per hour)
NH3_summary <- ddply(ammonia, .(trial, treatment, chamber, daynumber, hour), summarize, 
                     mean_NH3= mean(NH3_prod_lday, na.rm=TRUE), 
                     se_NH3= std.error(NH3_prod_lday, na.rm=TRUE))


NH3_summary$chamber <- ifelse(NH3_summary$chamber=="7", "C7", "C8")



#calculate emissions per nh3/hour/kg initial dm manure
NH3_summary <- merge(NH3_summary, drymanure_g, by="trial")
NH3_summary$newvalue <- ((NH3_summary$mean_NH3 *1000)/NH3_summary$dryweight_grams)/24
NH3_summary$variable <- "mean_NH3"

#Make the first letter of day capital
##function to make first letter capital
firstup <- function(x) {
  x <- tolower(x)
  substr(x, 1, 1) <- toupper(substr(x, 1, 1))
  x
}

NH3_summary$daynumber <- firstup(NH3_summary$daynumber)

NH3_summary <- NH3_summary[, c(1, 2, 3, 4, 5, 11, 10)] #this dataset will be bind to others to make the figure

Step18 - Bind gaseous with ammonia

#rbind ammonia with gas average
Gasses <- rbind(Gas_average, NH3_summary)

#change names and labels
#new names for figures
Gasses$name2 <- ifelse(Gasses$variable=="mean_CH4", 'CH4 (l/h)',
                       ifelse(Gasses$variable=="mean_CO2", 'CO2 (l/h)',
                              ifelse(Gasses$variable=="mean_heat", "Heat (KJ/h)", 
                                     ifelse(Gasses$variable== "mean_O2", "O2 (l/h)", 
                                            ifelse(Gasses$variable=="mean_RQ", "RQ", "NH3 (l/h)")))))

Step19 - Nitrous oxide - data manipulation

#Import dataset nitrous oxide
nitrous_oxide <- read.csv("nnitrous_oxide.csv")

#Name of columns
#[1] day - day number (from day1 and day10) in which the measurement took place
#[2] trial - trial number (T1, T2, T3 and T4)
#[3] roof - concentration of nitrous oxide (in ppm) outside the chamber
#[4] C7 - concentration of nitrous oxide (in ppm) inside chamber 7
#[5] C8 - concentration of nitrous oxide (in ppm) inside chamber 8


#Calculate the production inside the chamber and keep these columns
nitrous_oxide$Chamber7 <- nitrous_oxide$C7-nitrous_oxide$roof
nitrous_oxide$Chamber8 <- nitrous_oxide$C8-nitrous_oxide$roof
nitrous_oxide <- nitrous_oxide[, c(1, 2, 6, 7)]
nitrous_oxide <- gather(nitrous_oxide,  "chamber", "production", 3:4)

#name of the treatment
nitrous_oxide$treatment <-  ifelse(nitrous_oxide$chamber=="Chamber7" & nitrous_oxide$trial=="T1", "Without", 
                                ifelse(nitrous_oxide$chamber=="Chamber8" & nitrous_oxide$trial=="T1", "With_BSFL",
                                       ifelse(nitrous_oxide$chamber=="Chamber7" & nitrous_oxide$trial=="T2", "With_BSFL",
                                              ifelse(nitrous_oxide$chamber=="Chamber8" & nitrous_oxide$trial=="T2", "Without",
                                                     ifelse(nitrous_oxide$chamber=="Chamber7" & nitrous_oxide$trial=="T3", "Without",
                                                            ifelse(nitrous_oxide$chamber=="Chamber8" & nitrous_oxide$trial=="T3", "With_BSFL",
                                                                   ifelse(nitrous_oxide$chamber=="Chamber7" & nitrous_oxide$trial=="T4", "With_BSFL", "Without")))))))

#To calculate the production of N2O, we need to calculate the average ppm concentrations between two consecutive measurements
nitrous_oxide$con <- paste0(nitrous_oxide$trial, nitrous_oxide$treatment)
nitrous_oxide$value <- ifelse(nitrous_oxide$con== lead(nitrous_oxide$con), ((nitrous_oxide$production + lead(nitrous_oxide$production))/2), NA)
nitrous_oxide$daynumber <- ifelse(nitrous_oxide$con== lead(nitrous_oxide$con), paste0(nitrous_oxide$day, "-", lead(nitrous_oxide$day)), NA)
nitrous_oxide <- nitrous_oxide[, c(2, 5, 7, 8)]
nitrous_oxide <- na.omit(nitrous_oxide)


#Get ventilation per day
chambers_2 <- chambers
chambers_2$hour <- as.numeric(chambers_2$hour)

#ventilation is in l/min so to convert it to l/hour -> STPD*60
chambers_n2o <- ddply(chambers_2, .(trial, treatment, daynumber, hour), summarize, 
                      mean_ventilation = mean(STPD*60, na.rm=TRUE))

chambers_n2o$hour <- as.numeric(as.character(chambers_n2o$hour))

#new day classification based on the time of the day in which air samples were collected for N2O analysis
chambers_n2o$new <- ifelse((chambers_n2o$daynumber=="day1")|(chambers_n2o$daynumber=="day2" & chambers_n2o$hour<=8), "Day1-Day2", 
                         ifelse((chambers_n2o$daynumber=="day2" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day3" & chambers_n2o$hour<=8), "Day2-Day3",
                                ifelse((chambers_n2o$daynumber=="day3" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day4" & chambers_n2o$hour<=8), "Day3-Day4", 
                                       ifelse((chambers_n2o$daynumber=="day4" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day5" & chambers_n2o$hour<=8), "Day4-Day5",
                                              ifelse((chambers_n2o$daynumber=="day5" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day6" & chambers_n2o$hour<=8), "Day5-Day6",
                                                     ifelse((chambers_n2o$daynumber=="day6" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day7" & chambers_n2o$hour<=8), "Day6-Day7",
                                                            ifelse((chambers_n2o$daynumber=="day7" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day8" & chambers_n2o$hour<=8), "Day7-Day8",
                                                                   ifelse((chambers_n2o$daynumber=="day8" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day9" & chambers_n2o$hour<=8), "Day8-Day9",
                                                                          ifelse((chambers_n2o$daynumber=="day9" & chambers_n2o$hour>=9)|(chambers_n2o$daynumber=="day10" & chambers_n2o$hour<=8), "Day9-Day10", "nada"
                                                                          )))))))))
#calculate the total ventilation liters per new day number
chambers_n2o <- ddply(chambers_n2o, .(trial, treatment, new), summarize, 
                      total_stpd_vent = sum(mean_ventilation, na.rm=TRUE))


# merge both datasets (n2o production and ventilation)
chambers_n2o$concatenate <- paste0(chambers_n2o$trial, chambers_n2o$treatment, chambers_n2o$new)
nitrous_oxide$concatenate <- paste0(nitrous_oxide$trial, nitrous_oxide$treatment, nitrous_oxide$daynumber)
nitrous_oxide <- merge(nitrous_oxide, chambers_n2o, by="concatenate")

#apply equation to calculate production in mg per day
# 44 is the molar mass of n2o
# 22.4 is the molar volume of a gas
# *1000 is to convert it to miligrams
nitrous_oxide$n2o_mg <- ((nitrous_oxide$value*10^(-6))*nitrous_oxide$total_stpd_vent*(44/22.4))*1000

nitrous_oxide <- nitrous_oxide[, c(2, 3, 5, 10 )]
names(nitrous_oxide) <- c("trial", "treatment", "day", "n2o_mg")


#merge with drymatter manure to estimate n2o production per kg of dry matter manure
nitrous_oxide <- merge(nitrous_oxide, drymanure_g, by="trial")
nitrous_oxide$n2o_mg_perkg <- (1000*nitrous_oxide$n2o_mg)/nitrous_oxide$dryweight_grams
nitrous_oxide

Step20 - Statistics - Is N2O production different from zero?

#Test if N2O differs from zero
#t-test
t_n2o <- nitrous_oxide %>%
  group_by(treatment, day) %>%
  do(tidy(t.test(.$n2o_mg, alternative = "two.sided")))
t_n2o

#write.csv(t_n2o, "T_test_n2o.csv", row.names = F)

Step21 - Supplementary Figure 3 panel A

#Order factors for Supplementary Figure
Gasses$name2 <- factor(Gasses$name2, levels = c("CO2 (l/h)",
                                                            "CH4 (l/h)",
                                                            "Heat (KJ/h)",
                                                            "NH3 (l/h)",
                                                            "O2 (l/h)",
                                                            "RQ"))
Gasses$daynumber<- factor(Gasses$daynumber, levels = c("Day1", 
                                                                   "Day2", 
                                                                   "Day3", 
                                                                   "Day4", 
                                                                   "Day5", 
                                                                   "Day6", 
                                                                   "Day7", 
                                                                   "Day8", 
                                                                   "Day9", 
                                                                   "Day10"))

 ggplot(data=Gasses, aes(x=as.numeric(hour), y=newvalue, group=trial))+
  geom_line(size=1, aes(linetype=trial, col=trial))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(name2~treatment+daynumber, scales = "free_y")+
  theme_bw()+
  ylab("")+
  scale_x_continuous(name="Time of the day (hours)", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  ggtitle("")+
  theme(legend.position = "bottom")+
  theme(axis.title.y = element_text(size=12)) + 
  theme(axis.title.x  = element_text(size=12))+
  scale_linetype_manual(values=c("solid", "dotted", "twodash", "longdash"))

 
 #ggsave(path="outputs", file="FigureS3_A.pdf", width = 25, height = 20, units="cm", dpi=300) #saves g

Step21 Supplementary Figure 3 - panel (B)

#Supplementary figure N2O
ggplot(data=nitrous_oxide, aes(x=day, y=n2o_mg_perkg, group=trial))+
  geom_line(size=1, aes(linetype=trial, col=trial))+
  geom_point()+
  facet_wrap(~treatment)+
  theme(axis.text = element_text(angle = 0, hjust = 0.5, size=10)) + 
  ylab("N2O (mg)/day/kg DM manure")+
  xlab("")+
  theme_bw()+
  theme(panel.spacing = unit(0.1, "lines"))+
  ggtitle("")+
  theme(legend.position = "right")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=4))+
  scale_linetype_manual(values=c("solid", "dotted", "twodash", "longdash"))+
  scale_x_discrete(labels=function(x){sub("-", "\n", x)})

#ggsave(path="outputs", file="FigureS3_B.pdf", width = 25, height = 10, units="cm", dpi=300) #saves g

Step22 - Figure 4

#Calculate average and std_error per day per treatment (for all gasses except N2O)
figure4 <- ddply(Gasses, .(treatment, daynumber, hour, name2), summarize, 
                           average= mean(newvalue, na.rm=TRUE),
                           std_error= std.error(newvalue, na.rm=TRUE))


#order days for figure
figure4$daynumber<- factor(figure4$daynumber, levels = c("Day1", 
                                                                   "Day2", 
                                                                   "Day3", 
                                                                   "Day4", 
                                                                   "Day5", 
                                                                   "Day6", 
                                                                   "Day7", 
                                                                   "Day8", 
                                                                   "Day9", 
                                                                   "Day10"))
#subset to make independent figures of each
CO2 <- subset(figure4, (name2=="CO2 (l/h)"))
CH4 <- subset(figure4, (name2=="CH4 (l/h)"))
HEAT <- subset(figure4, (name2=="Heat (KJ/h)"))
NH3 <- subset(figure4, (name2=="NH3 (l/h)"))
RQ <- subset(figure4, (name2=="RQ"))

#Figure CO2
CO2_fig <- ggplot(data=CO2, aes(x=as.numeric(hour), y=average, col= treatment, group=treatment, fill=treatment))+
  geom_line(size=1, aes(linetype=treatment))+
  geom_ribbon(aes(ymin=(average-std_error), ymax=(average+std_error), linetype=treatment), alpha=0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(~daynumber, scales = "free_y")+
  theme_bw()+
  ylab(expression(CO[2]~(L/hour)))+
  scale_x_continuous(name="", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  theme(strip.background =element_rect(fill="white"))+
  scale_colour_manual(values = c("#5cb85c",'#E69F00'))+
  scale_fill_manual(values = c("#5cb85c",'#E69F00'))+
  ggtitle("")+
  theme(legend.position = "none")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=10))+
  scale_linetype_manual(values=c("solid", "dotted"))

#Figure CH4
CH4_fig <- ggplot(data=CH4, aes(x=as.numeric(hour), y=average*1000, col= treatment, group=treatment, fill=treatment))+
  geom_line(size=1, aes(linetype=treatment))+
  geom_ribbon(aes(ymin=(average*1000-std_error*1000), ymax=(average*1000+std_error*1000), linetype=treatment), alpha=0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(~daynumber, scales = "free_y")+
  theme_bw()+
  ylab(expression(CH[4]~(mL/hour)))+
  scale_x_continuous(name="Time of the day (hours)", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  theme(strip.background =element_rect(fill="white"))+
  scale_colour_manual(values = c("#5cb85c",'#E69F00'))+
  scale_fill_manual(values = c("#5cb85c",'#E69F00'))+
  ggtitle("")+
  theme(legend.position = "none")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=10))+
  scale_linetype_manual(values=c("solid", "dotted"))

#Figure NH3
NH3_fig <- ggplot(data=NH3, aes(x=as.numeric(hour), y=average*1000, col= treatment, group=treatment, fill=treatment))+
  geom_line(size=1, aes(linetype=treatment))+
  geom_ribbon(aes(ymin=(average*1000-std_error*1000), ymax=(average*1000+std_error*1000), linetype=treatment), alpha=0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(~daynumber, scales = "free_y")+
  theme_bw()+
  ylab(expression(NH[3]~(mL/hour)))+
  scale_x_continuous(name="Time of the day (hours)", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  theme(strip.background =element_rect(fill="white"))+
  scale_colour_manual(values = c("#5cb85c",'#E69F00'))+
  scale_fill_manual(values = c("#5cb85c",'#E69F00'))+
  ggtitle("")+
  theme(legend.position = "none")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=10))+
  scale_linetype_manual(values=c("solid", "dotted"))


#Figure NH3
RQ_fig <- ggplot(data=RQ, aes(x=as.numeric(hour), y=average, col= treatment, group=treatment, fill=treatment))+
  geom_line(size=1, aes(linetype=treatment))+
  geom_ribbon(aes(ymin=(average-std_error), ymax=(average+std_error), linetype=treatment), alpha=0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(~daynumber, scales = "free_y")+
  theme_bw()+
  ylab("Respiratory quotient")+
  scale_x_continuous(name="Time of the day (hours)", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  theme(strip.background =element_rect(fill="white"))+
  scale_colour_manual(values = c("#5cb85c",'#E69F00'))+
  scale_fill_manual(values = c("#5cb85c",'#E69F00'))+
  ggtitle("")+
  theme(legend.position = "none")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=10))+
  scale_linetype_manual(values=c("solid", "dotted"))


#Figure N2O
#Given that we didn't take hourly measurements of N2O (only once per day), we need to make some data manipulation in order to
#get the same figure as the other gases.

#Extract same format of dataset from the dataset gas_average and adapt it for N2O
fig_n2o <- Gas_average
fig_n2o <- subset(fig_n2o, (variable=="mean_CH4"))
fig_n2o <- fig_n2o[, c(1, 2, 4, 5)]
fig_n2o$variable <- "mean_N2O"
fig_n2o$n2o_mg_perkg <- NA
fig_n2o <- subset(fig_n2o, (hour!="9"))
fig_n2o <- subset(fig_n2o, ((daynumber!="day1") & (hour!=17)))
fig_n2o <- subset(fig_n2o, ((daynumber!="day10") & (hour!=8)))


#Modify the name nomenclature
nitrous_oxide_modified <- nitrous_oxide
nitrous_oxide_modified$daynumber <- ifelse(nitrous_oxide_modified$day=="Day1-Day2", "Day2",
                          ifelse(nitrous_oxide_modified$day=="Day2-Day3", "Day3",
                                 ifelse(nitrous_oxide_modified$day=="Day3-Day4", "Day4",
                                        ifelse(nitrous_oxide_modified$day=="Day4-Day5", "Day5",
                                               ifelse(nitrous_oxide_modified$day=="Day5-Day6", "Day6",
                                                      ifelse(nitrous_oxide_modified$day=="Day6-Day7", "Day7",
                                                             ifelse(nitrous_oxide_modified$day=="Day7-Day8", "Day8",
                                                                    ifelse(nitrous_oxide_modified$day=="Day8-Day9", "Day9", "Day10"))))))))

#Indicate at what hour of the day N2O measurements were taken
nitrous_oxide_modified$hour <- ifelse(nitrous_oxide_modified$daynumber=="Day10", "8", "9")
nitrous_oxide_modified2 <- subset(nitrous_oxide_modified, (day=="Day1-Day2"))
nitrous_oxide_modified2$daynumber <- "Day1"
nitrous_oxide_modified2$hour <- 17
nitrous_oxide_modified2$n2o_mg_perkg <- 0
nitrous_oxide_modified <- rbind(nitrous_oxide_modified, nitrous_oxide_modified2)
nitrous_oxide_modified$variable <- "mean_N2O"
nitrous_oxide_modified <- nitrous_oxide_modified[, c(1, 2, 8, 9, 10, 7)]


#Bind the dataset with the structure of days and hours, with the one that contain the n2o measurements. 
fig_n2o <- rbind(fig_n2o, nitrous_oxide_modified)
names(fig_n2o)[5] <- "average"
fig_n2o$hour <- as.numeric(fig_n2o$hour)

#Make figure
fig_n2o$daynumber <- factor(fig_n2o$daynumber, levels = c("Day1",
                                                    "Day2",
                                                    "Day3",
                                                    "Day4",
                                                    "Day5",
                                                    "Day6",
                                                    "Day7",
                                                    "Day8",
                                                    "Day9",
                                                    "Day10"))


N2O_fig <- ggplot(data=fig_n2o, aes(x=as.numeric(hour), y=n2o_mg_perkg, fill=treatment, linetype=treatment))+
  geom_boxplot(width=8)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  facet_grid(~daynumber, scales = "free_y")+
  theme_bw()+
  ylab(expression(N[2]~O~(mL/day)))+
  scale_x_continuous(name="Time of the day (hours)", expand = c(0, 0), breaks = c(0, 12, 0))+
  theme(panel.spacing = unit(0.1, "lines"))+
  theme(strip.background =element_rect(fill="white"))+
  scale_colour_manual(values = c("#5cb85c",'#E69F00'))+
  scale_fill_manual(values = c("#5cb85c",'#E69F00'))+
  ggtitle("")+
  theme(legend.position = "none")+
  theme(axis.title.y = element_text(size=10)) + 
  theme(axis.title.x  = element_text(size=10))+
  scale_linetype_manual(values=c("solid", "dotted"))


allgraph <- arrangeGrob(CO2_fig, CH4_fig, NH3_fig, N2O_fig, RQ_fig , nrow=3) #generates g

Removed 1456 rows containing non-finite values (stat_boxplot).

#ggsave(path="outputs", file="Figure4.pdf", width = 25, height = 20, units="cm", dpi=300, allgraph) #saves g

knitr::include_graphics("Figure4.png")

Step 22 - Table 2 - emissions per kg of dry manure and per kg of fresh and dry larvae

#Get overall emissions (CO2, CH4, Energy, Oxygen, N and N2O) per replicate, per kg of manure, per kg of fresh and dry larvae
#Calculate N2O in miligrams

nitrous_oxide_mg <- ddply(nitrous_oxide, .(trial,treatment), summarize, 
                 total = sum(n2o_mg))
nitrous_oxide_mg $item <- "N2O (mg)"
names(nitrous_oxide_mg )[3] <- "value"
nitrous_oxide_mg <- nitrous_oxide_mg[, c(1, 2, 4, 3)]


#get total N (sum condensed water + acid)
total_n$item <- "N (g)"
names(total_n)[3] <- "value"
total_n<- total_n[, c(1,2,4,3)]


#get CO2, CH4, Oxygen and Heat
total_CO2_CH4_HEAT_O2

total_CO2_CH4_HEAT_O2 <- total_CO2_CH4_HEAT_O2[, c(1, 2, 6:9)]
names(total_CO2_CH4_HEAT_O2)[3] <- "CO2 (g)"
names(total_CO2_CH4_HEAT_O2)[4] <- "O2 (g)"
names(total_CO2_CH4_HEAT_O2)[5] <- "CH4 (g)"
names(total_CO2_CH4_HEAT_O2)[6] <- "Heat (KJ)"
total_CO2_CH4_HEAT_O2 <- gather(total_CO2_CH4_HEAT_O2, "item", "value", 3:6)


##bind datasets (nitrogen, chambers_summary and nitrous oxide)
Air_emissions <- rbind(total_n, total_CO2_CH4_HEAT_O2, nitrous_oxide_mg)


### Calculate emissions expressed in different units (per kg of dry manure, per kg of fresh larvae, per kg of dry larvae)

#Per kg of fresh larvae
#Add dry matter manure (grams)
Air_emissions <- merge(Air_emissions, drymanure_g, by="trial")
Air_emissions$kg_dry_manure <- (Air_emissions$value*1000)/Air_emissions$dryweight_grams
names(Air_emissions)[2] <- "treatment"


#Per kg of fresh larvae
#Fresh larvae yield per trial
f_l_yield <- subset(fresh_biomass, (part=="Larvae"))
f_l_yield <- f_l_yield[, c(1, 2, 5)]
f_l_yield$treatment <- "Fresh larvae yield (g)"
Air_emissions <- merge(Air_emissions, f_l_yield, by="trial")
Air_emissions$kgfreshlarvae <- (Air_emissions$value.x*1000)/Air_emissions$value.y
Air_emissions <- Air_emissions[, c(1, 2, 3, 4, 7, 10)]

#Per kg of dry larvae
#dry larvae yield per trial
d_l_yield <- subset(NB_outputs, (part2=="Larvae"))
d_l_yield <- d_l_yield[, c(2, 7)]
d_l_yield$treatment <- "Dry larvae yield (g)"
colnames(d_l_yield)[2] <- "value"
Air_emissions <- merge(Air_emissions, d_l_yield, by="trial")
Air_emissions$kgdrylarvae <-  (Air_emissions$value.x*1000)/Air_emissions$value
Air_emissions <- Air_emissions[, c(1:3, 5, 6, 9)]


#Remove the values per kg of fresh and dry larvae for the treatment Without (as there was no larvae)
Air_emissions <- gather(Air_emissions, "variable", "value", 4:6)
Air_emissions$concatenate <- paste0(Air_emissions$treat, Air_emissions$variable)
Air_emissions <- subset(Air_emissions, (concatenate!="Withoutkgfreshlarvae" & concatenate!= "Withoutkgdrylarvae")) #This contains the emissions per trial. This dataset will be used in step 24
names(Air_emissions)[2] <- "treatment"
#Use air emissions table to calculate averages (of all trials)

#Create supplementary table 2
TableS2 <- Air_emissions %>% select(1, 3, 4, 5)
TableS2 <- subset(TableS2, (variable!="kg_dry_manure"))
names(TableS2) <- c("Trial", "Item", "Variable", "Value")
TableS2 <- spread(TableS2, Variable, Value)

write.csv(TableS2, file = "TableS2.csv", row.names = F)

Table4 <- ddply(Air_emissions, .(treatment, item,  variable), summarize, 
                   mean = mean(value),
                   std_error =std.error(value))

#This will be done with standard deviation  to compare our emissions with other studies
Table_emissions_perkg_larvae <- ddply(Air_emissions, .(treatment, item,  variable), summarize, 
                   mean = mean(value),
                   std_dev =sd(value))

Table4$meanstderror <- paste0(round(Table4$mean, digits = 2), "±", round(Table4$std_error, digits = 2))
Table4 <- Table4[, c(1, 2, 3, 6)]
Table4 <- spread(Table4, treatment, meanstderror)

Table_emissions_perkg_larvae$meanstderror <- paste0(round(Table_emissions_perkg_larvae$mean, digits = 2), "±", round(Table_emissions_perkg_larvae$std_dev, digits = 2))
Table_emissions_perkg_larvae <- Table_emissions_perkg_larvae[, c(1, 2, 3, 6)]
Table_emissions_perkg_larvae <- spread(Table_emissions_perkg_larvae, treatment, meanstderror)

#Emissions per kg of fresh and dry larvae
em_dry_fresh <- subset(Table_emissions_perkg_larvae, variable=="kgdrylarvae" | variable=="kgfreshlarvae")
em_dry_fresh <- spread(em_dry_fresh, variable, With_BSFL)

#emissions per kg of dm manure
Table4 <- na.omit(Table4)
Table4 <- Table4[, c(1, 4, 3)]


#Calculate emissions in CO2 eq (combination of CH4 and N2O) to add it to Table4
CO2eq <- Air_emissions
CO2eq$CO2eq <- ifelse(CO2eq$item=="N2O (mg)", (((CO2eq$value)/1000)*298), 
                      ifelse(CO2eq$item=="CH4 (g)", (((CO2eq$value))*34), NA))

CO2eq <- na.omit(CO2eq)
CO2eq <- subset(CO2eq, (variable=="kg_dry_manure" | variable=="kgdrylarvae"))

CO2eq <-  ddply(CO2eq, .(trial, treatment, concatenate), summarize, 
                sum = sum(CO2eq))

CO2eq_mean <- ddply(CO2eq, .(treatment, concatenate), summarize, 
                    mean = mean(sum),
                    std_error = sd(sum))

CO2eq_mean$meanstderror <- paste0(round(CO2eq_mean$mean, digits = 1), "±", round(CO2eq_mean$std_error, digits = 1))#here you can see CO2eq per kg DM larvae 343.99±43.14 grams
CO2eq_mean <- CO2eq_mean[c(1, 3),]
CO2eq_mean <- CO2eq_mean[, c(1, 5)]


CO2eq_mean <- spread(CO2eq_mean, treatment, meanstderror)
CO2eq_mean$item <- "CO2 eq (g)"
CO2eq_mean <- CO2eq_mean[, c(3, 2, 1)]


#Finish table 4 (bind Table 4 with CO2eq_mean)
Table4 <- rbind(Table4, CO2eq_mean)
Table4

Step 24 - Statistics do overall emissions differ per treatment?

#Merge air emissions with CO2_eq before doing the multiple regression
#Mofify dataframe CO2eq to bind it to Air emissions
CO2eq <- subset(CO2eq, (concatenate!="With_BSFLkgdrylarvae")) 
CO2eq$item <- "CO2eq (g)"
CO2eq$variable <- "kg_dry_manure"
names(CO2eq)[4] <- "value"
CO2eq <- CO2eq[, c(1, 2, 5, 6, 4, 3)]

#Bind Air emissions and CO2eq and format dataset for statistical test
Air_emissions2 <- rbind(Air_emissions, CO2eq)
Air_emissions2 <- Air_emissions2[, c(1:5)] #remove column concatenate
#Select only the values expressed per kg of dry manure
Air_emissions2 <- subset(Air_emissions2, (variable=="kg_dry_manure"))

#Multiple regression

stats_gases <- Air_emissions2%>% 
  group_by(item) %>%
  group_modify(~ broom::tidy(aov(value~treatment + trial, data = .x)))

stats_gases$significance <- ifelse(stats_gases$p.value<0.05 & stats_gases$p.value>0.01, "*",
                              ifelse(stats_gases$p.value<0.01 & stats_gases$p.value>0.001, "**",
                                     ifelse(stats_gases$p.value<0.001, "***", "")))


stats_gases

#write.csv(stats_gases, file="ANOVA_Table5.csv", row.names = F)

Part 5 - Bioconversion parameters

Step 25 - Individual starter larvae size per batch (For Table S4)

#Estimate individual size of starter larvae per batch
#Weight (grams) group of 50 starters
starters_batch1 <- c(0.67,0.74,0.65,0.61,0.62,0.61,0.62,0.63,0.61,0.65)
starters_batch2 <- c(0.79,0.82,0.76,0.82,0.7,0.78,0.82,0.76,0.74,0.77)
starters_batch3 <- c(0.34,0.36,0.4,0.36,0.38,0.31,0.43,0.38,0.38,0.36)
starters_batch4 <- c(0.29,0.28,0.29,0.36,0.27,0.27,0.26,0.31,0.25,0.29)

batch_weight <- c(starters_batch1, starters_batch2, starters_batch3, starters_batch4)
trial <- rep(c("T1", "T2", "T3", "T4"), each=10, len=40)

#Create dataset starters_weight
starters_weight <- data.frame(batch_weight, trial)

#Calculate weight per individual larva (mg)
starters_weight$individual_weight <- (starters_weight$batch_weight/50)*1000

#Mean and standard error
starters1<- starters_weight %>%
  group_by(trial) %>%
  summarise(
    mean = mean(individual_weight),
    std_error = std.error(individual_weight))

starters1

Step 26 - Bioconversion parameters

#######Supplementary Table with detailed data of the bioconversion process.
#Fresh manure provided
trial <- c("T1", "T2", "T3", "T4")
fresh_m_provided<- c(20814.06, 16649.35, 18404.91, 18407.05)
fresh_manure <- data.frame(trial, fresh_m_provided)
fresh_manure$treatment <- "Fresh manure provided (g)"
fresh_manure <- fresh_manure[, c(1, 3, 2)]
colnames(fresh_manure)[3] <- "value"

#manure dry matter
dm_manure <- c(21.35, 25.39, 23.21, 25.65)
dry_manure <- data.frame(trial, dm_manure)
dry_manure$treatment <- "Manure dry matter (%)"
dry_manure <- dry_manure[, c(1, 3, 2)]
colnames(dry_manure)[3] <- "value"

#dry manure provided
drymanure <- merge(fresh_manure, dry_manure, by="trial")
drymanure$value <- drymanure$value.x*(drymanure$value.y/100)
drymanure <- drymanure[, c(1, 2, 6)]
colnames(drymanure)[2] <- "treatment"
drymanure$treatment <- "Dry manure provided (g)"


#Fresh larvae yield per trial
f_l_yield <- subset(fresh_biomass, (part=="Larvae"))
f_l_yield <- f_l_yield[, c(1, 2, 5)]
f_l_yield$treatment <- "Fresh larvae yield (g)"

#dry larvae yield per trial
d_l_yield <- subset(NB_outputs, (part2=="Larvae"))
d_l_yield <- d_l_yield[, c(2, 7)]
d_l_yield$treatment <- "Dry larvae yield (g)"
d_l_yield <- d_l_yield[, c(1, 3, 2)]
colnames(d_l_yield)[2] <- "treatment"
colnames(d_l_yield)[3] <- "value"

#DM percentage larvae
dm_l <- subset(chemical_analysis, (part=="Larvae"))
dm_l <- dm_l[, c(1, 2, 5)]
dm_l$treatment <- "Larvae dry matter content(%)"
colnames(dm_l)[3] <- "value"
dm_l$value <- dm_l$value/10

#fresh larvae gain
f_lg <- subset(fresh_biomass, (part=="Starters"))
f_lg <- f_lg[, c(1, 4, 5)]
f_lg$part <- "Fresh larvae gain weight (g)"
colnames(f_lg)[2] <- "treatment"
f_lg <- merge(f_lg, f_l_yield, by="trial")
f_lg$value <- f_lg$value.y-f_lg$value.x
f_lg <- f_lg[, c(1, 4, 6)]
colnames(f_lg) [2] <- "treatment"
f_lg$treatment<- "Fresh larvae gain weight (g)"

#dry larvae gain
d_lg <- subset(fresh_biomass, (part=="Starters"))
d_lg <- d_lg[, c(1, 4, 5)]
startersdm <- subset(chemical_analysis, (part=="Starters"))
startersdm <- subset(startersdm, (treatment=="With_BSFL"))
startersdm <- startersdm[, c(1, 2, 5)]
d_lg <- merge(d_lg, startersdm, by="trial")
d_lg$starters_dry<- d_lg$value*(d_lg$DM/1000)
d_lg <- d_lg[, c(1, 4, 6)]
d_lg <- merge(d_lg, d_l_yield, by="trial")
d_lg$drygain <- d_lg$value-d_lg$starters_dry
d_lg <- d_lg[, c(1, 2, 6)]
colnames(d_lg)[2] <- "treatment"
d_lg$treatment <- "Dry larvae gain weight (g)"
colnames(d_lg)[3] <- "value"


#kg of fresh larvae per kg of fresh
fly_xfreshmanure <- merge(fresh_manure, f_l_yield, by="trial")
fly_xfreshmanure$value2 <-(fly_xfreshmanure$value.y*1000)/fly_xfreshmanure$value.x
fly_xfreshmanure <- fly_xfreshmanure[, c(1, 2, 6)]
colnames(fly_xfreshmanure)[2] <- "treatment"
fly_xfreshmanure$treatment <- "Fresh larvae yield (g) per kg of fresh manure"
colnames(fly_xfreshmanure)[3] <- "value"  

#kg dry larvae per kf of dry manure
fly_xdrymanure <- merge(drymanure, d_l_yield, by="trial")
fly_xdrymanure$value2 <-(fly_xdrymanure$value.y*1000)/fly_xdrymanure$value.x
fly_xdrymanure <- fly_xdrymanure[, c(1, 2, 6)]
colnames(fly_xdrymanure)[2] <- "treatment"
fly_xdrymanure$treatment <- "Dry larvae yield (g) per kg of dry manure"
colnames(fly_xdrymanure)[3] <- "value"  


#Individual mature larvae weight
maturelarvae_trial1 <- c(4.24, 4.29, 4.12, 4.26, 3.45, 4.12, 4.27, 3.4, 4.18, 4.14)
maturelarvae_trial2 <- c(4.96, 5.07, 4.87, 4.69, 4.9, 4.7, 5.07, 5.07, 4.76, 4.74)
maturelarvae_trial3 <- c(3.94, 3.91, 4.03, 3.98, 3.98, 3.84, 4.12, 3.76, 4, 3.98)
maturelarvae_trial4 <- c(3.94, 3.91, 4.03, 3.98, 3.98, 3.84, 4.12, 3.76, 4, 3.98)

weight_mature_trial <- c(maturelarvae_trial1, maturelarvae_trial2, maturelarvae_trial3, maturelarvae_trial4)
trial_name <- rep(c("trial1", "trial2", "trial3", "trial4"), each=10, len=40)

#df weight mature larvae
maturelarvae_weight <- data.frame(trial_name, weight_mature_trial)
maturelarvae_weight$value <- (maturelarvae_weight$weight_mature_trial/50)*1000
colnames(maturelarvae_weight)[1] <- "trial"
colnames(maturelarvae_weight)[2] <- "treatment"
maturelarvae_weight$treatment <- "Final larvae individual weight (mg)"

#Weight starter larvae

starterlarvae_weight <- starters1
starterlarvae_weight$treatment <- "Starter larvae individual weight (mg)"
names(starterlarvae_weight)[2] <- "value"
starterlarvae_weight <- starterlarvae_weight[, c(1, 4, 2)]

#Join all datasets and create one with average and standard error
table_bioconversion_details <- rbind(fresh_manure, dry_manure, drymanure, f_l_yield, d_l_yield, dm_l, f_lg, d_lg, fly_xfreshmanure, fly_xdrymanure, maturelarvae_weight, starterlarvae_weight)
average_table_bioconversion_details <- ddply(table_bioconversion_details, .(treatment),
                                            summarize,
                                            mean=mean(value, na.rm=T),
                                            std_error = std.error(value, na.rm=T))


average_table_bioconversion_details$value2 <- paste0(round(average_table_bioconversion_details$mean, digits = 0), " ± ", 
                                                     (round(average_table_bioconversion_details$std_error, digits = 1)))


average_table_bioconversion_details <- average_table_bioconversion_details[, c(1, 4)]
colnames(average_table_bioconversion_details)[1] <- "Parameter"
colnames(average_table_bioconversion_details)[2] <- "Value"

average_table_bioconversion_details <- average_table_bioconversion_details[c(8, 4, 6, 2, 10, 11, 5, 1, 7, 3, 9),]
average_table_bioconversion_details

#write.csv(average_table_bioconversion_details,file = "Bioconversion_parameters.csv", row.names = F)