Manuscript-figures.Rmd

---
title: "FireSim2019_ManuscriptFigures"
author: "Dana Johnson"
date: "9/9/2021"
output: html_document
---
```{r setup, include=FALSE}

library(ggplot2)
library(dplyr)
library(tidyr)
library(phyloseq)
library(cowplot) 
library(vegan)
library(tidyverse)

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(out.width = '100%')

```


Defaults = theme_bw(), black = O, "organic"; grey40 - M, "mineral"; shapes = circle, triangle

# Metadata
```{r, include=FALSE}
df.meta <- read.csv('../../data/all-sample-metadata.csv')

df.meta <- subset(df.meta, core.id != '19UW-WB-11-02' &
                         Full.id != '19UW-WB-06-08-A-SI-duplicate' &
                         Full.id != '19UW-WB-07-02-O-SI-duplicate' &
                         Full.id != '19UW-WB-11-03-O-SI-duplicate' &
                         Full.id != '19UW-WB-19-07-A-SI-duplicate' &
                         Full.id != '19UW-WB-08-10-O-SI')

# Assign degree hours by horizon:
for (i in 1:nrow(df.meta)) {
  if (df.meta$PRE.O.hor.thickness.cm[i] == 10) {
    df.meta$hrzn.degree.hrs[i] = (df.meta$Low.degree.C.hours[i] + df.meta$Mid.degree.C.hours[i])/2
  } else if (df.meta$horizon[i] == 'A') {
    df.meta$hrzn.degree.hrs[i] = df.meta$Low.degree.C.hours[i]
  } else {
    df.meta$hrzn.degree.hrs[i] = df.meta$Mid.degree.C.hours[i]
  }
}

# Create labels and palettes for plots
HorizonLabels = c('Organic','Mineral')
names(HorizonLabels) = c('O','A')

DNALabels   = c('RNA','DNA')
names(DNALabels) = c('cDNA','gDNA')

ComparisonLabels = c('dry vs. unburned','wet vs. unburned')
names(ComparisonLabels) = c("DvC",'WvC')

# Edit names of Leading Species
VegLabels_Condensed = c("Picea spp.","Pinus banksiana","Populus \n tremuloides")
names(VegLabels_Condensed) <- c("Picea","Pinus_banksiana","Populus_tremuloides")

VegLabels_Full = c("Picea glauca","Picea mariana","Pinus banksiana","Pinus banksiana \n overstory","Populus \n tremuloides")
names(VegLabels_Full) <- c("Picea_glauca","Picea_mariana","Pinus_banksiana","Pinus_banksiana overstory piceglauc understory","Populus_tremula")

ThermoLabels = c('Organic-mineral \ninterface','Core base')
names(ThermoLabels) <- c('Mid','Low')

BurnLabels = c('Dry soil burn','Wet soil burn', 'Control')
names(BurnLabels) <- c('dry','wet', 'control')

# Edit names of Site
SiteLabels <- c("Site 01","Site 02","Site 03","Site 04","Site 05","Site 06",
               "Site 07","Site 08","Site 09","Site 10","Site 11","Site 12",
               "Site 13","Site 14","Site 15","Site 16","Site 17","Site 18",
               "Site 19")
names(SiteLabels) <- c("1","2","3","4","5","6","7","8","9","10","11",
                      "12","13","14","15","16","17","18","19")


BurnSeverityLabels = c('Wet burn','Dry burn')
names(BurnSeverityLabels) <- c("Low_Sev",'High_Sev')

IncubLabels = c('24 hours post-burn','5 weeks post-burn', '6 months post-burn', '6 months post-burn: \nAutoclaved','6 months post-burn: \nNot autoclaved')
names(IncubLabels) = c('pb','SI', 'LI', 'LIwA','LIn')


CoefsLabel = c('Slow pool (M2)', 'Fast pool (M1)', 'Slow pool decay rate (k2)','Fast pool decay rate (k1)')
names(CoefsLabel) = c('M2','M1', 'k2','k1')


# Reorder burn treatment levels
df.meta$horizon =factor(df.meta$horizon, levels = c('O','A'))
df.meta$incub.trtmt = factor(df.meta$incub.trtmt, levels = c('pb','SI','LIwA','LIn'))
df.meta$burn.trtmt = factor(df.meta$burn.trtmt, levels = c('control','wet','dry'))
```


# Load phyloseq objects
```{r, message = FALSE}
# Starting phyloseq:
### Normalized + tree
ps.norm.full <- readRDS('../../data/sequence-data/LibCombined/phyloseq-objects/ps.norm.full')
ps.norm.full <- prune_taxa(taxa_sums(ps.norm.full)>0, ps.norm.full)

ps.norm.full <- prune_samples(sample_data(ps.norm.full)$Full.id != '19UW-WB-06-08-A-SI-duplicate' &
                    sample_data(ps.norm.full)$Full.id != '19UW-WB-07-02-O-SI-duplicate' &
                    sample_data(ps.norm.full)$Full.id != '19UW-WB-11-03-O-SI-duplicate'&
                    sample_data(ps.norm.full)$Full.id != '19UW-WB-19-07-A-SI-duplicate' &
                    sample_data(ps.norm.full)$Full.id != '19UW-WB-08-10-O-SI', ps.norm.full)

ps.raw.full <- readRDS('../../data/sequence-data/LibCombined/phyloseq-objects/ps.raw.full')

ps.raw.full <- prune_samples(sample_data(ps.raw.full)$Full.id != '19UW-WB-06-08-A-SI-duplicate' &
                    sample_data(ps.raw.full)$Full.id != '19UW-WB-07-02-O-SI-duplicate' &
                    sample_data(ps.raw.full)$Full.id != '19UW-WB-11-03-O-SI-duplicate'&
                    sample_data(ps.raw.full)$Full.id != '19UW-WB-19-07-A-SI-duplicate' &
                    sample_data(ps.raw.full)$Full.id != '19UW-WB-08-10-O-SI', ps.raw.full)
ps.raw.full <- prune_taxa(taxa_sums(ps.raw.full)>0, ps.raw.full)

# Reorder vegetation:
sample_data(ps.norm.full)$Veg.type = factor(sample_data(ps.norm.full)$Veg.type, levels = c('Picea','Populus_tremuloides','Pinus_banksiana'))
sample_data(ps.norm.full)$horizon = factor(sample_data(ps.norm.full)$horizon, levels = c('O','A'))
sample_data(ps.norm.full)$incub.trtmt = factor(sample_data(ps.norm.full)$incub.trtmt, levels = c('pb', 'SI','LIwA','LIn'))

#Pull out just the gDNA samples:
ps.norm.gDNA <- prune_samples(sample_data(ps.norm.full)$DNA.type == 'gDNA', ps.norm.full)
ps.raw.gDNA <- prune_samples(sample_data(ps.raw.full)$DNA.type == 'gDNA', ps.raw.full)

ps.norm.cDNA <- prune_samples(sample_data(ps.norm.full)$DNA.type == 'cDNA', ps.norm.full)
ps.raw.cDNA <- prune_samples(sample_data(ps.raw.full)$DNA.type == 'cDNA', ps.raw.full)

```



