#!/usr/local/bin/Rscript

###Instalar librerias

install.packages('pacman')
BiocManager::install("phyloseq")
remotes::install_github("MadsAlbertsen/ampvis2")
remotes::install_github("adw96/breakaway")
BiocManager::install("limma")
BiocManager::install("DECIPHER")
install.packages("phangorn")
BiocManager::install("tidyverse")

####Cargar librerias

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
if (!requireNamespace("dada2", quietly = TRUE))
  BiocManager::install("dada2", version = "3.15")
library("dada2")
library(pacman)
packages <- c('phyloseq', 'Biostrings', 'ggplot2', 'tidyverse', 'vegan', 'decontam', 'ape', 'DESeq2', 'microbiome', 'metagenomeSeq', 'remotes')
pacman::p_load(char = packages)
library(tidyverse)
library(limma)
library(phyloseq); packageVersion("phyloseq")
library(phangorn)
library(dada2)
library(DECIPHER)
library(ggplot2); packageVersion("ggplot2")
theme_set(theme_bw())

####DADA2 PIPELINE 2.0 - DATASET PROYECTO FUSARIUM

#Establecer ruta a archivos limpios (usar buscar y reemplazar para cambiar la ruta)
path <- "/Users/esteb/My Documents/Universidad/atbr/tesis"
list.files(path)

#Los archivos forward y reverse tienen el formato: _1.fastq.gz y _2.fastq.gz
fnFs <- sort(list.files(path, pattern="_R1_001.fastq.gz", full.names = TRUE)) #para correr despues de trimmomatic
fnRs <- sort(list.files(path, pattern="_R2_001.fastq.gz", full.names = TRUE))

#Extraer los nombres de las muestras, asumiendo que tienen el formato: NOMBRE_XXX.fastq.gz
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

#Inspect read quality profiles
plotQualityProfile(fnFs[1:9]) 
plotQualityProfile(fnRs[1:9])

## Trimming
# Place filtered files in filtered/ subdirectory
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE)                  #Si va a cortar las muestras, usar truncLen=(240,230)

head(out)

## Trimming2
# Place filtered files in filtered/ subdirectory
fnFs2 <- sort(list.files(path, pattern="_R1_001.fastq.gz", full.names = TRUE)) #para correr despues de trimmomatic
fnRs2 <- sort(list.files(path, pattern="_R2_001.fastq.gz", full.names = TRUE))

#Extraer los nombres de las muestras, asumiendo que tienen el formato: NOMBRE_XXX.fastq.gz
sample.names2 <- sapply(strsplit(basename(fnFs2), "_"), `[`, 1)
filtFs2 <- file.path(path, "filtered2", paste0(sample.names2, "_F_filt.fastq.gz"))
filtRs2 <- file.path(path, "filtered2", paste0(sample.names2, "_R_filt.fastq.gz"))
names(filtFs2) <- sample.names2
names(filtRs2) <- sample.names2

out <- filterAndTrim(fnFs2, filtFs2, fnRs2, filtRs2, 
                     maxN=0, truncLen=240,230, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE)                  #Si va a cortar las muestras, usar truncLen=(240,230)

head(out)

#Learn the Error Rates (Toma 10 min aprox)
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

#Plot error
plotErrors(errF, nominalQ=TRUE)
plotErrors(errR, nominalQ=TRUE)

#Guardar datos
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/err.rda")
#Cargar datos guardados
load(file="/Users/esteb/My Documents/Universidad/atbr/tesis/err.rda")

##Dereplicate (this step does not be included in version 1.16 but it is in version 1.8)
derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

##Name the derep-class objects by the sample names
names(derepFs) <- sample.names
names(derepRs) <- sample.names

#Sample inference
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)  #continuando el proceso luego de la dereplicacion añadir pool=TRUE
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

dadaPFs <- dada(derepFs, err=errF, multithread=TRUE, pool=TRUE)  
dadaPRs <- dada(derepRs, err=errR, multithread=TRUE, pool=TRUE)

#Inspect returned class object (Todo lo siguiente con P se refiere a pool)
dadaFs[[1]]
dadaRs[[1]]
dadaPFs[[1]]
dadaPRs[[1]]

#Merge paired reads to obtain the full denoised sequences
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)  #con dereplicacion
mergersP <- mergePairs(dadaPFs, derepFs, dadaPRs, derepRs, verbose=TRUE)

# Inspect the merger data.frame from the first sample
head(mergers[[1]])
head(mergersP[[1]])

