Tornando as matrizes comparativas

NOTA: essa compilação mistura sequencias de individuos diferentes (ou ao menos não sabe se é ou não o mesmo) para a mesma espécie, portanto, a análise desses dados pode ter resultados errados, considerando que a taxonomia de amostras depositadas em GeneBank não é padronizada. Neste contexto de dados secundários, pode até fazer sentido fazer isso, mas idealmente sequencias de diferentes marcadores devem vir do mesmo indivíduo.

O script abaixo pega teus arquivos ALINHADOS em fasta para cada marcador e gera arquivos (PARACOMPARACAO.fasta). Ele pega o nome das espécies dos nomes das sequencias (o nome da espécie deve ser o primeiro elemento antes de "_" ou mude o script para ajustar outros formatos) e cria:

• Um arquivo para cada marcador chamado marcador_PARACOMPARACAO.fasta com uma linha por especie, todos os marcadores com as mesmas espécies e na mesma ordem. Se para o marcador não houver dado para uma espécie ela é incluida como dado faltante "-".
• O arquivo meusdadosconcatenados.nex com os marcadores concatenados na mesma linha, e indicação nesse nexus das partições (genes, marcadores) amostrados.
• O arquivo tabelacomGENBANKidsDasSequenciasUsadas.csv - tabela com os identificadores originais das sequencias utilizadas para cada espécie.
• Você deve definir seus arquivos e outras opções no inicio do script.

#install.packages("ape", dependencies = T) #se nao tiver instalado esse pacote
library(ape)
#pasta que você está
getwd()
#configure se não estiver na pasta onde estao os arquivos abaixo
#adapta isso para os seus dados, um objeto para cada marcador DEVE ESTAR ALINHADO 
its = read.FASTA("ITS.fasta")
ets = read.FASTA("ETS.fasta")
phyc = read.FASTA("PhyC.fasta")
rps16 = read.FASTA("rps16.fasta")
trnlf = read.FASTA("trnL-trnF.fasta")
 
#junte seus marcadores num vetor, atribuindo a eles o nome dos objetos
meusmarcadores = list(its=its,ets=ets,phyc=phyc,rps16=rps16,trnlf=trnlf)
 
#defina o numero minimo de marcadores para incluir a especie
nmin = 2
 
#defina uma subpasta da pasta atual, onde quer salvar os arquivos gerados
subpasta = "dados"
#cria a subpasta se nao houver
dir.create(subpasta,showWarnings=F)
 
########DAQUI PARA BAIXO NAO PRECISA MUDAR NADA ###########
 
#qual marcador tem o maior numero de especies?
pegaespecie = function(x) {strsplit(x,"_")[[1]][1]}
pegaespecies <- function(x) { 
  spp = sapply(names(x),pegaespecie)
  length(unique(spp))
}
tt = sapply(meusmarcadores,pegaespecies)
#numero de especie em cada marcador
tt
 
#cria uma tabela com os identificadores de cada especie
asespecies = NULL
for(n in 1:length(meusmarcadores)) {
  asespecies = c(asespecies,sapply(names(meusmarcadores[[n]]),pegaespecie))
}  
#funcoes uteis
mytrim <- function (x) {gsub("^\\s+|\\s+$", "", x)}
esogaps<-function(x) { sum(x=="-")==length(x)}
 