# Figure 2A-C (burn temperature, pH and C:N change, respiration coefficients):
```{r}
# Temperature, Fig. 2A ----
df.Temp <-read.csv('../../data/burn-thermocouple-temp-logs/Temp-all-burns-clean.csv')
df.Temp$Thermo_position <- factor(df.Temp$Thermo_position, levels = c('Mid','Low'))

pTemp = ggplot(df.Temp, aes(x = time_hr, 
                          y = Temp)) + 
  geom_point(aes(color = burn.trtmt), size=1) +
  theme_bw()+
  facet_wrap(~Thermo_position,scales = 'free_y', ncol=1,
             labeller = labeller(Thermo_position = ThermoLabels,
                                 burn.trtmt = BurnLabels)) +
  labs(x = "Time (h)",
       y = expression(paste("Thermocouple temperature ( ",degree ~ C, ")")),
       color = "Burn treatment") + 
  #scale_color_manual(values = c("black","grey70"),
   #                  labels = c('O horizon/mineral \n interface','9 cm'),
    #                 breaks = c('wet','dry')) + 
  scale_color_manual(values = c('orange','red'),
                     limits = c('wet','dry'),
                     labels = c('Moist soil burn','Dry soil burn')) +
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size=14),
        strip.text.y = element_text(size=14), 
        strip.text.x = element_text(size=14), 
        legend.text = element_text(size=12),
        legend.title = element_text(size=12), 
        legend.position = 'right',
       strip.background = element_rect(colour="black", fill="grey92")) +
  geom_hline(yintercept=100)
pTemp

# pH and C:N change, Fig. 2B ----
pChange_pH <- ggplot(subset(df.meta, incub.trtmt == 'SI')) + 
  geom_boxplot(aes(x=horizon, y = pH, fill = burn.trtmt), varwidth = FALSE) + 
    theme_bw()+
  labs(x= 'Soil horizon',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Soil pH',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  #scale_fill_manual(values = c("deepskyblue3","saddlebrown"),
  # scale_fill_manual(values = c("slateblue4","palegreen3"),
  # #scale_fill_manual(values = c("darkgreen","saddlebrown"),
  #                    limits = c('O','A'),
  #                    labels = c('Organic','Mineral')) +
  # scale_x_discrete(limits = c('control','wet','dry'),
  #                    labels = c('Unburned\ncontrol','Wet soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")


pChange_CN  <- ggplot(subset(df.meta, incub.trtmt == 'SI' &site !=10)) + 
  geom_boxplot(aes(x=horizon, y = C.N.ratio, fill = burn.trtmt), varwidth = FALSE) + 
    theme_bw()+
  labs(x= 'Soil horizon',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Soil C:N',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  #scale_fill_manual(values = c("deepskyblue3","saddlebrown"),
  # scale_fill_manual(values = c("slateblue4","palegreen3"),
  # #scale_fill_manual(values = c("darkgreen","saddlebrown"),
  #                    limits = c('O','A'),
  #                    labels = c('Organic','Mineral')) +
  # scale_x_discrete(limits = c('control','wet','dry'),
  #                    labels = c('Unburned\ncontrol','Wet soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=12),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=12),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

cowplot::plot_grid(pChange_pH, pChange_CN, ncol=2)

# Stats ----
df.stats <- subset(df.meta, incub.trtmt == 'SI' )

test = aov(C.N.ratio ~ burn.trtmt*horizon, subset(df.stats))
summary(test)
TukeyHSD(test)




# Respiration coefficients, Fig. 2C ----
#   Decay model coefficients at end of 5 week incubation (Days=34)
#   after adjusting for total C
df.coefs <- read.csv('../../data/2-pool-decay-model-output-SI-adjusted-for-total-C.csv') %>%
  subset(t==34) 
# df.coefs <- read.csv('../../data/2-pool-decay-model-output-LIwA-adjusted-for-total-C.csv') %>%
#   subset(t==180)

df.stack <- gather(df.coefs, 
                   key = 'coefs', 
                   value = value,M1,M2,k1,k2)

df.stack$coefs <- factor(df.stack$coefs, levels = c('M1', 'M2', 'k1', 'k2'))

pResp2 = ggplot(subset(df.stack, coefs %in% c('M1','M2')), aes(x=burn.trtmt)) + 
  geom_boxplot(aes(y=value, fill=burn.trtmt),varwidth = FALSE)+
  #geom_boxplot(aes(y=log(M2), fill=burn.trtmt))+
  theme_bw()+
  labs(x= 'Burn treatment',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Fractional pool size',
       fill = 'Burn treatment',
       color = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet burn','dry burn'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('control','wet burn','dry burn'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  #scale_x_discrete(limits = c('O','A'), 
  #                 labels = c('Organic','Mineral'))
  facet_wrap(~coefs, ncol=1,scales = 'free_y',
             labeller = labeller(coefs = CoefsLabel))+
  theme(axis.text.y = element_text(size = 12),
        legend.title = element_text(size=12),
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 12, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=25, l=10)),
        strip.text = element_text(size=12),
        legend.text = element_text(size=12),
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical", 
       legend.position = '',
       strip.background = element_rect(colour="black", fill="grey92"))


pResp1 = ggplot(subset(df.stack, coefs %in% c('k1','k2')), aes(x=burn.trtmt)) + 
  geom_boxplot(aes(y=value, fill=burn.trtmt),varwidth = FALSE)+
  #geom_boxplot(aes(y=log(M2), fill=burn.trtmt))+
  theme_bw()+
  labs(x= 'Burn treatment',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Decay rate',
       fill = 'Burn treatment',
       color = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet burn','dry burn'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('control','wet burn','dry burn'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  #scale_x_discrete(limits = c('O','A'), 
  #                 labels = c('Organic','Mineral'))
  facet_wrap(~coefs, ncol=1,scales = 'free_y',
             labeller = labeller(coefs = CoefsLabel))+
  theme(axis.text.y = element_text(size = 12),
        legend.title = element_text(size=12),
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 12, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=12),
        legend.text = element_text(size=14),
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical", 
       legend.position = '',
       strip.background = element_rect(colour="black", fill="grey92"))

cowplot::plot_grid(pResp2, pResp1, ncol=2)

# Stats ----
df.stats <- subset(df.meta, incub.trtmt != 'pb' & horizon == 'O')

summary(subset(df.stats, incub.trtmt == 'LIwA' & texture != 'organic' &
                 horizon == 'O' & burn.trtmt == 'wet')$k1)
sd(subset(df.stats, incub.trtmt == 'LIwA' & burn.trtmt == 'wet')$k1)


test = aov(k2 ~ burn.trtmt, subset(df.stats, incub.trtmt == 'LIwA' & texture != 'organic'))
summary(test)
TukeyHSD(test)



test = wilcox.test(M1~burn.trtmt, data = subset(df.stats, incub.trtmt == 'SI' & burn.trtmt != 'control'))
test

test = aov(M1 ~ texture, data = subset(df.meta, incub.trtmt %in% c('LIwA','LIn') & horizon == 'O'))

test

p = ggplot(df.stats, aes(x=burn.trtmt, y =k1)) + 
  geom_boxplot() + 
  facet_grid(~incub.trtmt) +
  theme_bw()

```