#REVISION (Si se pierden muchas secuencias probalmente hay que revisar el corte y filtro)
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/merged.rda")
load(file="/Users/esteb/My Documents/Universidad/atbr/tesis/merged.rda")

#Construct sequence table
seqtab <- makeSequenceTable(mergers)
seqtabP <- makeSequenceTable(mergersP)

## The sequences being tabled vary in length.
dim(seqtab)
dim(seqtabP)

# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))
table(nchar(getSequences(seqtabP)))

#Remove chimeras
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(seqtab)

seqtab.nochimP <- removeBimeraDenovo(seqtabP, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochimP)
sum(seqtab.nochimP)/sum(seqtabP)

##REVISIÓN FINAL del PROCESAMIENTO
#Track reads trough pipeline 
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))

getN <- function(x) sum(getUniques(x))
trackP <- cbind(out, sapply(dadaPFs, getN), sapply(dadaPRs, getN), sapply(mergersP, getN), rowSums(seqtab.nochimP))

# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)

colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)
track

write.csv(track, file.path('/Users/esteb/My Documents/Universidad/atbr/tesis/load/track.csv')) 

colnames(trackP) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(trackP) <- sample.names
head(trackP)
trackP

saveRDS(track, "/Users/esteb/My Documents/Universidad/atbr/tesis/summary_table.rds")
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/summary_table.rda")

###ASSIGN TAXONOMY###
taxa <- assignTaxonomy(seqtab.nochimP, "/Users/esteb/My Documents/Universidad/atbr/tesis/silva_nr_v138_train_set.fa.gz", multithread=TRUE)

#ASSIGN SPECIES STEP
taxa <- addSpecies(taxa, "/Users/esteb/My Documents/Universidad/atbr/tesis/silva_species_assignment_v138.fa.gz", verbose=TRUE)  #this function Only those genus-species binomials which are consistent with the genus assigned in the provided taxonomy table are retained in the output

#Inspect taxonomic assignments
taxa.print <- taxa #Crea una copia de la tabla de taxa
taxa.print <- genus.species  #Usando assignSpecies (no suele funcionar)
rownames(taxa.print) <- NULL # Removing sequence rownames for display only

#Si se quiere retomar los nombres de taxa en taxa.print:
rownames(taxa.print) <- rownames(taxa)

head(taxa.print)

##FINAL de DADA2
saveRDS(seqtab.nochimP, "/Users/esteb/My Documents/Universidad/atbr/tesis/seqtab.nochimP.rds")
saveRDS(taxa, "/Users/esteb/My Documents/Universidad/atbr/tesis/taxa.rds")
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/dada2.rda")

# PHYLOSEQ PIPELINE

#Extraer los metadatos de la codificación de muestras

samples.out <- rownames(seqtab.nochimP)
sitio <- sapply(strsplit(samples.out, "M"), `[`, 1)
tratamiento <- substr(sitio,1,1)
muestra <- sapply(strsplit(samples.out, "M"), `[`, 2)
tejido <- substr(muestra,1,1)
replica <- substr(muestra,2,2)
samdf <- data.frame(Tratamiento=tratamiento, Tejido=tejido, Réplica=replica)
samdf$Análisis <- "Deflagración"
samdf$Análisis[samdf$Tratamiento<2] <- "Remanente"
samdf$Análisis[samdf$Tratamiento>2] <- "Predetonación"
rownames(samdf) <- samples.out

#Almacenar tabla de metadatos
write.csv(samdf, file = file.path(path, 'metadatos.csv'))
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/tablamet.rda") #Punto de recuperación antes del arbol

#Crear un archivo fasta con los ASVs

#Preparar secuencias y ancabezados
secuencias <- colnames(seqtab.nochimP)
encabezados <- paste('>ASV', 1:ncol(seqtab.nochimP), sep = '_')

#Interlinear
asv_fasta <- c(rbind(encabezados, secuencias))

#Almacenar archivo fasta en el outdir
write(asv_fasta, file = file.path(path, 'ASVs.fa'))

#(PROCESAR ARBOL) Generar un arbol filogenetico en cluster. Instalar el docker de Bioconductor e ingresar, ver archivo Docker.txt donse einstala el container con el nombre biocpersonal. Instalar y cargar librerias necesarias en R ingresando desde linux con el codigo "R"

#Preparar secuencias
seqs <- getSequences(seqtab.nochimP)
names(seqs) <- seqs

