library(phyloseq)
library(ggplot2)
library(pals)

metadata <- read.table("metadata.txt", header = T,dec=",")
metadata$Type<-as.factor(metadata$Type)
metadata$Size<-factor(metadata$Size, levels=c("U","S","M","L"))
metadata$Sample<-factor(metadata$Sample, levels=c("gDNA_U", "gDNA_S", "gDNA_M", "gDNA_L", "cDNA_U", "cDNA_S", "cDNA_M", "cDNA_L"))
rownames(metadata)<-metadata[,1]

#Los nombres de las muestras de los archivos metadata y de abundancia de ASV (seqtab.nochim) deben tener el mismo codigo, para eso:
identical(rownames(metadata), rownames(seqtab.nochim))

#Para crear un objeto phyloseq necesito mínimo archivo de abundancias de ASV y taxonomía, pero puedo sumarle metadata (variables que me interesen para graficar) y un árbol filogenético.
ps<- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE), tax_table(taxa2), 
              sample_data(metadata),phy_tree(fitGTR$tree))

#En este caso trabajo sin arbol. Lo siguiente solo se puede hacer si no sumo arbol al phyloseq.
#Para que las secuencias de las ASV queden guardadas en phyloseq primero hay que copiar los nombres de las filas del taxa2 en donde estan y unirle una columna que sean las secuencias (en el archivo taxa2 las secuencias de los ASV están como los row.names).
taxa3<-cbind(taxa2,row.names(taxa2))
colnames(taxa3)[8]<-"Secuencias"
rownames(taxa3)<-paste0("ASV_", seq(nrow(taxa3))) #nombro a las ASV como ASV_1, ASV_2,..; en lugar de las secuencias nucleotidicas.
saveRDS(taxa3, "taxa3.rds")
colnames(seqtab.nochim_b)<-rownames(taxa3) #le cambio al archivo seqtab.nochim los nombres de las columnas que eran las secuencias por ASV_1, ASV_2,...

## Crear el objeto phyloseq ##
ps<- phyloseq(otu_table(seqtab.nochim_b, taxa_are_rows=FALSE),
              sample_data(metadata),tax_table(taxa3))
saveRDS(ps, "ps.rds")

#### FILTRADO PRIMARIO ####
##FILTRO 1: Sacar Archeas, y ADN de Cloroplastos y Mitocondria
##Cloroplastos (filtrar de order) y Mitocondria (filtrar de family)
ps1=subset_taxa(ps, !Kingdom %in% "Archaea" & !Order %in% "Chloroplast" & !Family %in% "Mitochondria")

## FILTRO 2: Se puede descartar los que no tengan mas de 10 conteos(ej) dado que hay muuchas ASV con 0 en todas las muestras:
ps2 <- prune_taxa(taxa_sums(ps1) > 10, ps1) # Eliminar las que tienen menos de 10 conteos en general

#Luego veo las prevalencias de los filos (se puede hacer a otros niveles taxonómicos) para ver cual filtar
# Create table, number of features for each phyla
table(tax_table(ps2)[, "Phylum"], exclude = NULL)
# Compute prevalence of each feature, store as data.frame
prevdf = apply(X = otu_table(ps2),
               MARGIN = ifelse(taxa_are_rows(ps2), yes = 1, no = 2),
               FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf = data.frame(Prevalence = prevdf,
                    TotalAbundance = taxa_sums(ps2),
                    tax_table(ps2))
plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})

ggplot(prevdf, aes(TotalAbundance, Prevalence / nsamples(ps2),color=Phylum)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) + geom_point(size = 2, alpha = 0.7) +
  scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_wrap(~Phylum) + theme(legend.position="none")

# FILTRO 3: Sacar aquellos phyla con menor a 10 reads (segun tabla de prevalencia)
filterPhyla = c("Abditibacteriota", "Caldisericota", "Campylobacterota","Dependentiae", "DTB120", "Elusimicrobiota", "FCPU426", "Hydrogenedentes","Latescibacterota", "Methylomirabilota","NB1-j","Nitrospinota","Poribacteria","Sumerlaeota", "Sva0485","WPS-2","WS4", "Zixibacteria")
# Procedemos a filtrar
ps3 = subset_taxa(ps2, !Phylum %in% filterPhyla)

#Nota: también se pueden quitar los filos que no se lograron asignar a ninguno ("NA"), pero dependerá del objetivo del trabajo, quizás es interesante incluirlos. En tal caso, una forma es:
ps0 <- subset_taxa(ps, !is.na(Phylum) & !Phylum %in% c("", "uncharacterized"))

#Se puede ver cuántas ASV se eliminaron de cada muestra, comparando entre el phyloseq sin filtrar y el último
sort(sample_sums(ps))
sort(sample_sums(ps3))

#Problema: Silva pone a los Burkholderiales como Gamma y son Beta, los modifique:
tax<-as(tax_table(ps3),"matrix")
write.csv2(tax,"tax.csv") #lo modifico en excel a mano
tax<-as.matrix(read.csv2("tax.csv",header=T, row.names = 1,sep=",")) #subo el modificado
tax_table(ps3)<-tax #lo ingreso al phyloseq
#lo chequeo
tax_b<-as(tax_table(ps3),"matrix")

# Para unir asignaciones taxonómicas y abundancias de las ASV por muestra:
#Extraer tabla con abundancia de seq por muestra
otu<-as(otu_table(ps3), "matrix")
#Extraer tabla con seq asiganadas taxonomicamente
tax<-as(tax_table(ps3), "matrix")
#UNIR TAXA Y OTU EN UN DATA FRAME#
otu<-t(otu)
identical(row.names(tax),row.names(otu)) #[1] TRUE
asv_table<-cbind(otu,tax)
write.csv2(asv_table,"asv_table.csv")