# Figure 2D: RNA and DNA concentration (ng/g soil)
```{r}
df <- read.csv('../../data/raw-data/pb-RNA-and-DNA-conc.csv')

colnames(df)[1] <- 'core.id.hor.incub'

df.sub <- df.meta %>%
  subset(incub.trtmt == 'pb') %>%
  subset(select = c(core.id.hor.incub, DNA.type, horizon, burn.trtmt))

df.DNA <- subset(df, DNA.type == 'gDNA') %>%
  mutate(Conc.DNA = Concentration_ng.uL) %>%
  subset(select = c('core.id.hor.incub', 'Conc.DNA'))
df.RNA <- subset(df, DNA.type == 'cDNA') %>%
  mutate(Conc.RNA = Concentration_ng.uL) %>%
  subset(select = c('core.id.hor.incub', 'Conc.RNA'))

df <- merge(df, df.sub)

NaLabels = c('RNA','DNA')
names(NaLabels) <- c('cDNA','gDNA')


p1 = ggplot(df) + 
  geom_boxplot(aes(x=horizon, y = nucleic.acid.conc.in.soil.ng.per.g/1000, fill = burn.trtmt),varwidth = FALSE) + 
  theme_bw()+
  facet_wrap(~DNA.type, ncol=2, labeller = labeller(DNA.type = NaLabels))+
  labs(x= 'Soil horizon',
       #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Concentration\n(ug/g soil)',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  # scale_fill_manual(values = c("darkgreen","grey70"),
  #                   limits = c('O','A'),
  #                   labels = c('Organic','Mineral')) +
  # scale_x_discrete(limits = c('control','wet','dry'),
  #                  labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
        # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")
p1

p_DNA = ggplot(subset(df, DNA.type == 'gDNA')) + 
  geom_boxplot(aes(x=horizon, y = nucleic.acid.conc.in.soil.ng.per.g/1000, fill = burn.trtmt),varwidth = FALSE) + 
  theme_bw()+
  labs(x= 'Soil horizon',
       #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'DNA concentration\n(ug/g soil)',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  # scale_fill_manual(values = c("darkgreen","grey70"),
  #                   limits = c('O','A'),
  #                   labels = c('Organic','Mineral')) +
  # scale_x_discrete(limits = c('control','wet','dry'),
  #                  labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
        # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

p_RNA = ggplot(subset(df, DNA.type == 'cDNA')) + 
  geom_boxplot(aes(x=horizon, y = nucleic.acid.conc.in.soil.ng.per.g/1000, fill = burn.trtmt),varwidth = FALSE) + 
  theme_bw()+
  labs(x= 'Soil horizon',
       #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'RNA concentration\n(ug/g soil)',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned control','Moist soil burn','Dry soil burn')) +
  scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  # scale_fill_manual(values = c("darkgreen","grey70"),
  #                   limits = c('O','A'),
  #                   labels = c('Organic','Mineral')) +
  # scale_x_discrete(limits = c('control','wet','dry'),
  #                  labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
        # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

cowplot::plot_grid(p_DNA, p_RNA, ncol=2)


test <- aov(nucleic.acid.conc.in.soil.ng.per.g ~ burn.trtmt*horizon, data = subset(df, DNA.type == 'gDNA'))
summary(test)
TukeyHSD(test)


```

# Lab and field ordination, Figure 3 -----
See "Merged_paper_figures_Revision_Github.R" script

# Trait abundance in field data, Figure 4: ----
See "Merged_paper_figures_Revision_GitHub.R" script

# Rel abundance adjusted by copy num, Figure 5 -----
See "Merged_paper_figures_Revision_GitHub.R" script

# Total trait abun by BC dissimilarity, Figure 6 -----
See "Merged_paper_figures_Revision_GitHub.R" script

# Supplementary
```{r}
# Carbon and nitrogen loss, Figure S2 -----
pChange_C <- ggplot(subset(df.meta, incub.trtmt == 'SI')) + 
  geom_boxplot(aes(x=burn.trtmt, y = percent.total.C, fill = horizon), alpha = 0.7, varwidth = FALSE) + 
    theme_bw()+
  labs(x= 'Burn treatment',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Total C (%)',
       fill = 'Soil horizon') + 
  scale_fill_manual(values = c('black','grey90'),
                     limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  scale_x_discrete(limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Wet soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

pChange_N <- ggplot(subset(df.meta, incub.trtmt == 'SI')) + 
  geom_boxplot(aes(x=burn.trtmt, y = percent.total.N, fill = horizon), alpha = 0.7, varwidth = FALSE) + 
    theme_bw()+
  labs(x= 'Burn treatment',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Total N (%)',
       fill = 'Soil horizon') + 
  scale_fill_manual(values = c('black','grey90'),
                     limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  scale_x_discrete(limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Wet soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

cowplot::plot_grid(pChange_C, pChange_N, ncol=2)

test <- aov(percent.total.N ~ burn.trtmt, data = subset(df.meta, incub.trtmt == 'SI' & horizon == 'O'))
summary(test)
TukeyHSD(test)

# C remaining as fraction of total, Figure S3, S4 -----
#   See incubations_microbial_respiration.R
# Full ordination, Fig. S6 ----
ps.full.dataset <-readRDS('../../../CombinedFireSim1519/ps.1519FireSim.merged')
tree <- ape::read.tree('../../../CombinedFireSim1519/tree/merged-exported-tree/tree.nwk')

tree <- phy_tree(tree)
tax <- tax_table(ps.full.dataset)
otu <- otu_table(ps.full.dataset)
sam <- sample_data(ps.full.dataset)

ps.full.dataset = phyloseq::phyloseq(otu,tax, sam,tree)
ps.full.dataset

# Setting some plot parameters
year.labs = c("One year post-fire","Five years post-fire")
names(year.labs) = c(1,5)
horizon.labs = c("Mineral","Organic")

# Playing with new datasets
ps = ps.full.dataset

# Add column with pairs info
sample_data(ps)$Pairs = paste(sample_data(ps)$Site_ID,sample_data(ps)$Org_or_Min,sep="_")
head(sample_data(ps)$Pairs)

# Normalize sample counts
ps.norm = transform_sample_counts(ps,function(x) x/sum(x))
sample_data(ps.norm)$Study = sample_data(ps.norm)$Years_Since_Fire
sample_data(ps.norm)$Study[is.na(sample_data(ps.norm)$Study)]="Lab"


ord.nmds = ordinate(ps.norm,method = "NMDS",ddistance="wunifrac",k=3, weighted=TRUE)

plot_ordination(ps.norm,ord.pcoa,color="burn.trtmt",shape="incub.trtmt")
plot_ordination(ps.norm,ord.nmds,color="pH",shape="Study")
# Neat, they all fall out on top of each other. That's nice.
colnames(sample_data(ps.norm))


ps.norm.gDNA <- prune_samples(sample_data(ps.norm.gDNA)$incub.trtmt %in% c('pb','SI','LIwA'), ps.norm.gDNA)
ps.norm.gDNA <- prune_taxa(taxa_sums(ps.norm.gDNA)>0,ps.norm.gDNA)

MyNMDS.full.gDNA <- ordinate(ps.norm.gDNA, method = 'NMDS', k=3, ddistance="wunifrac", weighted=TRUE)

# Ordination plot with horizons and deg.hrs
p1 = plot_ordination(ps.norm.gDNA, MyNMDS.full.gDNA, axes = c(1,2), 
                      color = "burn.trtmt",
                      shape = 'incub.trtmt') + 
  facet_grid(~horizon, labeller = labeller(horizon = HorizonLabels,
                                                   DNA.type = DNALabels,
                                                   incub.trtmt = IncubLabels))+
  #facet_grid(~horizon, labeller = labeller(horizon = HorizonLabels))+
  geom_point(size=3, alpha=0.8) +
  #geom_line(aes(group = core.id.hor)) +
  labs(color = 'Burn treatment',
       #color = expression("Degree ("*degree*C*") hours"), 
       shape = 'experiment') +
  #scale_color_gradientn(colours = c('gold2','red','red4','black')) +
  scale_color_manual(values = c('black','orange','red'),
                     labels = c('Unburned control', 'Wet soil burn','Wry soil burn'),
                     breaks = c('control','wet','dry'))+
  scale_shape_manual(values = c('circle','square','triangle'),
                     labels = c('Fire survival', 'Fast growth','Post-fire envir. affinity'),
                     breaks = c('pb','SI','LIwA')) +
  #scale_shape_discrete(labels = c('Picea' = 'Picea spp.', 'Populus_tremuloides' = 'Populus tremuloides','Pinus_banksiana' = 'Pinus banksiana')) +
  #scale_shape_manual(values = c('circle','triangle'),
   #                  labels = c('O horizon','Mineral'),
    #                 breaks = c('O','A')) + 
  theme(legend.text = element_text(face = 'italic',size = 16),
        axis.title = element_text(size=16),
        axis.text = element_text(size=16),
        legend.title = element_text(size=16)) + 
  theme_bw()

p1

# BC dissimilarity and pH change, Figure S7-----
df.distances <- read.csv('../../data/bray-curtis-dissimilarites.csv')

df.distances <- merge(df.distances, subset(df.meta, burn.trtmt != 'control'))    

df.distances$incub.trtmt <- factor(df.distances$incub.trtmt, levels = c('pb','SI','LIwA','LIn'))
df.distances$horizon <- factor(df.distances$horizon, levels = c('O','A')) 
 

# plot: Dry soil burns
p3=ggplot(subset(df.distances, incub.trtmt == 'LIwA' & burn.trtmt == 'dry'),
         aes(x=change.pH, y = BC.distance.btw.burned.and.control, color = horizon, shape = horizon)) + 
  geom_point(size=3, alpha = 0.8) +
  theme_bw() +
  facet_wrap(~horizon,ncol=2,
             labeller = labeller(horizon = HorizonLabels),
             scales = 'fixed') + 
  labs(x = expression(Delta~'pH with burning'),
       y = 'Bray-Curtis dissimilarity \ncompared to control', 
       color = 'Soil horizon',
       shape = 'Soil horizon') + 
  scale_color_manual(values = c("darkgreen","grey70"),
                     limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  scale_shape_manual(values = c('circle','triangle'),
                     limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  theme(axis.title.x = element_text(size=14, margin = margin(t=10)),
        axis.title.y = element_text(size=14, margin = margin(r=10)),
        axis.text = element_text(size=12),
        legend.title = element_text(size=14),
        legend.text = element_text(size=14),
        strip.text = element_text(size=14),
        panel.spacing.x = unit(2, "lines"),
       legend.position = 'right') +
  stat_smooth(data = subset(df.distances, incub.trtmt == 'LIwA' & burn.trtmt == 'dry' & horizon == 'O'),
              method=lm,se = FALSE,color='black') +
    ylim(0,1)
  


# Respiration rate figures, S8, S9, S10-----
#   See "incubations_microbial_respiration.R" script

# Traits as fraction of total reads in lab data, Figure S11-----
df <- read.csv('../../data/sequence-data/LibCombined/corncob-output/Manuscript/Trait-responder-OTUs.csv')

ps.responders <- prune_taxa(taxa_names(ps.norm.full) %in% df$OTU, ps.norm.full)

df.melt <- psmelt(ps.responders)

responders <- df.melt %>% group_by(Sample) %>% summarize(Abundance) 
max(responders$Abundance)
mean(responders$Abundance)

dry <- subset(df.melt, burn.trtmt == 'dry') %>% subset(select = c(Sample, horizon, burn.trtmt, Abundance, incub.trtmt, Thermo.mid.max)) %>% group_by(Sample) %>% mutate(total.Abun = sum(Abundance))

summary(dry)

# Total gDNA reads in lab experiments: 39115
# Total reads in lab experiments: 43450
# Number responders = 117
# max abundance in lab data = 60%
# mean abundance, dry burns = 0.15% +/- 14%


df.frac.of.total <- df.melt %>% subset(select = c(Sample, horizon, burn.trtmt, Veg.type, Abundance, DNA.type, incub.trtmt, Thermo.mid.max, site)) %>% group_by(Sample) %>% mutate(total.Abun = sum(Abundance))


p_FractionTotal <- ggplot(subset(df.frac.of.total, DNA.type == 'cDNA'), aes(x=site, y=total.Abun, fill=burn.trtmt)) +
  geom_bar(stat='identity',  position='dodge') + 
  facet_wrap(~horizon+incub.trtmt, labeller = labeller(horizon = HorizonLabels)) + 
  theme_bw()+
  labs(x= 'Sampling site ID',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Fraction of total reads',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")

df.frac.of.total$burn.trtmt <- factor(df.frac.of.total$burn.trtmt, levels = c('control','wet','dry'))

p_FractionTotal <- ggplot(subset(df.frac.of.total,incub.trtmt!='LIn'), aes(x=horizon, y=total.Abun, fill=burn.trtmt)) +
  geom_boxplot()+
  facet_wrap(~incub.trtmt+DNA.type, labeller = labeller(incub.trtmt = IncubLabels, DNA.type = DNALabels)) + 
  theme_bw()+
  labs(x= 'Soil horizon',
     #  y = 'Decay rate constant for fast C pool \nlog(k1)',
       y = 'Fraction of total reads',
       fill = 'Burn treatment') + 
  scale_fill_manual(values = c('grey50','orange','red'),
                     limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
    scale_x_discrete(limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  theme(axis.text.y = element_text(size = 14),
        legend.title = element_text(size=14),
        axis.text.x = element_text(size = 14),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        strip.text = element_text(size=14),
        legend.text = element_text(size=14),
        legend.position = '',
       # legend.background = element_rect(fill = "darkgray"),
        legend.box = "vertical")




```