#Alinear
alignment <- AlignSeqs(DNAStringSet(seqs),
                       anchor = NA,
                       iterations = 5,
                       refinements = 5)
phang.align <- phyDat(as(alignment, 'matrix'), type = 'DNA')
dm <- dist.ml(phang.align)
treeNJ <- NJ(dm) 
fit = pml(treeNJ, data = phang.align)
fitGTR <- update(fit, k = 4, inv = 0.2)
fitGTR <- optim.pml(fitGTR,
                    model = 'GTR',
                    optInv = TRUE,
                    optGamma = TRUE,
                    rearrangement = 'stochastic',
                    control = pml.control(trace = 0)) #Este paso es el más largo y puede no finalizar (, epsilon =0.001, maxit = 1000)

#Una opción para acelerar la optimización del arbol es correrlo con todos los nodos del CPU (Seleccionar todo y correr)

library(parallel)

cl <- makeCluster(detectCores()) # Crea un cluster con el número de núcleos disponibles
clusterEvalQ(cl, library(phangorn))
clusterExport(cl, "fitGTR") # Exporta la variable fitGTR al cluster

# Ejecuta la optimización en paralelo
fitGTR <- parLapply(cl, 1:9, function(i) {
  optim.pml(fitGTR,
            model = 'GTR',
            optInv = TRUE,
            optGamma = TRUE,
            rearrangement = 'stochastic',
            control = pml.control(trace = 1, epsilon =0.001, maxit = 1000))
})

stopCluster(cl) # Detiene el cluster

#Exportar los datos a un objeto de phyloseq
ps <- phyloseq(otu_table(seqtab.nochimP, taxa_are_rows=FALSE), sample_data(samdf), tax_table(taxa), phy_tree(fitGTR$tree)) #probar con taxa.print en vez de taxa (utilizar phy_tree aunque no se tenga el arbol para evitar problemas después)

#Guardar imagen posterior al arbol y phyloseq
save.image(file="/home/esteban.mena/esteban/atbr/secuencias/arbol.rda")
saveRDS(ps, file = file.path(path, 'ps.rds'))
saveRDS(fitGTR, '/Users/esteb/Documents/Universidad/atbr/tesis/load/fitGTR2.rds')

#Tambien se puede guardar individualmente
saveRDS(alignment, '/home/esteban.mena/esteban/atbr/secuencias/alignment.rds')
saveRDS(phang.align, '/Users/esteb/My Documents/Universidad/atbr/tesis/phang.align.rds')
saveRDS(treeNJ, '/Users/esteb/My Documents/Universidad/atbr/tesis/treeNJ.rds')
saveRDS(fit, '/Users/esteb/My Documents/Universidad/atbr/tesis/fit.rds')
saveRDS(dm, '/Users/esteb/My Documents/Universidad/atbr/tesis/dm.rds')

#Si se abre una nueva sesión, cargar archivos
fitGTR <- readRDS(file.path(path, 'fitGTR.rds'))
ps <- readRDS(file.path(path, 'ps.rds'))

#Cambiar el nombre de los ASVs en el objeto ps
dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- taxa_names(ps)

# Juntar el objeto de secuencias al objeto de phyloseq:
ps<- merge_phyloseq(ps, dna)

#Renombrar ASVs
taxa_names(ps) <- paste("ASV", 1:ntaxa(ps), sep = "_")

#Observar el cambio en los primeros 3 ASVs
taxa_names(ps)[1:3]

#En caso de tener información de la concentración de ADN de cada muestra postPRC se puede remover posibles contaminantes con el paquete "decontam"

##Eliminar taxones eucariotas

# Crear subsets para cloroplast, mitocondria, and eukaryotes:
chlr <- subset_taxa(ps, Order == "Chloroplast")
mit <- subset_taxa(ps, Family == "Mitochondria")
euk <- subset_taxa(ps, Kingdom == "Eukaryota")

# Generar una variablew con todos los taxones a remover:
elim_tax <- c(taxa_names(chlr), taxa_names(mit), taxa_names(euk))

# Mantener taxones deseados:
taxa <- taxa_names(ps)
des_tax <- taxa[!(taxa %in% elim_tax)]

# Objeto ps reducido conteniendo solo subset deseado:
ps_red <- prune_taxa(des_tax, ps)

#Determinar proporción que se mantuvo en p_red
pre <- sum(sample_sums(ps))
post <- sum(sample_sums(ps_red))
post / pre