#para cada especie nos dados tabula os ids dos tips dos alinhamentos
aspp = unique(sort(mytrim(asespecies)))
resultado = data.frame(ESPECIE=aspp)
rownames(resultado) = aspp
jafoi = vector(mode='list',length=length(meusmarcadores))
names(jafoi) = names(meusmarcadores)
for(n in 1:length(aspp)) {
  sp = tolower(aspp[n])
  for(o in 1:length(meusmarcadores)) {
    omarc  = meusmarcadores[[o]]
    if (n==1) {
      resultado[,names(meusmarcadores)[o]] = NA
    }
    spps = sapply(names(omarc),pegaespecie)
    spps = tolower(mytrim(spps))
    if (sum(spps==sp)>0) {
      valores = names(spps)[spps==sp]
      if (!is.null(jafoi[[names(meusmarcadores)[o]]])) {
        jf = jafoi[[names(meusmarcadores)[o]]]
        if (sum(!valores%in%jf)>0) {
          valor = valores[!valores%in%jf][1]
        } else { valor = NA}
      } else { valor = valores[1]}
      jafoi[[names(meusmarcadores)[o]]]  = c(jafoi[[names(meusmarcadores)[o]]] ,valor)
      resultado[aspp[n],names(meusmarcadores)[o]] = valor
    }
  }
}  
 
#quantos marcadores?
contana <- function(x) { sum(!is.na(x)) }
gp = apply(resultado,1,contana)
#quais tem pelo menos nmin marcadores
quais = gp>=nmin
#para essas sequencias cria uma tabela para cada conjunto de dados
tbb= resultado[quais,]
 
#agora cria um conjunto de dados para cada marcador
for(t in 2:ncol(tbb)) {
  spps = tbb$ESPECIE
  ids = tbb[,t]
  mc = colnames(tbb)[t]
  qvl = which(!is.na(ids))
  dad = meusmarcadores[[mc]]
  aids = ids[qvl]
  dad1 = dad[aids]
  names(dad1) = as.vector(spps)[qvl]
  mids = as.vector(spps)[-qvl]
  mmm = matrix("-",nrow=length(mids),ncol=length(dad1[[1]]),dimnames = list(mids,NULL))
  write.dna(dad1,file='temp.fasta',format='fasta')
  dad1 = read.dna(file='temp.fasta',format='fasta',as.character=T,as.matrix = T)
  dad2 = rbind(dad1,mmm)
  xx = sort(rownames(dad2))
  dad2 = dad2[xx,]
  rownames(dad2) = gsub("[ |.\\]","_",rownames(dad2))
  tt = apply(dad2,2,esogaps)
  dad2 = dad2[,!tt]
  dad2 = as.DNAbin(dad2)
  fn = paste(subpasta,"/",mc,"_PARACOMPARACAO.fasta",sep="")
  write.dna(dad2,file=fn,format='fasta')
}    
file.remove("temp.fasta")
 
#arquivos ja comparativos (mesma ordem mesmo nomes)
arqs = dir(pattern="_PARACOMPARACAO.fasta")
 
#concatenando os dados num unico arquivo
particoes = NULL
final = NULL
pos1=0
esogaps<-function(x) { sum(x=="-")==length(x)}
 
for(a in 1:length(arqs)) {
  marcador = strsplit(arqs[a],"_")[[1]][1]
  dad = read.dna(arqs[a],format='fasta',as.character=T,as.matrix = T) 
  tt = apply(dad,2,esogaps)
  dad = dad[,!tt]
  ll = ncol(dad)
  crset = paste(pos1+1,"-",ll+pos1,sep="")
  names(crset) = marcador
  pos1 = ll+pos1
  particoes = c(particoes,crset)
  #concatena
  if (!is.null(final)) { final = cbind(final,dad)} else { final= dad}
}
final = as.DNAbin(final)
txt = c("BEGIN SETS;",paste("charset ",names(particoes)," = ",particoes,";",sep=""),"END;")
write.nexus.data(final,file='temp.nex', interleaved = F)
tt = scan(file='temp.nex',what='raw',sep="\n")
tt = c(tt,txt)
fn = paste(subpasta,"/meusdadosconcatenados.nex",sep="")
write(tt,sep='\n',file=fn)
 
tbb2 = tbb
for(x in 1:nrow(tbb2)) {
  tbb2[x,-1] = gsub(paste(tbb2[x,1],"_",sep=""),"",tbb2[x,-1])
}
write.table(tbb2,file="tabelacomGENBANKidsDasSequenciasUsadas.csv",sep='\t',row.names = F,na='')