Passo 2

• O exemplo a seguir está baseado no Query ou Busca realizada no GeneBank utilizando a seguinte expressão: rps16[Gene Name] AND RUBIACEAE[Organism], ou seja, todas as Rubiaceae sequenciadas para o marcador de cloroplasto rps16 intron1. Esta busca retornou este arquivo.
• Este passo do exercício é feito no R e exige o pacote Ape. Assumimos que você já tem isso instalado e rodando no seu computador.
• O arquivo sequence.fasta deve estar na sua pasta de trabalho (Working Directory).

#chama o pacote
library(ape)
#le o arquivo
dna = read.dna("sequence.fasta",as.character=T,format='fasta')
#checa o numero de sequencias
length(dna)
#2900 foi o que o GeneBank disse que continha, então está ok
class(dna)
#veja estrutura do objeto
str(dna)
#note que cada elemento da lista é uma sequencia
dna[[1]]
#as sequencias tem nome
names(dna)
#veja o nome das duas primeiras sequencias
names(dna)[idx]
#o nome é comprido e contém muita informação e tem certa lógica (é ruim trabalhar com um nome tão comprido)
#vamos simplificar
#note que cada nome tem lógica: começa com gi, tem coisas separadas por | e depois o nome da especie. Vamos separar tudo isso e colocar numa tabela, depois criar um nome novo mais simples
#pega os nomes
nomes = names(dna)
#cria uma funcao
pegaid <- function(x,posicao=4) {
  xx = strsplit(x,"\\|")[[1]] #quebra o nome onde tem |
  return(xx[posicao]) #retorna a parte do nome da posicao
}
x = nomes[1]
#extrai o identificador do GeneBank (id) - quarta posicao no nome
ll = lapply(nomes,pegaid,posicao=4)
ncbiid = as.vector(ll, mode='character')
 
#pega o nome da especie (posicao > 5)
ll = lapply(nomes,pegaid,posicao=5)
onome = as.vector(ll, mode='character')
 
#as duas primeiras palavras do nome são geralmente especie e genero (para o exercicio vamos ficar só com isso)
 
#cria outra funcao
pegataxon <- function(x, posicao=1) {
  #apaga espacos duplos
  xx = gsub(" "," ",xx)
  xx = gsub(" "," ",xx)
  #quebra o nome onde tem espaco
  xx = strsplit(x," ")[[1]] 
  #apaga objeto nulos
  xx = xx[!xx==""]
  #pega o genero
  return(xx[posicao]) #retorna a parte do nome da posicao
}
 
#pega o genero
ll = lapply(nomes,pegataxon,posicao=2)
genero = as.vector(ll, mode='character')
 
#pega a especie
ll = lapply(nomes,pegataxon,posicao=3)
especie = as.vector(ll, mode='character')
#tira o simbolo de ponto '.' do nome da especie
especie = gsub("\\.","",especie)
 
#cria um novo nome para a sequencia, juntando os dados separados (e usando um separador lógico que nao deve estar no nome)
#limpa do nome se houver
especie = gsub("_","-",especie)
nn = paste(genero,especie,ncbiid,sep="_")
 
#repoem os nomes ao arquivo de dna e salve no seu computador
names(dna) = nn
write.dna(dna,file='rps16_Rubiaceae_205-06-22_ncbi.fasta',format='fasta')
 
#este arquivo está muito grande, vamos selecionar  1 especie de cada gênero (aquela com a sequencia mais longa)
newdad = NULL
id=1
for(g in 1:length(unique(genero))) {
  idx = which(genero==unique(genero)[g] & especie!='sp')
  if (length(idx)>0) {
    #pega a maior sequencia por genero
    lel = NULL
    for(j in 1:length(idx)) {
      lel = c(lel,length(dna[[idx[j]]]))
    }
    umsp = idx[which(lel==max(lel))[1]]
    newdad[[id]] = dna[[umsp]]
    names(newdad)[id] = names(dna)[umsp]
    id = id+1
  }
}
#salva
write.dna(newdad,file='rps16_Rubiaceae_205-06-22_ncbiUmPorGenero.fasta',format='fasta')
 
#ainda está grande
length(newdad)
#vamos diminuir para o exercício e amostrar aleatoriamente 20 generos (nao vai querer fazer isso normalmente)
nz = sample(names(newdad),20,replace=F)
ndd = NULL
for(n in 1:length(nz)) {
  ndd[[nz[n]]] = newdad[[nz[n]]]
}
names(ndd)
write.dna(ndd,file='rps16_Rubiaceae_205-06-22_ncbi_So20.fasta',format='fasta')