# PERMANOVAs, Supplementary:
```{r}
ps.E1 <- prune_samples(sample_data(ps.norm.full)$incub.trtmt == 'pb' &
                         sample_data(ps.norm.full)$DNA.type == 'cDNA', ps.norm.full)
ps.E1 <- prune_taxa(taxa_sums(ps.E1)>0, ps.E1)
ps.E1

ps.E1.DNA <- prune_samples(sample_data(ps.norm.full)$incub.trtmt == 'pb' &
                         sample_data(ps.norm.full)$DNA.type == 'gDNA', ps.norm.full)
ps.E1.DNA <- prune_taxa(taxa_sums(ps.E1.DNA)>0, ps.E1.DNA)
ps.E1.DNA

ps.E2 <- prune_samples(sample_data(ps.norm.full)$incub.trtmt == 'SI' &
                         sample_data(ps.norm.full)$DNA.type == 'gDNA', ps.norm.full)
ps.E2 <- prune_taxa(taxa_sums(ps.E2)>0, ps.E2)
ps.E2

ps.E3  <- prune_samples(sample_data(ps.norm.full)$incub.trtmt == 'LIwA' &
                         sample_data(ps.norm.full)$DNA.type == 'gDNA', ps.norm.full)
ps.E3 <- prune_taxa(taxa_sums(ps.E3)>0, ps.E3)
ps.E3


MyWunifrac_E1 <- distance(ps.E1, method = 'wunifrac')
MyWunifrac_E1.DNA <- distance(ps.E1.DNA, method = 'wunifrac')
MyWunifrac_E2 <- distance(ps.E2, method = 'wunifrac')
MyWunifrac_E3 <- distance(ps.E3, method = 'wunifrac')
MyWunifrac_All <- distance(ps.norm.gDNA, method = 'wunifrac')

SamDat_E1 = data.frame(sample_data(ps.E1))
SamDat_E1.DNA = data.frame(sample_data(ps.E1.DNA))
SamDat_E2 = data.frame(sample_data(ps.E2))
SamDat_E3 = data.frame(sample_data(ps.E3))
SamDat_All = data.frame(sample_data(ps.norm.gDNA))

test_E1 <- adonis(MyWunifrac_E1 ~ Veg.type+PRE.hor.thickness.cm+pH+percent.total.C+
                    percent.total.N+texture+burn.trtmt+horizon, data = SamDat_E1)
test_E1.DNA <- adonis(MyWunifrac_E1.DNA ~ Veg.type+PRE.hor.thickness.cm+pH+percent.total.C+
                    percent.total.N+texture+burn.trtmt+horizon, data = SamDat_E1.DNA)
test_E2 <- adonis(MyWunifrac_E2 ~ Veg.type+PRE.hor.thickness.cm+pH+percent.total.C+
                    percent.total.N+texture+burn.trtmt+horizon, data = SamDat_E2)
test_E3 <- adonis(MyWunifrac_E3 ~ Veg.type+PRE.hor.thickness.cm+pH+percent.total.C+
                    percent.total.N+texture+burn.trtmt+horizon, data = SamDat_E3)
test_All <- adonis(MyWunifrac_All ~ incub.trtmt+burn.trtmt+pH+horizon+Veg.type, data = SamDat_All)

test_E1
test_E1.DNA
test_E2
test_E3
test_All


test_E1_sub <- adonis(MyWunifrac_E1.DNA ~ percent.total.N, data = SamDat_E1.DNA)
test_E1_sub

# Fire survival
adonis(MyWunifrac_E1 ~ Veg.type, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ PRE.hor.thickness.cm, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ pH, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ percent.total.C, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ percent.total.N, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ texture, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ burn.trtmt, data = SamDat_E1)
adonis(MyWunifrac_E1 ~ horizon, data = SamDat_E1)

adonis(MyWunifrac_E1.DNA ~ Veg.type, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ PRE.hor.thickness.cm, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ pH, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ percent.total.C, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ percent.total.N, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ texture, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ burn.trtmt, data = SamDat_E1.DNA)
adonis(MyWunifrac_E1.DNA ~ horizon, data = SamDat_E1.DNA)


# Fast growth
adonis(MyWunifrac_E2 ~ Veg.type, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ PRE.hor.thickness.cm, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ pH, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ percent.total.C, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ percent.total.N, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ texture, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ burn.trtmt, data = SamDat_E2)
adonis(MyWunifrac_E2 ~ horizon, data = SamDat_E2)


# Test of autoclaving
ps.E3.autoclave  <- prune_samples(sample_data(ps.norm.full)$incub.trtmt %in% c('LIwA','LIn'), ps.norm.full)
ps.E3.autoclave <- prune_taxa(taxa_sums(ps.E3.autoclave)>0, ps.E3.autoclave)
ps.E3.autoclave
SamDat <- data.frame(sample_data(ps.E3.autoclave))
head(SamDat,5)

MyWunifrac_E3.autoclave <- distance(ps.E3.autoclave, method = 'wunifrac')

SamDat_E3.autoclaving = data.frame(sample_data(ps.E3.autoclave))

adonis(MyWunifrac_E3.autoclave ~ Veg.type+PRE.hor.thickness.cm+pH+percent.total.C+
                    percent.total.N+texture+burn.trtmt+horizon + incub.trtmt, data = SamDat_E3.autoclaving)

adonis(MyWunifrac_E3.autoclave ~ Veg.type, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ PRE.hor.thickness.cm, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ pH, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ percent.total.C, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ percent.total.N, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ texture, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ burn.trtmt, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ horizon, data = SamDat_E3.autoclaving)
adonis(MyWunifrac_E3.autoclave ~ incub.trtmt, data = SamDat_E3.autoclaving)



# All
adonis(MyWunifrac_All ~ Veg.type, data = SamDat_All)
adonis(MyWunifrac_All ~ pH, data = SamDat_All)
adonis(MyWunifrac_All ~ incub.trtmt, data = SamDat_All)
adonis(MyWunifrac_All ~ burn.trtmt, data = SamDat_All)
adonis(MyWunifrac_All ~ horizon, data = SamDat_All)
```