##Filtrar muestras poco informativas

#Conteo por muestra
sums <- sample_sums(ps)
sums[order(sums)]

#Conteo más bajo
min(sums)

#Muestra con conteo más bajo
names(sums[order(sums)][1])

#Remover muestras con menos de 1000 conteos (Si se va a almacenar en otra variable, considerar utilizar esa para los respectivos procesos)
ps <- subset_samples(ps_red, sample_sums(ps_red) > 1000)

#Observar datos en el objeto
sample_data(ps)[1:9]
head(t(otu_table(ps)[1:9]))
head(tax_table(ps))

#Almacenar ol objeto de phyloseq reducido
saveRDS(ps, file = file.path(path, 'phyloseq_red.rds'))

##Normalizacion

# Normalizacion por proporcion:
library(microbiome)
ps_prop <- transform(ps, "compositional")

# Observar tabla de OTUs:
otu_table(ps_prop)[1:9, 1:9]

# Observar el conteo por muestra, todos deben sumar 1:
head(sample_sums(ps_prop))

#Graficar abundancia SIN NORMALIZAR a nivel de phylum
r_phylum_sn <- tax_glom(ps, taxrank = "Phylum", NArm = FALSE)
r_phylum_sn <- subset_taxa(r_phylum_sn, taxa_sums(r_phylum_sn) > 0.1)
plot_bar(r_phylum_sn, fill = "Phylum") +
  facet_wrap(~Análisis, scales = "free_x", nrow = 1)

#Graficar abundancia NORMALIZADA a nivel de phylum
r_phylum <- tax_glom(ps_prop, taxrank = "Phylum", NArm = FALSE)
r_phylum <- subset_taxa(r_phylum, taxa_sums(r_phylum) > 0.1)
plot_bar(r_phylum, fill = "Phylum") +
  facet_wrap(~Análisis, scales = "free_x", nrow = 1)

#Abundancia NORMALIZADA a nivel de genero 
r_genus <- tax_glom(ps_prop, taxrank = "Genus", NArm = FALSE)
r_genus <- subset_taxa(r_genus, taxa_sums(r_genus) > 0.1)
plot_bar(r_genus, fill = "Genus") +
  facet_wrap(~Análisis, scales = "free_x", nrow = 1)

#Abundancia SIN NORMALIZAR(?) a nivel de genero (siempre sale solo la leyenda)
r_genus_sn <- tax_glom(ps_red, taxrank = "Genus", NArm = FALSE)
r_genus_sn <- subset_taxa(r_genus_sn, taxa_sums(r_genus_sn) > 0.1)
plot_bar(r_genus_sn, fill = "Genus") +
  facet_wrap(~Análisis, scales = "free_x", nrow = 1)


#Diversidad Alfa

#Curvas de rarefaccion
library(vegan)
rarecurve(otu_table(ps), step = 500, xlab = "Sample Size", ylab = "Taxa") #No suele funcionar

#Curvas de rarefaccion por grupo
library(ampvis2)
metadata <- data.frame(sample_data(ps), check.names = TRUE)
asv_table <- data.frame(t(otu_table(ps)),
                        tax_table(ps),
                        check.names = FALSE)
library(data.table)
metadata2=setDT(metadata, keep.rownames = "SampleID")
asv_table2=setDT(asv_table, keep.rownames = "ASV")

#Cargar objeto de ampvis2
ps3 <- amp_load(otutable=asv_table2, metadata=metadata2)
rar_plot <- amp_rarecurve(ps3,
                          stepsize = 100,
                          facet_by = "Análisis",
                          color_by = "Réplica") +
  ylab("ASVs observados")
rar_plot

#Indices de diversidad, usar ps (antes de normalizar, no ps_prop) pues singletons son importantes para estos indices
plot_richness(ps,
              x = "Análisis",
              shape = "Réplica",
              color = "Réplica", 
              measures = c("Observed", "Shannon", "Simpson", "InvSimpson"))

#Punto de guardado
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/divalfa.rda")
saveRDS(ps, '/Users/esteb/My Documents/Universidad/atbr/tesis/psdivalfa.rds') #Ya está reducido
saveRDS(ps_red, '/Users/esteb/My Documents/Universidad/atbr/tesis/ps_red_divalfa.rds')
saveRDS(ps_prop, '/Users/esteb/My Documents/Universidad/atbr/tesis/ps_prop_divalfa.rds')

