Estudante de Doutorado em Ecologia no IB-USP.

O meu projeto aborda biogeografia e macroecologia da maior família de serpentes neotropical e propõe elucidar como distribuição, riqueza e diversidade filogenética atual refletem processos históricos e ecológicos, gerando também informação para ações de conservação.

Orientado por Cristiano Nogueira

Link para exercícios resolvidos: Exercícios

A unidade da biodiversidade é geralmente o número de espécies num local, a Riqueza (SR). No entanto as espécies são todas iguais “mas algumas são mais iguais do que outras”, pois refletem processos evolutivos diferentes. A Diversidade Filogenética (PD)¹⁾ visa medir quanta história evolutiva está contida num só lugar, somando o comprimento do ramo entre cada espécie presente na filogenia.

PD é altamente correlacionada com SR e por isso, modelos nulos são necessários para entender se a PD de um local é superior ou inferior ao esperado com as espécies aí presentes. Até à data os modelos nulos disponíveis randomizam espécies da filogenia mantendo constante a SR do local. No entanto, estes modelos não representam um verdadeiro modelo nulo da PD do local pois realisticamente muitas espécies presentes na filogenia não ocorreriam no local por restrições abióticas (e.g. temperatura, humidade, etc), bióticas (tipo de presas) ou biogeográficas (barreiras ou contingência histórica). Desta forma é necessário especificar uma escala regional da qual randomizar as espécies potencialmente presentes à escala local.

A minha proposta é criar uma função que calcule um modelo nulo de PD de um local com composição e SR conhecidas e cuja escala regional possa ser determinada (pela distância ao local focal) para que apenas as espécies nela presente possam ser randomizadas.

Input:

1. Filogenia em formato .tre
2. Dataframe de cada local (1 local/linha) com colunas de:

• local = nome do local (class: string)
• x = longitude (classe: numeric com 2 casas decimais).
• y = latitude (classe: numeric om 2 casas decimais).
• A:k = presença/ausência de cada espécie A,B,k…(classe: binary)

Pseudo-código:

1. require(caper), calcula PD de uma filogenia
2. Cria coluna SR no dataframe, somando todas as presenças de cada local
3. Cria coluna PD no dataframe, usando caper
4. Plota todos os locais num só mapa, com valores de SR como label
5. Calcula distancia euclideana entre os locais utilizando coordenadas x e y
6. Pergunta qual o local focal e que raio de distância o utilizador quer selecionar para escala regional e seleciona todos locais.reg cuja distancia ao local focal seja igual ou menor que o raio
7. Lista todas as espécies species.reg presentes em locais.reg e calcula SR.reg, a soma de todas as espécies únicas
8. Cria objeto simulacao com N NAs.
9. Atribui ao primeiro valor de simulacao o PD de focal
10. Entra em um ciclo for com contador i de 2 até N.
    1. Amostra número de espécies igual aSRdas SR.reg sem reposição e calcula PD para esse subset.
    2. Repete N-1 vezes e coloca valores em simulacao.
11. Cria histograma de simulacao com linha vertical para PD observado em focal
12. Calcula o p-valor para verificar se PD observado é diferente do esperado pelo modelo nulo realista

Verificando os parâmetros:

• N é um número inteiro e maior que 0? Se não, retorna: “N precisa ser um numero inteiro e > 0.”
• focal é um string existent na coluna local? Se não, retorna: “focal precisa ser um local.”
• Informação nas colunas A:k é binária (1 ou 0)? Se não, retorna: “Para cada espécie A:k, colocar 1(presença) ou 0(ausência).”

Output:

• Histograma com a distribuição dos valores de PD do modelo nulo.
• Valor de PD observado plotado no histograma.
• Significância (p-value): probabilidade de valor de PD observado ser diferente ou superior/inferior (dependendo se teste bilateral ou unilateral for selecionado, respetivamente).

Input:

1. Dataframe de coordenadas de ocorrências de espécies provindas de grid

• especie = nome do especie (class: string)
• x = longitude (classe: numeric com 2 casas decimais).
• y = latitude (classe: numeric om 2 casas decimais).

2. Dataframe de coordenadas de habitat provindas de grid

• habitat = nome de classe geográfica (e.g. ecorregião, bioma, fitofisionomia) (class: string)
• x = longitude (classe: numeric com 2 casas decimais).
• y = latitude (classe: numeric om 2 casas decimais).

