.Rhistory

phylum  <- strsplit(taxonom, split=";")[[1]][2]
gramneg <- levels(read.table(gramnegpath)[[1]])
grampos <- levels(read.table(grampospath)[[1]])
if (phylum %in% gramneg) { #FIXTHIS
result <- "--verbose -u gramneg -g M9[cit] --mediadb ../my_media.tsv"
}
else if (phylum %in% grampos) {
result <- "--verbose -u grampos -g M9[cit] --mediadb ../my_media.tsv"
}
else {
result <- ""
}
return(result)
}
carve = function(line, taxonom, outputpath, db_protein_folder) {
# Dado un string del tipo "<ID de hoja> <ID de genoma anotado>", crea un
# modelo metabólico de dicha hoja a partir del archivo correspondiente al
# genoma dado.
leaf = strsplit(line, split=" ")[[1]][1]
file = strsplit(line, split=" ")[[1]][2]
outf = paste0(outputpath,leaf,".xml")
if (file.exists(outf)) {
print(paste0(outf," model already exists. Moving to the next one..."))
} else {
system(paste0("carve ",db_protein_folder,file,"_protein.faa ",
gram(taxonom)," -o ",outf)) # nombres de archivo = nombre de hoja
print(paste0(outf," model created"))
}
}
smetana = function(pair, nodos, modelfilepath="models/", output, coupling=TRUE, output_coupling=NULL, generated_pairs_filename="generated_
pairs.txt") {
# Dado un string del tipo <hoja1> <hoja2>, lanza smetana para las hojas dadas,
# si existen sus archivos. Si no existe alguno de los dos archivos o ninguno,
# devuelve sus nombres. Si coupling==TRUE, se computan también los análisis sin
# la opción de smetana "--no-coupling".
# File creation
if (coupling==TRUE & is.null(output_coupling)) {
stop("A output name for the coupling results must be specified.")
}
filepath1 = paste0(modelfilepath,nodos[1],"/")
filepath2 = paste0(modelfilepath,nodos[2],"/")
m1 = pair[[1]]
m2 = pair[[2]]
if (file_test("-f",paste0(filepath1,m1,".xml")) & file_test("-f",paste0(filepath2,m2,".xml"))) {
for (i in c("global","detailed")){
output_filename=paste0(output,i,"/",m1,"_",m2,"_",medium)
if (!file.exists(paste0(output_filename,"_",i,".tsv"))) { #solo crea el archivo si no existía ya
system(paste0("smetana --",i," ",filepath1,m1,".xml"," ",filepath2,m2,".xml"," --flavor bigg -m ",medium,
" --mediadb ",mediadb," --molweight --no-coupling -o ",output_filename))
write(paste(m1,m2),file=paste0(output, generated_pairs_filename),append=TRUE) # lista de parejas que se han analizado con Smetana
} else {
print(paste0(output_filename,"_",i,".tsv already exists. Moving to the next pair..."))
}
}
if (coupling==TRUE) {  # se hace una ejecución más, pero solo para detailed.
output_filename_c=paste0(output_coupling,i,"/",m1,"_",m2,"_",medium)
if (!file.exists(paste0(output_filename,"_",i,".tsv"))) { #solo crea el archivo si no existía ya
system(paste0("smetana --detailed"," ",filepath1,m1,".xml"," ",filepath2,m2,".xml"," --flavor bigg -m ",medium,
" --mediadb ",mediadb," --molweight -o ",output_filename_c))
} else {
print(paste0(output_filename_c,"_detailed.tsv already exists. Moving to the next pair..."))
}
}
} else {
print("estamos aqui!!!!!!!!!!!!!!!")
return(pair) # se devuelve
}
}
find_alignment_hits = function(filepath, node_16S, nucmer_path, db_16S, showcoords, cores) {
# Runs nucmer on each of the leaves of a given node (input is 16S sequences).
# Creates 4 different files:
#    - nucmer_unfiltered: table including all the best hits for every leaf before the
#                         97% identity and 90% coverage filter.
#    - nucmer_filtered: table including the best hits with over 97% identity and 90%
#                       coverage. Some leaves may be filtered out at this point.
#    - queries_v_hits: list of every leaf that passed the filter vs its best database
#                      hit according to Nucmer. It's a two-column table
#    - genomes: list of every database hit that passed the filter. One-column table.
#
# Returns a character vector with the leaf names and hits (the same list of the
# queries_v_hits file).
# Creo un archivo temporal en formato fasta que contenga las secuencias de ese nodo
write(
simplify2array(mclapply(colnames(node_16S), FUN=function(hoja){
paste(paste(">",hoja,sep=""),as.character(node_16S[hoja][,]),sep="\n")},mc.cores=cores)),
file="sec_temp.fasta")
# Alineo cada una de sus hojas con la base de datos
nucmer(nucmer_path, db_16S, "sec_temp.fasta")
system(paste("mv out.delta",filepath))
# Descarto hits con identidad por debajo de 97% con show-coords
nucmer_res <- system(paste(showcoords," -c -l -I 97 ",filepath,
"out.delta",sep=""), intern = TRUE)
# -c  Include percent coverage information
# -l  Include the sequence length information
# -I float    Set minimum percent identity to display
# Creo archivo de mejores resultados sin filtrar por coverage
header = nucmer_res[4:5]
nucmer_res = gsub("|", "", nucmer_res[-(0:5)], fixed=T) # elimino separador
nucmer_res = gsub("(?<=[\\s])\\s*|^\\s+|\\s+$", "", nucmer_res, perl=TRUE) # fusiono espacios
write(header, paste(filepath,"nucmer_unfiltered",sep=""))
write(nucmer_res, paste(filepath,"nucmer_unfiltered",sep=""),append = TRUE) # guardo
# Creo archivo de mejores resultados filtrados por %ID y coverage y lo guardo
write(header, paste(filepath,"nucmer_filtered",sep=""))
system(paste("awk '{ if ($10>=90) { print } }' ", # filtro por cobertura
filepath,"nucmer_unfiltered > ",
filepath,"nucmer_temp1",sep=""))
system(paste("sort -r -k13,13 -k7,7 ", # ordeno por query y por %ID
filepath,"nucmer_temp1 > ",
filepath,"nucmer_temp2",sep=""))
system(paste("awk -F ' ' -v q='' '{if ($13!=q) {q=$13; print}}' ", # escojo los mejores resultados
filepath,"nucmer_temp2 >>",
filepath,"nucmer_filtered",sep="")) #guardo
# Lista filtrada de hits a modelar
queries_v_hits = system(paste("awk -F ' ' -v q='' '{if ($13!=q) {q=$13; print $13,$12}}' ",
filepath,"nucmer_temp2",sep=""), intern = TRUE)
# guardo parejas como archivo (a modo de índice informativo)
write(queries_v_hits, file=paste0(filepath, "queries_v_hits"))
# guardo solo los genomas en un archivo (útil para anotar posteriormente por ejemplo)
system(paste("awk -F ' ' -v q='' '{if ($13!=q) {q=$13; print $12}}' ",
filepath,"nucmer_temp2 > ",filepath,"genomes",sep=""), intern = TRUE)
# Borrado de archivos temporales
system(paste("rm sec_temp.fasta ", filepath,"nucmer_temp*",sep=""))
return(queries_v_hits)
}
emapper = function(input_fa, db_protein_folder, outputname, outputdir, emapper_path, cores) {
system(paste0(emapper_path," -m diamond --cpu ",cores," --no_annot --no_file_comments -i ",
db_protein_folder, input_fa,"_protein.faa"," -o ",outputname," --output_dir ",
outputdir," --temp_dir /dev/shm --dmnd_db $PWD/",dmnd_db," --override"))
returned <- system(paste0(emapper_path," --annotate_hits_table ",outputdir,"/",outputname,
".emapper.seed_orthologs -o ",outputname," --output_dir ",outputdir," --override"))
if (returned != 0) {
stop("eggNOG-mapper returned non-zero status. If your Diamond version is different from the eggNOG-mapper one
(e.g. the CarveMe version is in the PATH) please create a new database with create_compatible_database.R
and select it with --dmnd_db when next running annotate.R")
}
}
annotate = function(genomes, outputdir, db_protein_folder, emapper_path, cores) {
if (!file.exists(outputdir)){
system(paste0("mkdir ",outputdir)) }
N = length(genomes) # número total de genomas
dump <- simplify2array(mclapply(genomes,FUN=function(genome) {
print(paste0("Annotating genome ",genome," (total: ",N,")"))
emapper(input_fa=genome,
db_protein_folder = db_protein_folder,
outputname=genome,
outputdir=outputdir,
emapper_path=emapper_path,
cores=cores)
},mc.cores=cores))
}
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
abuntable=exp
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
View(checked_tipl)
check_h
sum(!checked_tipl_h[[1==checked_tipl[[1]])
sum(!checked_tipl_h[[1]]==checked_tipl[[1]])
sum(!checked_tipl_h[[2]]==checked_tipl[[2]])
sum(!(checked_tipl_h[[1]]==checked_tipl[[1]]))
sum((checked_tipl_h[[1]]==checked_tipl[[1]]))
length(unique(checked_tipl$Var1))
length(unique(checked_tipl$Var2))
length(unique(checked_tipl_h$Var2))
length(unique(checked_tipl_h$Var1))
unique(checked_tipl_h$Var2) %in%
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
unique(checked_tipl_h$Var2) %in% unique(checked_tipl$Var2)
unique(checked_tipl_h$Var2) %in% unique(checked_tipl$Var)
unique(checked_tipl_h$Var2) %in% unique(checked_tipl$Var2)
unique(checked_tipl$Var2) %in% unique(checked_tipl_h$Var2)
unique(checked_tipl$Var2)[unique(checked_tipl$Var2) %in% unique(checked_tipl_h$Var2)]
unique(checked_tipl$Var2)[!unique(checked_tipl$Var2) %in% unique(checked_tipl_h$Var2)]
extra <- unique(checked_tipl$Var2)[!unique(checked_tipl$Var2) %in% unique(checked_tipl_h$Var2)]
nodos[2]
nodos <- mclapply(node_names, function(nodo) {ape::extract.clade(tn, nodo);}, mc.cores=cores)
nodos
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
function(nodos, abuntable, cutoff = 0, cores=4) {
# Dada una pareja de nodos ("nodos", formato ape::phylo) y una tabla de
# abundancias, determina y devuelve las combinaciones de OTUs inter-nodo que
# están presentes en una misma muestra experimental.
# Cutoff: minimum abundance necessary for considering an OTU to be "present"
# in a given sample
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable <- abuntable[1:(ncol(abuntable)-1)] # exclude taxonomy column
# Seleccionamos solo las hojas que estaban en el experimento
lista_nodos <- mclapply(1:length(nodos), function(i) {
nodo <- nodos[[i]]
hojas_nodo <- hojas_nodos[[i]]
nodo$tip.label[nodo$tip.label %in% hojas_nodo]}, mc.cores=cores)
# Obtengo todas las combinaciones de hojas para cada nodo
pairs = expand.grid(lista_nodos[[1]], lista_nodos[[2]], KEEP.OUT.ATTRS = F)
# Selecciono las parejas presentes en una misma muestra
pairs[,3] = vector(mode="logical", length=length(pairs[,2])) # indicador de si coinciden o no
for (muestra in colnames(abuntable)){
# anoto qué otus hay en cada muestra
otus=rownames(subset(abuntable[muestra], abuntable[muestra]>cutoff))
# indico para cada pareja si coinciden en esa muestra o en otra ya comprobada
pairs[,3] = pairs[,3] | pairs[,1] %in% otus & pairs[,2] %in% otus
}
filtered_pairs  = subset(pairs, pairs[,3]==TRUE)[,0:2]
return(filtered_pairs)
}
function(nodos, abuntable, hojas_nodos, cutoff = 0, cores=4) {
# Dada una pareja de nodos ("nodos", formato ape::phylo) y una tabla de
# abundancias, determina y devuelve las combinaciones de OTUs inter-nodo que
# están presentes en una misma muestra experimental.
# Cutoff: minimum abundance necessary for considering an OTU to be "present"
# in a given sample
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable <- abuntable[1:(ncol(abuntable)-1)] # exclude taxonomy column
# Seleccionamos solo las hojas que estaban en el experimento
lista_nodos <- mclapply(1:length(nodos), function(i) {
nodo <- nodos[[i]]
hojas_nodo <- hojas_nodos[[i]]
nodo$tip.label[nodo$tip.label %in% hojas_nodo]}, mc.cores=cores)
# Obtengo todas las combinaciones de hojas para cada nodo
pairs = expand.grid(lista_nodos[[1]], lista_nodos[[2]], KEEP.OUT.ATTRS = F)
# Selecciono las parejas presentes en una misma muestra
pairs[,3] = vector(mode="logical", length=length(pairs[,2])) # indicador de si coinciden o no
for (muestra in colnames(abuntable)){
# anoto qué otus hay en cada muestra
otus=rownames(subset(abuntable[muestra], abuntable[muestra]>cutoff))
# indico para cada pareja si coinciden en esa muestra o en otra ya comprobada
pairs[,3] = pairs[,3] | pairs[,1] %in% otus & pairs[,2] %in% otus
}
filtered_pairs  = subset(pairs, pairs[,3]==TRUE)[,0:2]
return(filtered_pairs)
}
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
checked_tipl_h<-function(nodos, abuntable, hojas_nodos, cutoff = 0, cores=4) {
# Dada una pareja de nodos ("nodos", formato ape::phylo) y una tabla de
# abundancias, determina y devuelve las combinaciones de OTUs inter-nodo que
# están presentes en una misma muestra experimental.
# Cutoff: minimum abundance necessary for considering an OTU to be "present"
# in a given sample
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable <- abuntable[1:(ncol(abuntable)-1)] # exclude taxonomy column
# Seleccionamos solo las hojas que estaban en el experimento
lista_nodos <- mclapply(1:length(nodos), function(i) {
nodo <- nodos[[i]]
hojas_nodo <- hojas_nodos[[i]]
nodo$tip.label[nodo$tip.label %in% hojas_nodo]}, mc.cores=cores)
# Obtengo todas las combinaciones de hojas para cada nodo
pairs = expand.grid(lista_nodos[[1]], lista_nodos[[2]], KEEP.OUT.ATTRS = F)
# Selecciono las parejas presentes en una misma muestra
pairs[,3] = vector(mode="logical", length=length(pairs[,2])) # indicador de si coinciden o no
for (muestra in colnames(abuntable)){
# anoto qué otus hay en cada muestra
otus=rownames(subset(abuntable[muestra], abuntable[muestra]>cutoff))
# indico para cada pareja si coinciden en esa muestra o en otra ya comprobada
pairs[,3] = pairs[,3] | pairs[,1] %in% otus & pairs[,2] %in% otus
}
filtered_pairs  = subset(pairs, pairs[,3]==TRUE)[,0:2]
return(filtered_pairs)
}
check_h <- checked_tipl_h()
check_h <- checked_tipl_h
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
}
View(checked_tipl_h)
View(checked_tipl)
all_leaves
check_h <- function(nodos, abuntable, hojas_nodos, cutoff = 0, cores=4) {
# Dada una pareja de nodos ("nodos", formato ape::phylo) y una tabla de
# abundancias, determina y devuelve las combinaciones de OTUs inter-nodo que
# están presentes en una misma muestra experimental.
# Cutoff: minimum abundance necessary for considering an OTU to be "present"
# in a given sample
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable <- abuntable[1:(ncol(abuntable)-1)] # exclude taxonomy column
# Seleccionamos solo las hojas que estaban en el experimento
lista_nodos <- mclapply(1:length(nodos), function(i) {
nodo <- nodos[[i]]
hojas_nodo <- hojas_nodos[[i]]
nodo$tip.label[nodo$tip.label %in% hojas_nodo]}, mc.cores=cores)
# Obtengo todas las combinaciones de hojas para cada nodo
pairs = expand.grid(lista_nodos[[1]], lista_nodos[[2]], KEEP.OUT.ATTRS = F)
# Selecciono las parejas presentes en una misma muestra
pairs[,3] = vector(mode="logical", length=length(pairs[,2])) # indicador de si coinciden o no
for (muestra in colnames(abuntable)){
# anoto qué otus hay en cada muestra
otus=rownames(subset(abuntable[muestra], abuntable[muestra]>cutoff))
# indico para cada pareja si coinciden en esa muestra o en otra ya comprobada
pairs[,3] = pairs[,3] | pairs[,1] %in% otus & pairs[,2] %in% otus
}
filtered_pairs  = subset(pairs, pairs[,3]==TRUE)[,0:2]
return(filtered_pairs)
}
check_h <- function(nodos, abuntable, hojas_nodos, cutoff = 0, cores=4) {
# Dada una pareja de nodos ("nodos", formato ape::phylo) y una tabla de
# abundancias, determina y devuelve las combinaciones de OTUs inter-nodo que
# están presentes en una misma muestra experimental.
# Cutoff: minimum abundance necessary for considering an OTU to be "present"
# in a given sample
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable <- abuntable[1:(ncol(abuntable)-1)] # exclude taxonomy column
# Seleccionamos solo las hojas que estaban en el experimento
lista_nodos <- mclapply(1:length(nodos), function(i) {
nodo <- nodos[[i]]
hojas_nodo <- hojas_nodos[[i]]
print(length(hojas_nodo))
nodo$tip.label[nodo$tip.label %in% hojas_nodo]}, mc.cores=cores)
# Obtengo todas las combinaciones de hojas para cada nodo
pairs = expand.grid(lista_nodos[[1]], lista_nodos[[2]], KEEP.OUT.ATTRS = F)
# Selecciono las parejas presentes en una misma muestra
pairs[,3] = vector(mode="logical", length=length(pairs[,2])) # indicador de si coinciden o no
for (muestra in colnames(abuntable)){
# anoto qué otus hay en cada muestra
otus=rownames(subset(abuntable[muestra], abuntable[muestra]>cutoff))
# indico para cada pareja si coinciden en esa muestra o en otra ya comprobada
pairs[,3] = pairs[,3] | pairs[,1] %in% otus & pairs[,2] %in% otus
}
filtered_pairs  = subset(pairs, pairs[,3]==TRUE)[,0:2]
return(filtered_pairs)
}
checked_tipl_h <- check_h(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
hojas_nodos = all_leaves,
cores = cores)
node_names
nodo
nodo$tip.label==136777
nodo$tip.label[nodo$tip.label=="136777"]
nodos[1]$tip.label[nodos[1]$tip.label=="136777"]
nodos[[1]]$tip.label[nodos[[1]]$tip.label=="136777"]
nodos[[2]]$tip.label[nodos[[2]]$tip.label=="136777"]
nodos[[1]]
nodos[[2]]
exp
checked_tipl <- check(nodos = nodos,
cutoff = 0.005,
abuntable = abuntable,
cores = cores)
View(checked_tipl)
View(checked_tipl)
node_names
nodos[[2]]$node.label
nodos[[2]]$node.label==nodos[[1]]
for (n in nodos[[2]]$node.label) {if(n==nodos[[1]]) {print(n)}}
for (n in nodos[[2]]$node.label) {if(n==node_names[[1]]) {print(n)}}
for (n in nodos[[1]]$node.label) {if(n==node_names[[2]]) {print(n)}}
nodos_16S
checked_tipl
mclapply(checked_tipl, function(n) {print((unique(n))); fasta[n]}, mc.cores=cores)
nodos_16S    <- mclapply(checked_tipl, function(n) {fasta[unique(n)]}, mc.cores=cores)
nodos_16S
mclapply(nodos, function(nodo) {fasta[nodo$tip.label]}, mc.cores=cores)
# Para cada nodo creo una carpeta de resultados
for (i in c(1:length(nodos))) {
filepath   <- paste("models/",nodos[i],"/",sep="")
system(paste("mkdir", filepath)) # la carpeta tendrá el mismo nombre que el nodo
# Alineo cada hoja de cada nodo con Nucmer y obtengo un genoma adecuado para cada una,
# seleccionando solo las hojas cuyos hits pasan un filtro de calidad
nucmer_res_final <- find_alignment_hits(filepath, nodos_16S[[i]], nucmer_path, db_16S, showcoords, cores)
# Modelado con CarveMe de todas las hojas de cada nodo que pasan el filtro de Nucmer
print(paste0("Creating models for ",nodos[i],"..."))
dump <- simplify2array(mclapply(nucmer_res_final,
FUN = function(line) {carve(line, taxonom[i], filepath, db_protein_folder)},mc.cores=cores))
print(paste0("Finished modelling for ",nodos[i],"."))
}
getwd()
system("cd repos/TFM")
getwd()
setwd("repos/TFM/")
# Para cada nodo creo una carpeta de resultados
for (i in c(1:length(nodos))) {
filepath   <- paste("models/",nodos[i],"/",sep="")
system(paste("mkdir", filepath)) # la carpeta tendrá el mismo nombre que el nodo
# Alineo cada hoja de cada nodo con Nucmer y obtengo un genoma adecuado para cada una,
# seleccionando solo las hojas cuyos hits pasan un filtro de calidad
nucmer_res_final <- find_alignment_hits(filepath, nodos_16S[[i]], nucmer_path, db_16S, showcoords, cores)
# Modelado con CarveMe de todas las hojas de cada nodo que pasan el filtro de Nucmer
print(paste0("Creating models for ",nodos[i],"..."))
dump <- simplify2array(mclapply(nucmer_res_final,
FUN = function(line) {carve(line, taxonom[i], filepath, db_protein_folder)},mc.cores=cores))
print(paste0("Finished modelling for ",nodos[i],"."))
}
filepath
filepath   <- paste("models/",node_names[i],"/",sep="")
# Para cada nodo creo una carpeta de resultados
for (i in c(1:length(nodos))) {
filepath   <- paste("models/",node_names[i],"/",sep="")
system(paste("mkdir", filepath)) # la carpeta tendrá el mismo nombre que el nodo
# Alineo cada hoja de cada nodo con Nucmer y obtengo un genoma adecuado para cada una,
# seleccionando solo las hojas cuyos hits pasan un filtro de calidad
nucmer_res_final <- find_alignment_hits(filepath, nodos_16S[[i]], nucmer_path, db_16S, showcoords, cores)
# Modelado con CarveMe de todas las hojas de cada nodo que pasan el filtro de Nucmer
print(paste0("Creating models for ",node_names[i],"..."))
dump <- simplify2array(mclapply(nucmer_res_final,
FUN = function(line) {carve(line, taxonom[i], filepath, db_protein_folder)},mc.cores=cores))
print(paste0("Finished modelling for ",node_names[i],"."))
}
dump
nucmer_res_final
checked_tipl
nodos_16S
View(nodos_16S)
nodos
nodos
length(nodos)
tab <- na.omit(read.csv("temp",sep = "\t"))
tab
View(tab)
getwd()
exp2="temp"
abuntable="temp"
tab <- na.omit(read.csv("temp",sep = "\t"))
node_names <- tab[[1]]
taxonom    <- tab[[10]] # for Gram-type checking
taxonom
node_names
node_names <- tab[[1]]
taxonom    <- tab[[10]] # for Gram-type checking
nodos <- mclapply(node_names, function(nodo) {ape::extract.clade(tn, nodo)}, mc.cores=cores)
# ------------------------------------------------
# 3 --> alinear cada hoja, filtrar y hacer modelos
# ------------------------------------------------
# Asigno secuencias 16S a cada hoja
if (checking == TRUE) {
# De cada hoja que pase el checking, tomo la secuencia de 16S del fasta original.
# Son todas las hojas de cada nodo que estén presentes en la tabla de abundancia
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
nodos_16S    <- mclapply(checked_tipl, function(n) {fasta[unique(n)]}, mc.cores=cores)
} else {
# De cada hoja (sin checking), cojo la secuencia de 16S del fasta original
nodos_16S    <- mclapply(nodos, function(nodo) {fasta[nodo$tip.label]}, mc.cores=cores)
}
nodos
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
nodos = nodos
if (length(nodos) > 2) {
warning(paste("The node name list is longer than 2. ", length(nodos)," Only the first two nodes will be checked for coexistence"))
nodos <- nodos[1:2]
}
nodos
# Leo abundancias
abuntable <- read.csv(abuntable, sep="\t",skip = 1,row.names=1)
abuntable
abuntable=exp
checked_tipl <- check(nodos = nodos,
cutoff = 0,
abuntable = abuntable,
cores = cores)
checked_tipl
View(checked_tipl)
nodos_16S    <- mclapply(checked_tipl, function(n) {fasta[unique(n)]}, mc.cores=cores)
node_names
node_names[[2]] %in% nodos[[1]]$node.label
node_names[[1]] %in% nodos[[2]]$node.label
node_names[[2]] %in% nodos[[1]]$node.label
tab <- na.omit(read.csv("~/AAA/OTU_data/glucosa/results.txt",sep = "\t"))
node_names <- tab[[1]]
taxonom    <- tab[[10]] # for Gram-type checking
nodos <- mclapply(node_names, function(nodo) {ape::extract.clade(tn, nodo)}, mc.cores=cores)
node_names[[2]] %in% nodos[[1]]$node.label
node_names[[2]] %in% nodos[[3]]$node.label
node_names[[1]] %in% nodos[[2]]$node.label
tab <- na.omit(read.csv("~/AAA/2022-09-13_tomate/Final_results_221010/tomate_bc_result_100/Tree/results.txt",sep = "\t"))
tab
node_names <- tab[[1]]
taxonom    <- tab[[10]] # for Gram-type checking
node_names[[2]] %in% nodos[[3]]$node.label
node_names
nodos <- mclapply(node_names[1:2], function(nodo) {ape::extract.clade(tn, nodo)}, mc.cores=cores)
node_names[1:2]
node_names[[2]] %in% nodos[[3]]$node.label
node_names[[2]] %in% nodos[[1]]$node.label
node_names[[1]] %in% nodos[[2]]$node.label