#Diversidad Beta
#Calcular distancia wunifrac (requiere arbol) (PCoA)
set.seed(100)
wUF <- phyloseq::distance(ps_prop, method = "wunifrac")
as.matrix(wUF)[1:9, 1:9]

#Calcular distancia bray-curtis
bray <- vegan::vegdist(otu_table(ps_prop), method = "bray")

#Obsevar parte de la matriz
as.matrix(bray)[1:9,1:9]

#Análisis de coordinadas principales
do <- ordinate(ps_prop, method = "PCoA", distance = bray)
baseplot <- plot_ordination(ps_prop, do,
                            type = "samples",
                            color = "Réplica",
                            shape = "Análisis") +
  ggtitle("PCoA, distancia de Bray Curtis") +
  theme_bw()
baseplot +
  geom_point(size = 2)

# Por Tratamiento:
baseplot <- plot_ordination(ps_prop, do,
                            type = "samples",
                            color = "Réplica") +
  ggtitle("PCoA, distancia de Bray Curtis") +
  theme_bw()
baseplot +
  geom_point(size = 2) + 
  ggtitle("PCoA, distancia de Bray Curtis, Remanente") +
  facet_wrap(~Análisis, nrow = 1)

# Por Réplica:
baseplot <- plot_ordination(ps_prop, do,
                            type = "samples",
                            color = "Análisis") +
  ggtitle("PCoA, distancia de Bray Curtis") +
  theme_bw()
baseplot +
  geom_point(size = 2) + 
  ggtitle("PCoA, distancia de Bray Curtis, Réplica") +
  facet_wrap(~Réplica, nrow = 1)

# Por Tejido:
baseplot <- plot_ordination(ps_prop, do,
                            type = "samples",
                            color = "Análisis") +
  ggtitle("PCoA, distancia de Bray Curtis") +
  theme_bw()
baseplot +
  geom_point(size = 2) + 
  ggtitle("PCoA, distancia de Bray Curtis, Tejido") +
  facet_wrap(~Tejido, nrow = 1)

#Subset metadatos y estimar UniFrac para Raíz:
metadatosR <- as(sample_data(ps_R), "data.frame")
unifrac_dist <- UniFrac(ps_R, weighted = TRUE)

#PCoA con UniFrac
doUni <- ordinate(ps_R, method = "PCoA", distance = unifrac_dist)
baseplot <- plot_ordination(ps_R, doUni,
                            type = "samples",
                            color = "Análisis",
                            shape = "Réplica") +
  ggtitle("PCoA, distancia de UniFrac") +
  theme_bw()
baseplot +
  geom_point(size = 2)

#Mapa de calor. Los datos de ps3 serán normalizados pues amp_heatmap está por defecto con la opción normalise=true
ps3subset <- amp_subset_samples(ps3, Tejido %in% "R")
amp_heatmap(ps3subset, group_by = "Análisis")
amp_heatmap(ps3, facet_by = "Tejido", group_by="Análisis", tax_aggregate="Genus") 

#Mapa de calor para Tratamientos
ps3subset_a <- amp_subset_samples(ps3, Análisis %in% "Remanente")
(ps3subset_a, group_by = "Réplica")
amp_heatmap(ps3, facet_by = "Análisis", group_by="Réplica", tax_aggregate="Genus") 
amp_heatmap(ps3, facet_by = "Análisis", group_by="Réplica", tax_aggregate="Phylum") 

#Permanova, necesita paquete vegan
#Subset de acuerdo a tejido
ps_rem <- subset_samples(ps_prop, Tratamiento == "1")
ps_Defl <- subset_samples(ps_prop, Tratamiento == "2")
ps_Predet <- subset_samples(ps_prop, Tratamiento == "3")

ps_R <- subset_samples(ps_prop, Tejido == "R")

#Permanova utilizando adonis2 de vegan
set.seed(100)  # Establecer seed parsa hacer analisis reproducible
permanova <- adonis2(unifrac_dist ~ Análisis,
                     data = metadatosR,
                     permutations = 10000)
permanova

write.csv(permanova, file.path('/Users/esteb/My Documents/Universidad/atbr/tesis/load/permanova.csv'))

save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/permanova.rda")
saveRDS(ps_R, '/Users/esteb/My Documents/Universidad/atbr/tesis/ps_R.rds')

#Abundancia diferencial con deseq2

#Crear objeto deseq2 a partir de objeto de phyloseq no normalizado
#Cargar objeto de phyloseq de ser necesario
library(phyloseq)
library(DESeq2)
dds <- phyloseq_to_deseq2(ps, ~Análisis)
dds
colData(dds)