Pseudo-código:

1. Fazer merge dos dois dataframes num só
2. Cria matrix de contagem para número de pontos de cada espécie “única” em cada habitat
3. Calcula percentagem de ocorrência de cada espécie em cada habitat
4. Calcula o habitat com maior média, dom.habitat
5. Pergunta se utilizador quer fazer merge de um ou mais habitat, para coalescer categorias
6. Se sim, soma colunas de habitat indicadas
7. Calcula novamente o habitat com maior média, dom.habitat
8. Cria objeto simulacao com N+r NAs, com N igual ao número de espécie e r igual ao número de simulações
9. Aloca primeiros N valores em simulacao para dom.habitat de cada especie
10. Entra num ciclo for com contador i de N+1 até r+N.
    1. Simula, com reposição, cada possível par de coordenadas x e y encontrado segundo uma distribuição uniforme
    2. Calcula a percentagem de ocorrência dos pontos simulados para dom.habitat
    3. Guarda valores simulados em simulacao, começando em N+1
11. Cria histograma de simulacao com linha vertical para cada especie, identificada com label
12. Calcula o p-valor

  
regionull.PD=function(tre=NULL,com=NULL,focal=NULL,radius=NULL,nsim=2000,test=NULL,hist=FALSE,col="darkgreen"){
  
  com=as.data.frame(com)#força o objeto de composição de comunidades em dataframe
  com[4:length(com)]=(com[4:length(com)]>0)+0 #converte todas as entradas com números >1 de cada espécie para 1, de forma a que até um objeto com abundâncias seja utilizável
  com$SR=apply(com[4:length(com)],1,sum) #calcula species richness (SR), somando todas as presenças em cada local (desde a 4ª coluna até à penúltima)
  
  xy=data.frame(com[2:3]) #isola coordenadas noutro dataframe
  rownames(xy)=com$local #nomeia cada par de coordenadas com o local correspondente
  dist.com=as.matrix(dist(xy),upper=F) #calcula distância euclideana entre locais e guarda objeto numa matriz simétrica de distâncias
  SR.focal=com[com$local==focal,"SR"] #calcula SR de local focal
  
  if(SR.focal<2)  #se a riqueza do local focal for inferior a 2
    {stop ("local focal deve ter mais do que uma espécie presente!") #pára a função, avisando o utilizador que é impossível calcular o PD nulo da região se o local focal não tiver mais do que uma espécie
  }
  
  else { #se a riqueza do local for superior a 1
  dist.focal=dist.com[focal,] #retorna apenas distâncias entre focal e outros locais
  locais.reg=names(which(dist.focal<radius)) #retorna quais locais estão a menos de "radius" de distância
  if (length(locais.reg)<2) { #se o número de locais regionais for inferior a 2 (não abrangendo nenhum outro local além do focal)
    stop(paste("Nenhum local a menos de", radius, "de distância")) #pára a função e retorna a mensagem
  }
  
  com.reg=com[com$local %in% locais.reg, ]#retorna apenas informação dos locais a menos de "radius" de distância
  focal.range=4:(length(com.reg)-1) #seleciona só desde a 4ª coluna até à penúltima, onde estão as presenças/ausências
  species.focal=names(which(apply(com.reg[com.reg$local == focal,focal.range], 2, sum)>=1)) #retorna quais espécies existem no local focal
  species.reg=names(which(apply(com.reg[focal.range], 2, sum)>=1)) #retorna quais espécies únicas existem a "radius" de distância. Selecionam-se apenas espécies com 1(presença) em pelo menos um local
  
  if(!require(caper)){install.package("caper");library(caper)}#caso pacote "caper" não esteja carregado, instala e carrega pacote "caper". 
  clmat=clade.matrix(tre) #próprio do pacote "caper", manipula a árvore filogenética para análise
  pd.focal=pd.calc(clmat,tip.subset=species.focal) #calcula a distância filogenética (PD) do local focal, o nosso valor observado
  result=rep(NA,nsim) #cria objeto para armazenar nsim simulações
  result[1]=pd.focal #guarda o valor observado na primeira posição
    for(i in 2:nsim) {#loop de nsim vezes e armazena a partir da 2ª posição
    pd.reg=pd.calc(clmat,tip.subset=sample(species.reg,SR.focal)) #cálculo de PD com espécies aleatoriamente escolhidas das presentes na região (a menos de radius de distância) mas respeitando SR do local focal
    
    result[i]=pd.reg #armazena no objeto result
    res=mean(result) - result #calcula resíduos de PD
    res.maior=(sum(res > res[1])/nsim) #calcula probabilidade do resultado ser maior do que o esperado
    res.menor=(sum(res < res[1])/nsim) #calcula probabilidade do resultado ser menor do que o esperado
    res.dif=sum((abs(res) > abs(res[1])))/nsim #calcula probabilidade do resultado ser diferente do esperado
    }
  
  if(hist){ #se hist=TRUE
    hist(res,main="Calculando Diversidade filogenética nula realista",xlab="Resíduos de Diversidade Filogenética",ylab="Frequência",xlim = c(min(res),max(res))) #cria histograma de valores obtidos/esperados
    abline(v = res[1], col=col, lwd=2)#Assinala valor observado no histograma
  }
  if(test=="higher") #se for selecionado o parâmetro "higher" no argumento test
    { 
    if (res.maior>=0.05){ #se p-value for igual ou maior que 0.05, não significante
      res.maior.i="Local não apresenta divergência filogenética" #cria objeto com interpretação de não aceitar H0
    }
    else{ #se p-value for menor que 0.05, significante
      res.maior.i="Local apresenta divergência filogenética" #cria objeto com interpretação de não aceitar H0
    }
    list.maior=list("p-value" = res.maior, "PD do local" = pd.focal[1], "interpretação" = res.maior.i) #cria uma lista com os objetos que o utilizador receberá
    return(list.maior)} #devolve lista
  if(test=="lower") #se for selecionado o parâmetro "lower" no argumento test
  { 
    if (res.menor>=0.05){ #se p-value for igual ou maior que 0.05, não significante
      res.menor.i="Local não apresenta convergência filogenética" #cria objeto com interpretação de aceitar H0
    }
    else{ #se p-value for menor que 0.05, significante
      res.menor.i="Local apresenta convergência filogenética" #cria objeto com interpretação de não aceitar H0
    }
    list.menor=list("p-value" = res.menor, "PD do local" = pd.focal[1], "interpretação" = res.menor.i) #cria uma lista com os objetos que o utilizador receberá
    return(list.menor)} #devolve lista
  if(test=="different") #se for selecionado o parâmetro "different" no argumento test
  { 
    if (res.dif>=0.05){ #se p-value for igual ou maior que 0.05, não significante
      res.dif.i="Local não apresenta padrão filogenético" #cria objeto com interpretação de aceitar H0
    }
    else{ #se p-value for menor que 0.05, significante
      res.dif.i="Local apresenta padrão filogenético" #cria objeto com interpretação de não aceitar H0
    }
    list.dif=list("p-value" = res.dif, "PD do local" = pd.focal[1],"interpretação" = res.dif.i) #cria uma lista com os objetos que o utilizador receberá
    return(list.dif)} #devolve lista
    }
}  