Bray-Curtis dissimilarity? 
```{r, message = FALSE}
# This is a possibly convoluted piece of code to first create a Bray-Curtis
#     dissimilarity matrix and then pull out the dissimilarities between dry/wet
#     burns and corresponding control cores. 

# First: Calculate Bray-curtis dissimilarity matrix ----
ps.cDNA.raw <- prune_samples(sample_data(ps.raw.full)$DNA.type == 'cDNA', ps.raw.full)
ps.cDNA.raw <- prune_taxa(taxa_sums(ps.cDNA.raw) > 0, ps.cDNA.raw)

dist.bray.full <- as.matrix(distance(ps.raw.full, method = 'bray', type = 'samples'))

# Look at some figures
hist(dist.bray.full)
heatmap(dist.bray.full)

# Convert distance matrix into dataframe
df <- data.frame(dist.bray.full)

# Bring in meta data and subset for specific analysis
df.meta.sub <- subset(df.meta, select = c(Full.id, burn.trtmt,site, DNA.type,
                                          incub.trtmt,core, horizon,
                                          pH, change.pH))

# Create a dataframe with the colnames and Full.ID
x <- data.frame(str_split(colnames(df)[1:460], 'X', simplify = TRUE)[,2])

# Rename columes:
colnames(x)[1] <- 'ID.bray'

# Clean up x. 
x <- x %>%
  mutate(Full.id = ID.bray)

x$Full.id <- str_replace_all(x$Full.id, "\\.", "-")

x$X <- 'X'

x <- x %>%
  unite(Bray.ID, c('X',"ID.bray"), sep = '', remove = TRUE)

# Next, pull out the control cores from the meta data and assign them to "Control.id"
df.meta.control <- df.meta %>%
  subset(burn.trtmt == 'control') %>%
  mutate(Control.id = Full.id) %>%
  subset(select = c(Control.id, site, horizon, incub.trtmt, DNA.type))  

# Combine list of control cores with full df.meta
df.meta.list <- merge(subset(df.meta.sub, select = c(Full.id, site, horizon, 
                                                     DNA.type,burn.trtmt, incub.trtmt, 
                                                     pH, change.pH)), 
                      df.meta.control, by = c('site', 'horizon', 'incub.trtmt', 'DNA.type')) 
  
# ANd then combine the Full.id list with the control core meta data I just created.
#     df.list now contains Full.id with the corresponding control sample for each 
#     dry/wet burn and a colume ("Bray.ID") that has the sample ID formatted to 
#     match the distance matrix names.
df.list <- x %>%
  merge(df.meta.list, by = 'Full.id') %>%
  subset(burn.trtmt != 'control')

df.list$burn.trtmt <- as.character(df.list$burn.trtmt)
df.list$horizon <- as.character(df.list$horizon)
df.list$incub.trtmt <- as.character(df.list$incub.trtmt)


# Create and empty data frame to store distances
df.distances <- data.frame("Full.id" = df.list$Full.id, 
                           'control.core' = 'core',
                           #'burn.trtmt' = 'burn', 'horizon' = 'H',
                           #'DNA.type' = 'DNA.type',
                           'BC.distance.btw.burned.and.control' = 'distance'
                           #'incub.trtmt' = 'incub', 'pH' = 0
                           ) %>%
  unique()


# Now run a convoluted loop:
for (i in 1:nrow(df.list)) { # For each sample (row in df.list)
  for (j in 1:nrow(df)) { # For each row in the dissimilarity matrix
    for (k in 1:ncol(df)){ # And for each colume in the dissimilarity matrix
      if(colnames(df)[k] == df.list$Bray.ID[i] &  # If the distance matrix colume ID matches the sample ID in the list
         rownames(df)[j] == df.list$Control.id[i]) { # And the rowname ID matches the corresponding control core, then...
        df.distances$BC.distance.btw.burned.and.control[i] = df[j,k] # Pull out the distance and store it in new dataframe
        #df.distances$burn.trtmt[i] = df.list$burn.trtmt[i]
        df.distances$control.core[i] = df.list$Control.id[i]
        #df.distances$horizon[i] = df.list$horizon[i]
        #df.distances$incub.trtmt[i] = df.list$incub.trtmt[i]
        #df.distances$DNA.type[i] = df.list$DNA.type[i]
        #df.distances$pH[i] = df.list$pH[i]
        #df.distances$change.pH[i] = df.list$change.pH[i]
      }
    } 
  } 
}

df.distances$BC.distance.btw.burned.and.control <- as.numeric(df.distances$BC.distance.btw.burned.and.control)

# Save backup dataframe.
df.final <- df.distances


# -----
df.distances <- read.csv('../data/bray-curtis-dissimilarites.csv')

df.distances <- merge(df.distances, subset(df.meta, burn.trtmt != 'control'))                            
dim(df.distances)

# Plots ----
p=ggplot(subset(df.distances, DNA.type == 'cDNA')) + 
  geom_boxplot(aes(x=burn.trtmt, y = BC.distance.btw.burned.and.control, fill = horizon),
               alpha = 0.7) +
  theme_bw() +
  facet_wrap(~incub.trtmt, ncol=4,
             labeller = labeller(incub.trtmt = IncubLabels)) + 
  labs(x='Burn treatment',
       y = 'Bray-Curtis dissimilarity \ncompared to control', 
       fill = 'Soil horizon') + 
  scale_fill_manual(values = c('black','grey40'),
                     limits = c('O','A'),
                     labels = c('Organic','Mineral')) +
  scale_x_discrete(limits = c('wet','dry'),
                   labels = c('Moist soil \nburn','Dry soil \nburn')) +
  theme(axis.title = element_text(size=18),
        axis.text = element_text(size=18),
        legend.title = element_text(size=18),
        legend.text = element_text(size=18),
        strip.text = element_text(size=18)) + 
  ylim(0,1)


test <- aov(BC.distance.btw.burned.and.control~burn.trtmt*horizon, data = subset(df.distances, DNA.type == 'cDNA'))
summary(test)
TukeyHSD(test)
 
```