#Modelo Univariado para contraste de pares
#Agrupar los dos factores, Tratamiento (Análisis)
dds$Group <- factor(paste0(dds$Análisis))
design(dds) <- formula(~Group)
dds <- estimateSizeFactors(dds, type = 'poscounts')
dds <- DESeq(dds)
results(dds)

#Contrastar (quizá unicamente es necesario resD)
resD <- results(dds, contrast = c("Group", "Remanente", "Deflagración"))
summary(resD)
resP <- results(dds, contrast = c("Group", "Remanente", "Predetonación"))
summary(resP)
resDP <- results(dds, contrast = c("Group", "Deflagración", "Predetonación"))
summary(resDP)
results <- results(dds)
summary(results)

# Numero de taxa con abundancia significativamente superior en "RD":
sum(resD$padj < 0.1 & resD$log2FoldChange > 0, na.rm = TRUE)
# Numero de taxa con abundancia significativamente inferior en "RD":
sum(resD$padj < 0.1 & resD$log2FoldChange < 0, na.rm = TRUE)


# Numero de taxa con abundancia significativamente superior en "RD":
sum(results$padj < 0.1 & results$log2FoldChange > 0, na.rm = TRUE)
# Numero de taxa con abundancia significativamente inferior en "RD":
sum(results$padj < 0.1 & results$log2FoldChange < 0, na.rm = TRUE)

#Ordenar resultados por p-value y observar los primeros ASVs:
resD.order <- resD[order(resD$pvalue), ]
results.order <- results[order(results$pvalue), ]
head(resD.order)[c(1:2, 6)]
head(results.order)[c(1:2, 6)]
plotCounts(dds, gene = results.order@rownames[1], intgroup = "Group")
plotCounts(dds, gene = "ASV_42", intgroup = "Group")

#Mejorar gráfico con ggplot:
library(ggplot2)
d <- plotCounts(dds, gene = resD.order@rownames[1], intgroup = "Group", returnData = TRUE)
ggplot(d, aes(x = Group, y = count)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 2) +
  scale_y_log10() +
  ggtitle(expression(paste("ASVMax, log"[10], "-scale"))) +
  theme_bw()

#Escoger solo taxa con p-value < 0.05:
sig <- subset(resD.order, padj <= 0.05)
head(sig)

#Taxa associado con ASVs con abundancia diferencial
sig_asv <- rownames(sig)
sig_taxa <- tax_table(ps)[sig_asv,]
head(sig_taxa)
sig_write <- cbind(sig, sig_taxa)
write.csv(sig_write, file.path(path, 'DATaxa_PS_PA.csv'))  

# Numero de taxa con abundancia significativamente superior en "DP":
sum(resDP$padj < 0.1 & resDP$log2FoldChange > 0, na.rm = TRUE)
# Numero de taxa con abundancia significativamente inferior en "DP":
sum(resDP$padj < 0.1 & resDP$log2FoldChange < 0, na.rm = TRUE)

#Ordenar resultados por p-value y observar los primeros ASVs:
resDP.order <- resDP[order(resDP$pvalue), ]
head(resDP.order)[c(1:2, 6)]
plotCounts(dds, gene = resDP.order@rownames[1], intgroup = "Group")
plotCounts(dds, gene = "ASV_42", intgroup = "Group")

#Mejorar gráfico con ggplot:
library(ggplot2)
dDP <- plotCounts(dds, gene = resDP.order@rownames[1], intgroup = "Group", returnData = TRUE)
ggplot(dDP, aes(x = Group, y = count)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 2) +
  scale_y_log10() +
  ggtitle(expression(paste("ASVMax, log"[10], "-scale"))) +
  theme_bw()

#Escoger solo taxa con p-value < 0.05:
sigDP <- subset(resDP.order, padj <= 0.05)
head(sigDP)

#Taxa associado con ASVs con abundancia diferencial
sigDP_asv <- rownames(sigDP)
sigDP_taxa <- tax_table(ps)[sigDP_asv,]
head(sigDP_taxa)
sigDP_write <- cbind(sigDP, sigDP_taxa)
write.csv(sigDP_write, file.path(path, 'DATaxa_DP.csv'))

#Mapas de calor
library("RColorBrewer")
library("gplots")
genes <-  rownames(resD.order[1:20,]) heat
counts(dds[genes,])
mat <-  counts(dds[genes,])

