====== Estrutura Filogenética ======
===== Arquivo Phylo =====
A partir da lista de espécies da parcela foi gerada uma hipótese filogenética baseada na supertree do phylomatic [[http://www.phylodiversity.net/phylomatic/phylomatic.html]]. O arquivo gerado foi então lido no R utilizando-se a função read.tree do pacote APE. O arquivo original não consegue ser lido e a mensagem de erro fala que há mais do que uma raiz na árvore... isso é um problema que só consegui resolver tirando Podocarpus do arquivo!Não entendo bem porque, mas parece que o podocarpus + palmeiras + outras não consegue ter uma raiz. Entretanto, a saída gráfica do phylocom gera uma árvore o que significa que o problema pode ser no APE.
==== Atualização 25 setembro de 2009 ====
Não é preciso retirar o //Podocarpus// do arquivo de espécies para ler o arquivo newick no APE. Só é necessário tirar o último "()" e o taxon relacionado.
- {{phylomaticpeic.txt|}}
===== Nossa árvore =====
A partir do arquivo no R gerei a árvore com a função plot (precisa carregar o APE)!
{{peic_tree1.jpg?900|}}
===== Criando a matriz de distância filogenêtica =====
# 1. colocar o path do diretorio phylocom no sistema para rodar sem necessidade de indicar o caminho. Ver manual phylocom
## usando o comando system() para rodar o phylocom no R!
### rodando a distancia filogenetica a partir de um arquivo phylo newick (phylomaticPEIC.cgi) e gravando no arquivo (peic.filodist.txt)
## Deveria ser mais facil segundo o phylodist. Algum problema no comando system() não grava direto o arquivo de saída. Não consegui descobrir o porquê. Fiz um caminho alternativo!
write.table(system("phylocom phydist -f phylomaticPEIC.cgi", intern=TRUE), file="peic.filodist.txt", col.names=FALSE, row.names=FALSE, sep="\t", quote=FALSE )
### lendo o mesmo arquivo para o R! Essa redundância se faz necessária pela minha ignorância de fazê-la diretamente.
peic.filodist=read.table("peic.filodist.txt", row.names=1, header=T)
===== Criando o arquivo "sample" para o filocom =====
Criei uma função file filocom para formatar o arquivo file do filocom a partir do arquivo peic.
# CRIANDO ARQUIVO AMOSTRA (SAMPLE) FILOCOM #####
#### Para o filocom precisa de outro formato de arquivo (3 colunas: nome da parcela, abundancia, codigo spp)
file.filocom=function()
{
choose.files()->f
read.table(f, header=TRUE,sep="\t", as.is=TRUE)->dados
dados.uso=dados[dados$status0=="A" & dados$stemtag==1 & dados$sp !="Indet", ] ## retira dados não validos
file.phylocom=aggregate(dados.uso$stemtag,by=list(dados.uso$quadrat,dados.uso$sp),length)
file.phylocom[,2]<-sub(" +$", "", file.phylocom[,2]) ## retira o espaço no final do string sp se houver!!
file.phylocom[,2]<-gsub(" ", "_", file.phylocom[,2])
return(file.phylocom[order(file.phylocom[,1]),c(1,3,2)])
}
file.filocom()->peicfile.pc
write.table(peicfile.pc, "peicfile_pc.txt",quote=FALSE, col.names=FALSE, row.names=FALSE, sep="\t")
===== Rodando o Phylocom no R! =====
## RODANDO O PHYLOCOM NO R ######################
## MODELO 0 - suffle phylogenetic taxa label
peic.m0.pc=system("phylocom comstruct -f phylomaticPEIC.txt -s peicfile_pc.txt -m 0", intern=TRUE)[-1]
write.table(peic.m0.pc,file="peic_m0.pc", col.names=TRUE, row.names=FALSE, sep="\t", quote=FALSE)
peicm0.res=read.table("peic_m0.pc", header=TRUE, skip=1,sep="\t", as.is=T)
str(peicm0.res)
## investigando os dados significativos
peicm0.res[peicm0.res$NRI>=2 | peicm0.res$NRI<=-2,]
## MODELO 0a - suffle phylogenetic taxa label using abundance in the sample file
peic.m0a.pc=system("phylocom comstruct -f phylomaticPEIC.txt -s peicfile_pc.txt -m 0 -a", intern=TRUE)[-1]
write.table(peic.m0a.pc,file="peic_m0a.pc", col.names=TRUE, row.names=FALSE, sep="\t", quote=FALSE)
peicm0a.res=read.table("peic_m0a.pc", header=TRUE, skip=1,sep="\t", as.is=T)
str(peicm0a.res)
## investigando os dados significativos
peicm0a.res[peicm0a.res$NRI>=2 | peicm0a.res$NRI<=-2,]
## MODELO 3 - swap algoritm
peic.m3.pc=system("phylocom comstruct -f phylomaticPEIC.txt -s peicfile_pc.txt -m 3", intern=TRUE)[-1]
write.table(peic.m3.pc,file="peic_m3.pc", col.names=TRUE, row.names=FALSE, sep="\t", quote=FALSE)
peicm3.res=read.table("peic_m3.pc", header=TRUE, skip=1,sep="\t", as.is=T)
str(peicm3.res)
## investigando os dados significativos
peicm3.res[peicm3.res$NRI>=2 | peicm3.res$NRI<=-2,]
===== Juntando dados de Solo da Parcela =====
## Juntando dados de solo aos dados de estrutura filogenetica
### criando dados de distância cartesianas para as parcelas
letras=rep(LETTERS[1:16],each=16)
numeros=rep(c(paste(rep(0,10),0:9, sep=""),10:15),16)
codigo=paste(letras,numeros,sep="")
parc.x=rep(seq(0,300,by=20),each=16)
parc.y=rep(seq(0,300,by=20),16)
parc.peic=data.frame(parc=codigo, parc.x, parc.y, stringsAsFactors = FALSE)
str(parc.peic)
table(parc.peic$parc%in%peicm3.res$plot)## ver se todos os plot são parc
match(parc.peic$parc,peicm3.res$plot) ## ver a posição de cada elemento primeiro em relação ao segundo
match(parc.peic$parc,peic.solo$quad)
table(match(parc.peic$parc,peic.solo$quad)==1:256) ### confere posições
names(peicm3.res)
names(parc.peic)
names(peic.solo)
m3=data.frame(parc=peicm3.res$plot,parc.x=parc.peic$parc.x, parc.y=parc.peic$parc.y, solo.n=peic.solo$tipo, solo.fator=peic.solo$tipo2, nspp=peicm3.res$ntaxa, MPD=peicm3.res$MPD, NRI=peicm3.res$NRI, MNTD=peicm3.res$MNTD, NTI=peicm3.res$NTI, stringsAsFactors = FALSE)
#### Explorando os dados ##########
par(mfrow=c(2,4))
plot(m3$parc.x, m3$NRI)
plot(m3$parc.x, m3$MPD)
plot(m3$parc.x, m3$MNTD)
plot(m3$parc.x, m3$NTI)
plot(m3$parc.y, m3$NRI)
plot(m3$parc.y, m3$MPD)
plot(m3$parc.y, m3$MNTD)
plot(m3$parc.y, m3$NTI)
boxplot(m3$MPD~m3$solo.fator)
mod1.mpd=lm(MPD~solo.fator,data=m3)
summary(mod1.mpd)
plot(mod1.mpd)
===== Testando modelos de Autoregressão Espacial (SAR) =====
### Modelo GLS
library(nlme)
mod2.mpd=gls(MPD~solo.fator,data=m3)
summary(mod2.mpd)
m3.var=Variogram(mod2.mpd,form=~parc.x + parc.y)
mod3.mpd=update(mod2.mpd, corr=corSpher(200, form=~parc.x + parc.y))
summary(mod3.mpd)
anova(mod2.mpd,mod3.mpd)
### Modelo SAR do spdep
library(spdep)
nb16r=cell2nb(16,16)
nb16q=cell2nb(16,16, type="queen")
str(nb16q)
nblist.peic=nb2listw(nb16q)
## modelo errorsar MPD
mod1.spautolm=spautolm(MPD~solo.fator, data=m3,listw=nblist.peic)
summary(mod1.spautolm)
mod1.errorsar=errorsarlm(MPD~solo.fator, data=m3,listw=nblist.peic)
summary(mod1.errorsar)
mod1.sar=lagsarlm(MPD~solo.fator, data=m3,listw=nblist.peic)
summary(mod1.sar)
anova.(mod1.errorsar,mod1.sar)
### Autocorrelograma Moran I
choose.files()->f ### indicar o caminho da função FuncMoran_ale.r
##"C:\\Users\\Alexandre\\Ale_2009\\AleRfunctions\\AnalisEspacial\\funcMoran_ale.r"
source(f)
dados=matrix(t(m3[,6:10]),ncol=256)
names(m3)[6:10]->row.names(dados)
moran.dados=moran.ale(dados)
moran.nulo.m3=moran.nulo(dados=dados)
# gráfico para MPD
moran.mpd=moran.nulo.m3["MPD",,]
x.seq=seq(20,200, 20)
plot(x.seq, moran.mpd[,1],ylim=c(-0.2,0.4), xaxp=c(20,220,5),ylab="Moran's Index", xlab="Distance (m)")
lines(x.seq,moran.mpd[,2],lty=2)
lines(x.seq,moran.mpd[,3],lty=2)
x.pol=c(x.seq,rev(x.seq))
y.pol=c(moran.mpd[,2],moran.mpd[10:1,3])
polygon(x.pol, y.pol, col="gray", border=NA, density=20)
vf=moran.mpd[,1]moran.mpd[,3]
points(x.seq[vf], moran.mpd[,1][vf], pch=19)
points(x.seq[vf], moran.mpd[,1][vf], cex=1.5)
xy=locator(1)
legend(xy,"Mean Phylogenetic Distance", bty="n")
# gráfico para NRI
moran.nri=moran.nulo.m3["NRI",,]
x.seq=seq(20,200, 20)
plot(x.seq, moran.nri[,1],ylim=c(-0.2,0.4), xaxp=c(20,220,5),ylab="Moran's Index", xlab="Distance (m)")
lines(x.seq,moran.nri[,2],lty=2)
lines(x.seq,moran.nri[,3],lty=2)
x.pol=c(x.seq,rev(x.seq))
y.pol=c(moran.nri[,2],moran.nri[10:1,3])
polygon(x.pol, y.pol, col="gray", border=NA, density=20)
vf=moran.nri[,1]moran.nri[,3]
points(x.seq[vf], moran.nri[,1][vf], pch=19)
points(x.seq[vf], moran.nri[,1][vf], cex=1.5)
xy=locator(1)
legend(xy,"Net Relatedness Index", bty="n")
# gráfico para nspp
moran.nspp=moran.nulo.m3["nspp",,]
x.seq=seq(20,200, 20)
plot(x.seq, moran.nspp[,1],ylim=c(-0.2,0.4), xaxp=c(20,220,5),ylab="Moran's Index", xlab="Distance (m)")
lines(x.seq,moran.nspp[,2],lty=2)
lines(x.seq,moran.nspp[,3],lty=2)
x.pol=c(x.seq,rev(x.seq))
y.pol=c(moran.nspp[,2],moran.nspp[10:1,3])
polygon(x.pol, y.pol, col="gray", border=NA, density=20)
vf=moran.nspp[,1]moran.nspp[,3]
points(x.seq[vf], moran.nspp[,1][vf], pch=19)
points(x.seq[vf], moran.nspp[,1][vf], cex=1.5)
xy=locator(1)
legend(xy,"Net Relatedness Index", bty="n")
===== Gráficos de autocorrelação espacial filogenêtica =====
{{mfd.peic.jpg?400|}}{{nri.peci.jpg?400|}}
{{nssp.peic.jpg?800|}}
===== Investigando os resultados =====
* Parece não haver nenhum sinal filogenêtico associado aos habitats. Ou seja, independente da seleção que possa existir para os habitats de solo, não há assinatura filogenêtica. Os habitats portanto, não restingem as espécies por atributos de linhagens filogenêticas. Há entrentanto, uma maior similaridade florística associada aos habitats, o que nos leva a crer que há uma convergência em atributos selecionados.
* Apesar de não haver uma seleção de linhagens filogenêticas para habitas, há uma forte assinatura filogenética espacial.
* Concluindo: Há uma neutralidade filogenêtica com relação à seleção dos habitats. há seleção de espécies, mas independente da linhagem! YES!
===== Próximos Passos =====
A distância filognetica entre habitats é maior que a distância intra habitats? OU, se estou certo, a distância filogenética entre habitats é similar a intra hábitats (não há diferença significativa). Entre e intra hábitats é o mesmo! OU seja, é neutra a seleção das espécies evolutivamente!