Richness estimates?

If the community composition data of a single sample includes few singletons, we assume that sequencing captured most of the taxa present in the sample. If the data includes many singletons, we assume that sequencing missed many rare taxa. We can use breakaway to estiamte how many taxa were missed and use this to estimate richness. 
```{r}

# Load data and calculate richness ----
df.pb.gDNA.OTU <- read.csv('../data/sequence-data/LibCombined/phyloseq-objects/df.pb.gDNA.norm.csv')
df.pb.cDNA.OTU <- read.csv('../data/sequence-data/LibCombined/phyloseq-objects/df.pb.cDNA.norm.csv')
df.SI.OTU <- read.csv('../data/sequence-data/LibCombined/phyloseq-objects/df.SI.norm.csv')
df.LI.norm <- read.csv('../data/sequence-data/LibCombined/phyloseq-objects/df.LI.norm.csv')

head(df.pb.cDNA.OTU)

df.counts.pb.cDNA <- df.pb.cDNA.OTU %>%
  group_by(Full.id) %>%
  mutate(counts = length(unique(OTU))) %>%
  subset(select = c(Full.id, burn.trtmt, horizon, DNA.type, incub.trtmt, counts)) %>%
  unique()

df.counts.pb.gDNA <- df.pb.gDNA.OTU %>%
  group_by(Full.id) %>%
  mutate(counts = length(unique(OTU))) %>%
  subset(select = c(Full.id, burn.trtmt, horizon, DNA.type, incub.trtmt, counts)) %>%
  unique()

df.counts.SI <- df.SI.OTU %>%
  group_by(Full.id) %>%
  mutate(counts = length(unique(OTU))) %>%
  subset(select = c(Full.id, burn.trtmt, horizon, DNA.type, incub.trtmt, counts)) %>%
  unique()

df.counts.LI <- df.LI.norm %>%
  group_by(Full.id) %>%
  mutate(counts = length(unique(OTU))) %>%
  subset(select = c(Full.id, burn.trtmt, horizon, DNA.type, incub.trtmt, counts)) %>%
  unique()

df.counts <- rbind(df.counts.pb.cDNA, df.counts.pb.gDNA,
                   df.counts.SI, df.counts.LI)

df.counts$burn.trtmt <-  factor(df.counts$burn.trtmt, levels = c('control','wet','dry'))
df.counts$horizon <-  factor(df.counts$horizon, levels = c('O', 'A'))

df.counts.pb <- subset(df.counts, incub.trtmt == 'pb')



  

# Plots ----
p = ggplot(df.counts.pb, aes(x=burn.trtmt, y=counts)) +
  geom_boxplot(aes(fill = burn.trtmt), alpha = 0.7) + 
  geom_point(color = 'black', alpha = 0.4) + 
  theme_bw() + 
  facet_grid(horizon~incub.trtmt+DNA.type, labeller = labeller(incub.trtmt = IncubLabels,
                                                               horizon = HorizonLabels, 
                                                               DNA.type = DNALabels)) + 
  labs(y = 'unique non-zero OTU count')

test <- aov(counts ~ burn.trtmt, data = subset(df.counts.pb.gDNA, horizon =='O'))
summary(test)
TukeyHSD(test)



# Plot richness using "plot_richness" function in phyloseq: ----
ps <- prune_samples(sample_data(ps.raw.full)$incub.trtmt == 'pb',
                    ps.raw.full)

sample_data(ps)$burn.trtmt <- factor(sample_data(ps)$burn.trtmt, levels = c('control','wet','dry'))
sample_data(ps)$horizon <- factor(sample_data(ps)$horizon, levels = c('O','A'))

p = plot_richness(ps, x = "burn.trtmt", measures = 'Shannon') +
  theme_bw() + 
  geom_boxplot(aes(fill = horizon), alpha = 0.6) + 
  facet_grid(horizon~DNA.type+incub.trtmt, labeller = labeller(incub.trtmt = IncubLabels,
                                                   DNA.type = DNALabels,
                                                      horizon = HorizonLabels))
```