heatmap(t(mat))

#Mapa de calor tipo 2
genes1p <-  subset(resD.order, padj <= 0.01)
dim(genes1p)
select <- order(rowMeans(counts(dds,normalized=TRUE)),decreasing=TRUE)[1:30]
hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)

heatmap.2(counts(dds,normalized=TRUE)[select,], col = hmcol,
          Rowv = FALSE, Colv = FALSE, scale = "none",
          dendrogram="none", trace="none", margin=c(10,6))
dev.off()

##Final de pipeline Phyloseq
save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/heatmap.rda")
saveRDS(sig, '/Users/esteb/My Documents/Universidad/atbr/tesis/sig.rds')
saveRDS(resD, '/Users/esteb/My Documents/Universidad/atbr/tesis/resD.rds')

#seleccionar solo Pseudotallo (?)
matP=mat[, -grep("R", colnames(mat))]
matP=matP[, -grep("S", colnames(matP))]
heatmap(t(matP))
heatmap(t(matP), scale="column")
amp_heatmap(mat, group_by = "Estado")

##Predicción funcional (En terminal)
#Instalar picrust2
wget https://github.com/picrust/picrust2/archive/v2.4.1.tar.gz
tar xvzf  v2.4.1.tar.gz
cd picrust2-2.4.1/
  
  #Create and activate the environment (with requirements) and then install PICRUSt2 with pip.
  
  export PATH="/home/esteban.mena/miniconda3/bin:$PATH"
conda env create -f  picrust2-env.yaml
conda activate picrust2
pip install --editable .

#Run the tests to verify the installation (note that this wont work if you followed the bioconda install instructions):

pytest

#Crear una tabla de ASVs normalizada a partir de mi objeto de phyloseq
> otu_prop = otu_table(ps_prop)
> asv_prop = as.data.frame.matrix(otu_prop)

#Transponer tabla seqtab.nochim
pcrust = t(asv_prop)
#Crear archivo csv
write.csv(pcrust, "/Users/esteb/My Documents/Universidad/atbr/tesis/pcrust.csv")
#Crear un aschivo tsv
write.table(pcrust, file="/Users/esteb/My Documents/Universidad/atbr/tesis/pcrust.tsv", sep = "\t", row.names = TRUE, col.names = NA)

save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/befpicr.rda")

#Correr picrust2
conda activate picrust2
#export PATH=$PATH:/raid/home/francisco.flores/PROGRAMS/picrust2-2.4.1/scripts
picrust2_pipeline.py -s /home/francisco.flores/BANANO/16S_ASVs.fas -i /home/francisco.flores/BANANO/pcrust.tsv -o /home/francisco.flores/BANANO/picrust2_out_pipeline3 -p 1
metagenome_pipeline.py -f -i /home/francisco.flores/BANANO/pcrust.tsv -o /home/francisco.flores/BANANO/metagenome_predictions.tab

#Hacer un PCA en R
mydata <- read.table("C:/Users/franc/OneDrive - WULEL/IDGEN/Cursos/Metagenómica/rutas.tsv", header=TRUE, sep="\t", row.names = 1)

#Ver que columnas no tienen varinaza
names(mydata[, sapply(mydata, function(v) var(v, na.rm=TRUE)==0)])

#Remover columnas sin varianza
md2=mydata[,apply(mydata, 2, var, na.rm=TRUE) != 0]
#Transponer
tmd=t(md2)

#Ver que filas no tienen varianza
rownames(mydata)[apply(mydata, 2, var, na.rm=TRUE) == 0]

#Remover filas sin varianza
md3=mydata[apply(mydata, 2, var, na.rm=TRUE) < 0.02, ]
md3=mydata[apply(mydata, 2, var, na.rm=TRUE) != 0,]
tmd3=t(md3)
tmd.pca3 <- prcomp(tmd3, center = TRUE, scale = TRUE)

#PCA
tmd.pca <- prcomp(tmd, center = TRUE, scale. = TRUE)
tmd_pca2 <- princomp(tmd, cor = FALSE, scores = TRUE)
summary (tmd.pca)

#Plot PCA
library(devtools)
install_github("vqv/ggbiplot")
library(ggbiplot)
ggbiplot(tmd.pca)
ggbiplot(tmd.pca3, 
         obs.scale = 1, 
         var.scale = 1, 
         groups = dds$group,
         ellipse = TRUE, 
         circle = TRUE,
         varname.size = 2,
         varname.adjust = 2,
         alpha = 0.5,
         labels.size = 1,
         labels.shorten = 1) # Acortar los nombres de etiquetas a 10 caracteres