Ficheiro script da função:regionull.pd.r

regionull.PD             package:unknown                R Documentation

~~ Modelo nulo realista de Diversidade Filogenética ~~

Description:

     ~~ Esta função permite testar se a Diversidade Filogenética (PD, Faith 1992) de um local é maior ou menor do que a esperada ao acaso, usando um modelo nulo geograficamente realista para randomizar a presença no local de outras espécies presentes a certa distância.~

Usage:

     ~~ regionull.PD=function(tre,com,focal,radius,nsim,test,hist,col)~~

Arguments:

 tre 
	árvore filogenética a ser escolhida, no formato .tre

 com
	dataframe de locais, suas coordenadas e presença/ausência de espécies

 focal
	comunidade local focal

radius
	distância ao local focal

 nsim
	número de simulações, default=2000

 test
	"different", "lower" ou "higher" conforme padrão,convergência ou divergência filogenética for esperada, respetivamente, default="different"
	
 hist
	TRUE gera histograma, default=FALSE

 col
	cor de valor de PD assinalada no histograma, default="darkgreen"
	


Details:

     ~ A ordem das colunas do dataframe 'com' deverá ser: local, longitude, latitude, seguido de todas as espécies com presença/ausência~
	

Value:

     ~ A função devolve:
	- o p-value do teste escolhido ("different", "lower" ou "higher" ),
	- PD do local 'focal',
	- a possível interpretação do resultado e
	- histograma dos valores de PD simulados para o local, com o PD observado assinalado com a cor definida pelo argumento 'col'