RNA:DNA ratios?
```{r}
df.norm.ratios <- read.csv('../data/sequence-data/LibCombined/phyloseq-objects/RNA_DNA_norm_ratios.csv')
# RNA:DNA ratios were calculated by (RNA rel. abundance)/(DNA rel. abundance)
#     In instances of DNA rel. abundance = 0, I set DNA abundance = 0.0000038, which
#     is equal to the lowest measured rel. abundance of DNA across the dataset.
#     In instances of RNA re. abundance = 0, I set ratio = 0.

df.norm.ratios$horizon <- factor(df.norm.ratios$horizon, levels = c('O','A'))
df.norm.ratios$burn.trtmt <- factor(df.norm.ratios$burn.trtmt, levels = c('control','wet','dry'))

# Plot 
p = ggplot(data = subset(df.norm.ratios, DNA.abundance > 0.001), aes(x=burn.trtmt, y = RNA.DNA.ratio, fill = horizon)) + 
  geom_boxplot()+ 
  #stat_smooth(method=lm,se = FALSE,color='black') +
  #facet_wrap(~horizon, ncol=1,scale = 'free')+
  theme_bw() + 
  labs(title = 'RNA:DNA \n*only taxa present in paired burned/unburned trtmts') 



# What fraction of total reads for each sample fall below 0.001 relative abundance threshold? -----
df.DNA.tally <- df.norm.ratios %>%
  subset(RNA.abundance + DNA.abundance > 0) %>%
  #subset(DNA.abundance > 0) %>%
  group_by(core.id.hor.incub, burn.trtmt, horizon) %>%
  mutate(counts = length(OTU)) %>%
  group_by(core.id.hor.incub, burn.trtmt, horizon,counts) %>%
  tally(DNA.abundance > 0.001) %>%
  #mutate(DNA.fraction.above.cutoff = n/counts)
  mutate(non.zero.DNA.fraction = n/counts)
head(df.DNA.tally)

df.RNA.tally <- df.norm.ratios %>%
  subset(RNA.abundance + DNA.abundance > 0) %>%
  #subset(RNA.abundance > 0) %>%
  group_by(core.id.hor.incub, burn.trtmt, horizon) %>%
  mutate(counts = length(OTU)) %>%
  group_by(core.id.hor.incub, burn.trtmt, horizon,counts) %>%
  tally(RNA.abundance > 0.001)%>%
  mutate(non.zero.RNA.fraction = n/counts)
  #mutate(RNA.fraction.above.cutoff = n/counts)
head(df.RNA.tally)


df.DNA.tally$burn.trtmt <- factor(df.DNA.tally$burn.trtmt, levels = c('control','wet','dry'))
df.RNA.tally$burn.trtmt <- factor(df.RNA.tally$burn.trtmt, levels = c('control','wet','dry'))

p = ggplot(df.DNA.tally, aes(x=burn.trtmt, y =DNA.fraction.above.cutoff, fill = horizon)) + 
  geom_boxplot() + theme_bw()

test <- aov(DNA.fraction.above.cutoff ~ burn.trtmt*horizon, data = df.DNA.tally)
summary(test)
TukeyHSD(test)

test <- aov(non.zero.RNA.fraction ~ burn.trtmt*horizon, data = df.RNA.tally)
summary(test)
TukeyHSD(test)



# Re-calculate ratio including only taxa present in both burned and unburned cores: ----
df.burned <- subset(df.norm.ratios, burn.trtmt == 'dry') %>%
  subset(DNA.abundance > 0) %>%
  mutate(burned.DNA.abundance = DNA.abundance) %>%
  subset(select = c(site, horizon, OTU, burned.DNA.abundance))

df.control <- subset(df.norm.ratios, burn.trtmt == 'control') %>%
  subset(DNA.abundance > 0) %>%
  mutate(unburned.DNA.abundance = DNA.abundance) %>%
  subset(select = c(site, horizon, OTU, unburned.DNA.abundance))


df.shared.OTUs <- merge(df.burned, df.control)

df.common.taxa <- subset(df.norm.ratios, burn.trtmt != 'wet') %>%
  merge(df.shared.OTUs) %>%
  mutate(ratio = RNA.abundance/DNA.abundance) %>%
  group_by(core.id.hor.incub) %>%
  mutate(mean.ratio = mean(ratio))


df.counts <- df.common.taxa %>%
  group_by(core.id.hor.incub, burn.trtmt, horizon) %>%
  tally(DNA.abundance > 0)

dim(df.common.taxa)

# Plot 
p = ggplot(data = df.common.taxa, aes(x=burn.trtmt, y = ratio, fill = horizon)) + 
  geom_boxplot()+ 
  #stat_smooth(method=lm,se = FALSE,color='black') +
  #facet_wrap(~horizon, ncol=1,scale = 'free')+
  theme_bw() + 
  labs(title = 'RNA:DNA \n*only taxa present in paired burned/unburned trtmts') 

test <- aov(ratio ~ burn.trtmt*horizon, data = df.common.taxa)
summary(test)
TukeyHSD(test)




# Recalculating ratios using raw reads: -----
ps.raw.gDNA <- prune_samples(sample_data(ps.raw.pb)$DNA.type == 'gDNA', ps.raw.pb)
ps.raw.gDNA
ps.raw.gDNA <- prune_taxa(taxa_sums(ps.raw.gDNA)>0, ps.raw.gDNA)
ps.raw.gDNA

ps.raw.cDNA <- prune_taxa(taxa_sums(ps.raw.cDNA)>0, ps.raw.cDNA)

df.raw.pb.cDNA <- psmelt(ps.raw.cDNA)
df.raw.pb.gDNA <- psmelt(ps.raw.gDNA)

df.RNA <- df.raw.pb.cDNA %>%
  mutate(RNA.abundance = Abundance) %>%
  subset(select = c(core.id.hor.incub, OTU, RNA.abundance))

df.DNA <- df.raw.pb.gDNA %>%
  mutate(DNA.abundance = Abundance) %>%
  subset(select = c(core.id.hor.incub, OTU, DNA.abundance, horizon, site, pH, hrzn.degree.hrs,burn.trtmt))

rownames(df.RNA) <- NULL
rownames(df.DNA) <- NULL

df.merge <- merge(df.RNA, df.DNA, by = c("core.id.hor.incub", "OTU"), all = TRUE) 

df.pb <- df.merge %>%
  group_by(site, OTU) %>%
  mutate(min = min(DNA.abundance)) %>%
  subset(min > 0) 


df.pb[is.na(df.pb$RNA.abundance),3] <- 0

df.pb <- df.pb %>%
  mutate(ratio = RNA.abundance/DNA.abundance)

# Plot
df.pb$burn.trtmt <- factor(df.pb$burn.trtmt, levels = c('control','wet','dry'))

p = ggplot(data = subset(df.pb,horizon =='O'), aes(x=burn.trtmt, y = (ratio), color = burn.trtmt)) + 
  geom_boxplot()+ 
  #stat_smooth(method=lm,se = FALSE,color='black') +
  facet_wrap(~horizon, ncol=1,scale = 'free')+
  theme_bw() + 
  labs(title = 'RNA:DNA \nonly taxa present across all trtmts')

test <- aov((ratio) ~ burn.trtmt, data = subset(df.pb, horizon == 'O'))
summary(test)
TukeyHSD(test)




# Alternative calculations for RNA:DNA ratios: ----

### Set ratio = 100 if DNA rel. abundance = 0 (Bowsher et al. 2019)
for (i in 1:nrow(df.norm.ratios)) {
  if(df.norm.ratios$DNA.abundance[i] == 0) {
    df.norm.ratios$RNA.DNA.ratio[i] = 100
  } else {
    df.norm.ratios$RNA.DNA.ratio[i] = df.norm.ratios$RNA.DNA.ratio[i]
  }
}

### If DNA abundance = 0, change it to 1.  (Bowsher et al. 2019)
#     But this would only work for raw count data.
for (i in 1:nrow(df.norm.ratios)) {
  if(df.norm.ratios$DNA.abundance[i] == 0) {
    df.norm.ratios$DNA.abundance[i] = 0.000038
  } else {
    df.norm.ratios$DNA.abundance[i] = df.norm.ratios$DNA.abundance[i]
  }
}
df.norm.ratios$RNA.DNA.ratio = df.norm.ratios$RNA.abundance/df.norm.ratios$DNA.abundance

### Add 1 to all DNA and RNA abundances (Bowsher et al. 2019)
df.norm.ratios$DNA.abundance = df.norm.ratios$DNA.abundance + 0.000038

#df.norm.ratios$RNA.abundance = df.norm.ratios$RNA.abundance + 0.000038

df.norm.ratios <- subset(df.norm.ratios, DNA.abundance > 0.0001 || RNA.abundance >0.0001 )

df.norm.ratios$RNA.DNA.ratio = df.norm.ratios$RNA.abundance/df.norm.ratios$DNA.abundance
```