tmd.scale <- scale(tmd3, center = TRUE, scale = TRUE)
tmd.cor <- cor(tmd.scale)
# Instalar y cargar la librería ggcorrplot
install.packages("ggcorrplot")
library(ggcorrplot)

# Seleccionar las variables con una correlación mayor a 0.8
tmd.corr_high <- cor(tmd.scale)
tmd.corr_high[abs(tmd.corr_high) < 0.8] <- 0

# Graficar la matriz de correlación
ggcorrplot(tmd.corr_high, hc.order = TRUE, type = "lower", outline.color = "white",
           colors = c("#6D9EC1", "white", "#E46726"), legend.title = "Correlation", 
           lab = TRUE, lab_size = 3.5, method = "circle")


#Install factoextra
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

library(factoextra)

#Visualizar componentes principales
fviz_eig(tmd.pca)

#Ver contribucion de variables a componentes principales
fviz_contrib(tmd.pca, choice = "var", axes = 1, top = 10)
fviz_contrib(tmd.pca, choice = "var", axes = 2, top = 10)

#Ver PCA agrupado por sintoma
fviz_pca_ind(tmd.pca,
             label = "none", # hide individual labels
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE # Concentration ellipses
)

save.image(file="/Users/esteb/My Documents/Universidad/atbr/tesis/aftPCA.rda")

#Metacoder
#Utilizar archivos metadata, asv_table y taxa.table *para generar tax_table
asv_tablet <-otu_table(ps)
asv_tablet = as.data.frame.matrix(asv_tablet)
tax.table <-tax_table(ps)
tax.table = as.data.frame.matrix(tax.table) 
library(readr)
sample_data <- read_tsv(file.path(indir, 'metadatos.txt'),
                        col_types = "ccccc")
install.packages("metacoder")
library(metacoder)
mc_obj <- parse_dada2(
  asv_tablet,
  tax.table,
  class_key = "taxon_name",
  class_regex = "(.*)",
  include_match = TRUE
)
print(mc_obj$data$asv_table)
mc_obj$data$asv_props <- calc_obs_props(mc_obj, "asv_table", other_cols = TRUE)
mc_obj$data$tax_abund <- calc_taxon_abund(mc_obj, "asv_props")
mc_obj$data$type_abund <- calc_group_mean(mc_obj, "tax_abund",
                                          cols = sample_data$Sample_ID,
                                          groups = sample_data$Tejido)
#Subsample Pseudotallo
ps_P <- subset_samples(ps, Tejido == "P")
asv_tabletP <-otu_table(ps_P)
asv_tabletP = as.data.frame.matrix(asv_tablet)
tax.tableP <-tax_table(ps_P)
tax.tableP = as.data.frame.matrix(tax.tableP) 
mc_objP <- parse_dada2(
  asv_tabletP,
  tax.tableP,
  class_key = "taxon_name",
  class_regex = "(.*)",
  include_match = TRUE
)
print(mc_obj$data$asv_tableP)
mc_objP$data$asv_props <- calc_obs_props(mc_objP, "asv_table", other_cols = TRUE)
mc_objP$data$tax_abund <- calc_taxon_abund(mc_objP, "asv_props")
mc_objP$data$type_abund <- calc_group_mean(mc_objP, "tax_abund",
                                           cols = sample_data$Sample_ID,
                                           groups = sample_data$Estado)
mc_objP$data$diff_table <- compare_groups(mc_objP,
                                          dataset = "tax_abund",
                                          cols = sample_data$Sample_ID,
                                          groups = sample_data$Estado)
print(mc_objP$data$diff_table)	
set.seed(999)
heat_tree(mc_objP, 
          node_label = taxon_names,
          node_size = n_obs, # n_obs is a function that calculates, in this case, the number of ASVs per taxon
          node_color = log2_median_ratio, # A column from `mc_objP$data$diff_table`
          node_color_interval = c(-2, 2), # The range of `log2_median_ratio` to display
          node_color_range = c("cyan", "gray", "tan"), # The color palette used
          node_size_axis_label = "Conteo de ASVs",
          node_color_axis_label = "Radio Log 2 de proporciones de la media",
          layout = "davidson-harel", # The primary layout algorithm
          initial_layout = "reingold-tilford") # The layout algorithm that initializes node locations