Warning:

	~ A função retorna erro se:
	- Nome das espécies no dataframe 'com' não corresponde às espécies na 'tre' carregada
	- Número de espécies do 'focal' é menor do que 2 (não é possível calcular distância)
	- Número de locais a menos de 'radius' é inferior a 2 (só é selecionado o 'focal' e portanto não é possível simular PD com espécies que ocorram noutros locais)
     

Note:

     ~~Necessário o package "caper"~~
	~Função também aceita matrizes com valores de abundância, mas converte em presença/ausência~
	~Decidi usar este pacote para ler a árvore filogenética ao invés de calcular PD por uma matriz de distâncias pois o utilizador comum geralmente usa uma árvore já completa pois a matriz de distâncias é um jeito menos prático de visualizar relações filogenéticas~ 



Author:

     Filipe Serrano, IB-USP filipe.serrano@usp.br
	http://cnbiogeo.wixsite.com/cristiano-nogueira/people


References:

     ~Faith, D. P. (1992). Conservation evaluation and phylogenetic diversity. Biological conservation, 61(1), 1-10. ~

See Also:

     ~~ Funções: pd.calc(), dist() ~~~

Examples:

     ## Exemplo simples
	tre=(rtree(n=20)) #cria árvore filogenética randómica, com 20 tips ("espécies)
	tips.tre=tre$tip.label #cria objeto com nome das espécies
	com=data.frame(matrix(ncol=3+length(tips.tre),nrow=5)) #cria data.frame vazio

	colnames(com)=c("local","x","y",tips.tre) #atribui nomes às colunas do data.frame
	local=c("jaboticabal","matao","bebedouro","dobrada","ribeirao") # cria objeto com nomes dos locais
	com$local=local #atribui locais à coluna 'local' do data.frame

	x=runif(5,min=0, max=40) #simula valores de longitude de 0 a 40
	y=runif(5,min=0, max=40) #simula valores de latitude de 0 a 40
	com$x=x #atribui valores de longitude à coluna 'x' de com
	com$y=y #atribui valores de latitude à coluna 'x' de com

	for (i in 4:length(com)) #simula presença ou ausência para as espécies e armazena-as nas colunas de espécie do data.frame
		{
  		com[i]=rbinom(n=5,size=1,prob=0.5)
		}		
     ## Correndo a função regionull.PD para verificar convergência filogenética de Jaboticabal, simulando 1500 vezes o valor de PD ao incluir espécies presentes a menos de 80 de distância
	##
	regionull.PD(tre,com,focal="jaboticabal",radius=80,test="lower",hist=T,nsim=1500,col="darkblue") # executa a função


    ## Exemplo mais complexo, com lagartos Anolis
	tre2=read.tree("Anolis.tre") #ficheiro deve estar na pasta do working directory	
	tips.tre2=tre2[4]$tip.label
	df=data.frame(matrix(ncol=3+length(tips.tre2),nrow=5))
	colnames(df)=c("local","x","y",tips.tre2)
	local=c("jaboticabal","matao","bebedouro","dobrada","ribeirao")
	df$local=local
	x=runif(5,min=0, max=40)
	y=runif(5,min=0, max=40)
	df$x=x
	df$y=y

	for (i in 4:lenght(com))
		{
  		df[i]=rbinom(n=5,size=1,prob=0.5)
		}

## Correndo a função regionull.PD para verificar divergência filogenética de Ribeirao simulando 2500 vezes o valor de PD ao incluir espécies presentes a menos de 45 de distância
	##
	regionull.PD(tre=tre2,com=df,focal="ribeirao",radius=45,test="higher",hist=T,nsim=2500,col="firebrick") # executa a função

Example A - random tree Example B - Anolis Impossível carregar ficheiro .tre, fazer botão direito e “salvar como” anolis.tre Ficheiro propriedade de Latin American Macroevolution Workshop, co-taught by Ricardo Betancur, Luke Harmon, & Liam Revell, Página do Workshop