###Curva de rarefaccion###
library(vegan)
rarec <- function (x, step = 1, sample, xlab = "Sample Size", ylab = "Species", 
                   label = TRUE, cols = c(rep('#8968CD', nrow(x) / 2), rep('#43CD80', nrow(x) / 2)), ...) {
  tot <- rowSums(x)
  S <- specnumber(x)
  nr <- nrow(x)
  out <- lapply(seq_len(nr), function(i) {
    n <- seq(1, tot[i], by = step)
    if (n[length(n)] != tot[i]) 
      n <- c(n, tot[i])
    drop(rarefy(x[i, ], n))
  })
  Nmax <- sapply(out, function(x) max(attr(x, "Subsample")))
  Smax <- sapply(out, max)
  plot(c(1, max(Nmax)), c(1, max(Smax)), xlab = xlab, ylab = ylab, 
       type = "n", ...)
  if (!missing(sample)) {
    abline(v = sample)
    rare <- sapply(out, function(z) approx(x = attr(z, "Subsample"), 
                                           y = z, xout = sample, rule = 1)$y)
    abline(h = rare, lwd = 0.5)
  }
  for (ln in seq_len(length(out))) {
    N <- attr(out[[ln]], "Subsample")
    lines(N, out[[ln]], col = cols[ln], ...)
  }
  if (label) {
    ordilabel(cbind(tot, S), labels = rownames(x), ...)
  }
  invisible(out)
}
otu<-as(otu_table(ps3), "matrix")
rarec(otu, step = 20,  cex = 0.6, ylab="Number of observed ASVs",xlab="Sequencing depth (reads)")

#guardar imagen
jpeg(filename = "rarefvegan.jpeg",width = 8, height =6, units = "in",quality = 100,bg = "white", res = 300)
rarec(otu, step = 20,  cex = 0.6, ylab="Number of observed ASVs",xlab="Sequencing depth (reads)")
dev.off()

#### Alfa Diversidad ####
library(gridExtra) 
#Antes sacar Microcystis
ps3_sM<-subset_taxa(ps3,!Family %in% "Microcystaceae")
otu<-as(otu_table(ps3_sM), "matrix")

Shannon=diversity(otu, index= "shannon")
Simpson=diversity(otu, index= "simpson")
Richness=specnumber(otu)
Richness<-as.numeric(Richness)
Equit=Shannon/log(Richness)
#Le agrego al metadata los indices de diversidad
sample_data(ps3_sM)<-cbind(sample_data(ps3_sM),Shannon,Simpson,Richness,Equit)
saveRDS(ps3_sM, "ps3_sM.rds")

sam<-as(sample_data(ps3_sM), "data.frame")
sam$Size<-factor(sam$Size, levels=c("U","S","M","L")) #Así al graficar quedan en ese orden

#Incluyo solo muestras de ADN (fila 1:4)
a<-ggplot(sam[-c(5:10),],aes(x=Size, y=Richness))+ 
  geom_col(alpha = 0.3) #geom_col(alpha = 0.3,fill="blue")
b<-ggplot(sam[-c(5:10),],aes(x=Size, y=Shannon))+ 
  geom_col(alpha = 0.3)
grid.arrange(a,b,nrow=1)
c<-ggplot(sam[-c(5:10),],aes(x=Size, y=Simpson))+ 
  geom_col(alpha = 0.3)
d<-ggplot(sam[-c(5:10),],aes(x=Size, y=Equit))+ 
  geom_col(alpha = 0.3)
grid.arrange(a,b,c,d)

###Beta diiversidad###
library(ggrepel)
GP.ord <- ordinate(ps3_sM, "PCoA", "bray")
plot_ordination(ps3_sM, GP.ord)+
  geom_point(size=4, alpha=2,aes(color=Type)) +geom_text_repel(aes(label=Sample),vjust=0.3) +theme(legend.position = "none")

#### Grafico Barras phylum ####
sample_data(ps3)$Size<-factor(sample_data(ps3)$Size, levels=c("U","S","M","L")) #Así al graficar quedan en ese orden
ps3_100 = transform_sample_counts(ps3, function(x){(x/sum(x))*100})
phylumGlommed = tax_glom(ps3_100, "Phylum")
plot_bar(phylumGlommed,x="Size",fill="Phylum")+ scale_fill_manual(values=as.vector(stepped(24)))+ylab("Relative abundance") +
  facet_grid(~Type)+
  theme(axis.text.y =element_text(size=10),
        axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5,size=10))

#Repito lo mismo sin Microcystis:
ps3_100 = transform_sample_counts(ps3_sM, function(x){(x/sum(x))*100})
phylumGlommed = tax_glom(ps3_100, "Phylum")
plot_bar(phylumGlommed,x="Size",fill="Phylum")+ scale_fill_manual(values=as.vector(stepped(24)))+ylab("Relative abundance") +
  facet_grid(~Type)+
  theme(axis.text.y =element_text(size=10),
        axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5,size=10)) 


####Correlacion de ASVs entre ADN y ADNc de los distintos tamaños####
library(corrplot)
otuM<-as(otu_table(ps3_sM), "matrix")
M <- cor(otuM,method = "spearman") 
corrplot(M,method="number",tl.col = "black",type="lower")