# Weighted mean gene copy number
Which bacteria are fast-growing following burn treatments of varying severity and under a range of moisture conditions?
```{r}
df.WMCN.pb <- read.csv('../../data/WtMeanCopyNum_pb_gDNA.csv')
df.WMCN.SI <- read.csv('../../data/WtMeanCopyNum_SI.csv')
df.WMCN.LIwA <- read.csv('../../data/WtMeanCopyNum_LIwA.csv')

df.WMCN.SI <- subset(df.WMCN.SI, core.id != '19UW-WB-11-02' & 
                    Full.id != '19UW-WB-06-08-A-SI-duplicate' &
                    Full.id != '19UW-WB-07-02-O-SI-duplicate' &
                    Full.id != '19UW-WB-11-03-O-SI-duplicate'&
                    Full.id != '19UW-WB-19-07-A-SI-duplicate' &
                    Full.id != '19UW-WB-08-10-O-SI' )


summary(subset(df.WMCN.SI, burn.trtmt == 'control' &horizon == 'O')$WtMeanCopyNum)
sd(subset(df.WMCN.SI, burn.trtmt == 'control' &horizon == 'O')$WtMeanCopyNum)
summary(subset(df.WMCN.SI, burn.trtmt == 'control' &horizon == 'A')$WtMeanCopyNum)
sd(subset(df.WMCN.SI, burn.trtmt == 'control' &horizon == 'A')$WtMeanCopyNum)

test = aov(WtMeanCopyNum~ burn.trtmt*horizon, data = df.WMCN.SI)
summary(test)
TukeyHSD(test)

df.WMCN <- rbind(df.WMCN.SI, df.WMCN.LIwA, df.WMCN.pb)

df.WMCN <- df.WMCN %>%
  unite(core.id.hor, c(core.id, horizon), sep = '-', remove=FALSE)

df.WMCN <- merge(df.WMCN, df.meta)

# CHeck for normality:
p.QQ = qqplot(df.WMCN$hrzn.degree.hrs, df.WMCN$WtMeanCopyNum)


# Plot ----
df.WMCN$burn.trtmt <- factor(df.WMCN$burn.trtmt, levels = c('control','wet','dry'))
df.WMCN$horizon <- factor(df.WMCN$horizon, levels = c('O','A'))
df.WMCN$incub.trtmt <- factor(df.WMCN$incub.trtmt, levels = c('pb','SI','LIwA'))



  
  
p1 = ggplot(data = subset(df.WMCN, incub.trtmt != 'LIwA')) + 
  geom_boxplot(aes(x=burn.trtmt, y=WtMeanCopyNum, fill=horizon),alpha=0.7)+
  facet_wrap(~incub.trtmt, labeller = labeller(incub.trtmt = IncubLabels))+
  theme_bw() + 
  scale_fill_manual(values = c("darkgreen","grey70"),
                     labels = c('Organic','Mineral'),
                     limits = c('O','A'))+
  scale_x_discrete(limits = c('control','wet','dry'),
                     labels = c('Unburned\ncontrol','Moist soil\nburn','Dry soil\nburn')) +
  labs(x='Burn treatment',
       y='Weighted mean predicted \n16S rRNA gene copy number',
       fill = 'Soil horizon') +
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        strip.text = element_text(size=14),
        legend.position = 'right',
        axis.title.y = element_text(size = 14, margin = margin(r=10, l=10)),
        axis.title.x = element_text(size = 14, margin = margin(t=10, b=10))) +
  ylim(1,NA)


# Stats
df.WMCN.O <- subset(df.WMCN, horizon == 'O')
df.WMCN.A <- subset(df.WMCN, horizon == 'A')

fit.df.WMCN.O <- lm(WtMeanCopyNum~hrzn.degree.hrs, data=df.WMCN.O)
fit.df.WMCN.A <- lm(WtMeanCopyNum~hrzn.degree.hrs, data=df.WMCN.A)

summary(fit.df.WMCN.O)
# p = 2*10^-16, R^2 adj = 0.6081
summary(fit.df.WMCN.A)
# p = 2*10^-16, R^2 adj = 0.7883

# Weighted mean copy number STATS:
test <- aov(WtMeanCopyNum ~ burn.trtmt*horizon, data = subset(df.WMCN, incub.trtmt == 'SI'))
summary(test)
TukeyHSD(test)

test.sum <- df.WMCN %>%
  group_by(burn.trtmt, horizon) %>%
  summarize(mean.WMCN = mean(WtMeanCopyNum))





# Calculate WMCN change from pb to SI: ----
df.WMCN.initial <- subset(df.WMCN, incub.trtmt == 'pb') %>%
  mutate(WMCN.pb = WtMeanCopyNum) %>%
  subset(select = c(core.id, horizon, burn.trtmt, WMCN.pb))

df.WMCN.SI <- df.WMCN.SI %>%
  merge(df.WMCN.initial, all = TRUE) %>%
  mutate(change.WMCN = WtMeanCopyNum - WMCN.pb)


df.WMCN.SI$burn.trtmt <- factor(df.WMCN.SI$burn.trtmt, levels = c('control','wet','dry'))

p = ggplot(df.WMCN.SI, aes(x=burn.trtmt, y = change.WMCN, fill = horizon)) + 
  geom_boxplot() + 
  geom_hline(yintercept = 0)+
  theme_bw() + 
  labs(x = 'Comparison between samples 24 h and 5 wk post-burn',
       y = 'Change in predicted wt. mean gene copy number')

test <- aov(change.WMCN ~ burn.trtmt*horizon, data = df.WMCN.SI)
summary(test)
TukeyHSD(test)



# Compare copy number of positive/negative responders: ----
df.full.OTU <- read.csv('../../data/sequence-data/LibCombined/corncob-output/Manuscript/Trait-responder-OTUs.csv')
RDP <- read.csv('../../data/All-OTUs-CopyNum.csv')



df.fast <- subset(df.full.OTU, experiment == 'fast growth' &
                        Estimate > (mean(Estimate)-1*sd(Estimate))) %>%
  mutate(fast.est = Estimate)%>%
  dplyr::group_by(OTU) %>%
  dplyr::mutate(max.Estimate = max(Estimate)) %>%
  subset(Estimate == max.Estimate) %>%
  subset(select = c(OTU, Estimate, SE, t, p, p_fdr,Comparison, experiment,
                    horizon, Domain, Phylum, Class, Order, Family, Genus, 
                    Species, DNA.type, fast.est)) %>%
  unique()


ps.norm.SI <- prune_samples(sample_data(ps.norm.full)$incub.trtmt == 'SI', ps.norm.full)
ps.norm.SI <- prune_taxa(taxa_sums(ps.norm.SI)>0, ps.norm.SI)

df.not.responders <- data.frame(taxa_names(ps.norm.SI))


df.not.responders <- subset(df.not.responders, !(taxa_names.ps.norm.SI. %in% df.responders$OTU))

df.not.responders$Responders = 'no'
df.responders <- subset(df.fast, select = c(OTU)) %>%
  mutate(Responders = 'yes')
colnames(df.not.responders)[1] = 'OTU'

df.responders <- rbind(df.responders, df.not.responders) %>%
  merge(RDP)


p = ggplot(df.responders, aes(x=Responders, y = CopyNum)) + 
  geom_boxplot() + 
  theme_bw() 




df.SI.responders <- df.full.OTU %>%
  subset(experiment == 'fast growth') %>%
  merge(RDP) 

dim(df.SI.responders)

df.SI.positive = subset(df.SI.responders, Estimate > 0 ) %>%
  group_by(OTU) %>%
  mutate(Estimate = max(Estimate)) %>%
  mutate(Response = 'positive') %>%
  subset(select = c(OTU, Estimate, Comparison, horizon, Domain, Phylum, Class, Order, 
                    Family, Genus, Species, GenusRDP, CopyNum, Response)) %>%
  unique()

df.SI.negative = subset(df.SI.responders, Estimate < 0 ) %>%
  group_by(OTU) %>%
  mutate(Estimate = min(Estimate)) %>%
  mutate(Response = 'negative') %>%
  subset(select = c(OTU, Estimate, Comparison, horizon, Domain, Phylum, Class, Order, 
                    Family, Genus, Species, GenusRDP, CopyNum, Response)) %>%
  unique()

df.SI.responders <- rbind(df.SI.positive, df.SI.negative)


p = ggplot(df.SI.responders, aes(x=Estimate, y = Class, color = Phylum)) + 
  geom_point() + 
  theme_bw() + 
  facet_wrap(~Comparison) +
  geom_vline(xintercept = 0)


df.SI.dry.responders <- subset(df.SI.responders, Comparison == 'pb.dry.v.SI.dry') %>%
  subset(!(OTU %in% subset(df.SI.responders, Comparison %in% c('pb.wet.v.SI.wet', 'pb.control.v.SI.control')))) %>%
  subset(select = c(OTU, CopyNum, Response, horizon,Estimate)) %>%
  unique()
         
       
p = ggplot(subset(df.SI.dry.responders, is.na(CopyNum) == FALSE), aes(x=Response, y = CopyNum, fill = horizon)) + 
  geom_boxplot() + 
  theme_bw() + 
  facet_wrap(~horizon)

test = aov(CopyNum ~ Response+horizon, data = df.SI.dry.responders)
summary(test)
TukeyHSD(test)

p = ggplot(subset(df.SI.dry.responders, is.na(CopyNum) == FALSE), aes(x=Estimate, y = CopyNum)) + 
  geom_point()+
  facet_grid(Response~horizon) + 
  stat_smooth(method=lm,se = FALSE,color='black') +
  theme_bw() 
```