Rosalind – Data Formats

Como ninguém é uma ilha e existem mais dados do que podemos guardar em planilhas, o que alias é algo que eu acho que é bem ruim na biologia, sempre guardarmos as coisas em tabelas, e planilhas, o que faz com que, quando a base de dados começa a crescer, começamos a sofrer com organização e inconsistências, mesmo com poucas pessoas manipulando os dados, mas essa discussão fica para depois.

Esse problema é sobre como os dados são guardados, e como acessar um banco de dados.
O Genbank é um banco de dados de nucleotídeos (DNA) e possui uma interface para acesso remota para a maioria das linguagens de computador comuns, entre elas R e python.
O python tem um projeto chamado biopython, que traz um pacote com muita coisa para bioinformática, entre elas o pacote Entrez.

Aqui a gente so tem que aprender a olhar, e basicamente ler a documentação.

A gente recebe uma lista de ids

FJ817486 JX069768 JX469983

E tem que imprimir o registro em formato fasta do que tiver a menor sequência nessa lista.

Veja que se a gente digitar FJ817486 na entrade de busca do genbank, vamos cair aqui.

Dai é so contar o tamanho da sequência la em baixo, ou mesmo olhar ali na descrição que tem a quantidade de pares de base sequenciados, e ver quem ta lista tem a menor e imprimir no formato fasta assim:

>gi|408690371|gb|JX469983.1| Zea mays subsp. mays clone UT3343 G2-like transcription factor mRNA, partial cds ATGATGTATCATGCGAAGAATTTTTCTGTGCCCTTTGCTCCGCAGAGGGCACAGGATAATGAGCATGCAA GTAATATTGGAGGTATTGGTGGACCCAACATAAGCAACCCTGCTAATCCTGTAGGAAGTGGGAAACAACG GCTACGGTGGACATCGGATCTTCATAATCGCTTTGTGGATGCCATCGCCCAGCTTGGTGGACCAGACAGA GCTACACCTAAAGGGGTTCTCACTGTGATGGGTGTACCAGGGATCACAATTTATCATGTGAAGAGCCATC TGCAGAAGTATCGCCTTGCAAAGTATATACCCGACTCTCCTGCTGAAGGTTCCAAGGACGAAAAGAAAGA TTCGAGTGATTCCCTCTCGAACACGGATTCGGCACCAGGATTGCAAATCAATGAGGCACTAAAGATGCAA ATGGAGGTTCAGAAGCGACTACATGAGCAACTCGAGGTTCAAAGACAACTGCAACTAAGAATTGAAGCAC AAGGAAGATACTTGCAGATGATCATTGAGGAGCAACAAAAGCTTGGTGGATCAATTAAGGCTTCTGAGGA TCAGAAGCTTTCTGATTCACCTCCAAGCTTAGATGACTACCCAGAGAGCATGCAACCTTCTCCCAAGAAA CCAAGGATAGACGCATTATCACCAGATTCAGAGCGCGATACAACACAACCTGAATTCGAATCCCATTTGA TCGGTCCGTGGGATCACGGCATTGCATTCCCAGTGGAGGAGTTCAAAGCAGGCCCTGCTATGAGCAAGTC A

Claro que faremos isso usando um script do python, que vai ser assim:

from Bio import Entrez
from Bio import SeqIO

Primeiro importamos os pacotes que precisamos, o Entrez é para acessar o banco de dados e o SeqIO é para fazer o parse na saída do query que o Entrez faz.

arquivo = open('/home/augusto/Downloads/rosalind_frmt.txt')
 
ids=[]
for linha in arquivo:
	ids.extend(linha.rstrip().split())

Abrimos o arquivo que tem as chaves de acesso, para buscar no genbank, e colocamos ela dentro de uma lista, veja que o extend porque estamos adicionando na lista ids, se usássemos o append, teríamos listas dentro de listas.

Entrez.email = 'ribas.aca@gmail.com'
handle = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
records = list(SeqIO.parse(handle, 'fasta'))

Então com al ista de ids, pedimos pelo Entrez, para buscar no banco de nucleotídeos, a lista de ids, e retornar tudo no formato fasta, depois damos um parse para deixar num formato de lista, melhor para trabalhar no python.

menor=10**6
indice=-1
for i in range(len(records)):
	if len(records[i].seq) < menor:
		menor=len(records[i].seq)
		indice=i

Agora achamos o que tem o menor tamanho, veja que records é uma lista, que tem um atributo .seq para cada instancia, que é uma string, e com len vemos o tamanho dela, assumimos um número muito grande para o menor, cuidado que potencia no python é com ** e nao ^ como no R, eu errei aqui :), e depois é so achar o menor, podíamos ter considerado o primeiro da lista como o menor também, mas tanto faz.

sequencia=records[indice].seq
lista=[]
while sequencia:
    lista.append(sequencia[:70])
    sequencia = sequencia[70:]

Uma vez que achamos, lembre-se que no records temos uma mega string, então temos que formatar de novo para fasta, então basta quebrar em linhas de 70 caracteres, aqui uma lista com esses pedaços

print '>'+records[indice].description
for linha in lista:
	print linha

E agora é so imprimir, detalhe que eu errei umas 3 vezes esse problema, porque esqueci de adicionar o “>” antes do nome da sequencia, os detalhes são sempre importantes, para fazer as coisas certas afinal.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python
'''
 
'''
 
from Bio import Entrez
from Bio import SeqIO
 
arquivo = open('/home/augusto/Downloads/rosalind_frmt.txt')
 
ids=[]
for linha in arquivo:
	ids.extend(linha.rstrip().split())
 
Entrez.email = 'ribas.aca@gmail.com'
handle = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
records = list(SeqIO.parse(handle, 'fasta'))
 
menor=10**6
indice=-1
for i in range(len(records)):
	if len(records[i].seq) < menor:
		menor=len(records[i].seq)
		indice=i
 
print indice
 
sequencia=records[indice].seq
lista=[]
while sequencia:
    lista.append(sequencia[:70])
    sequencia = sequencia[70:]
 
print '>'+records[indice].description
for linha in lista:
	print linha

Rosalind – Introduction to Pattern Matching

Bem, é chegada a hora de voltar a ver os problemas do Rosalind. Faz um bom tempo que não resolvia nenhum exercício de la.
Então vamos começar com um fácil, ou não.

No exercício Introduction to Pattern Matching, a gente tem que implementar uma arvore de prefixos ou trie para as sequências.
Basicamente essa é uma arvore que permite guardar dados sem ter que repetir prefixos, e tem vários usos, por exemplo, se a gente falar um prefixo, a partir dela é fácil achar quais palavras tem aquele prefixo.

Então a diferença aqui é que ao invés de um alfabeto completo, a gente vai ter apenas o alfabeto das sequencias de dna, que é o ATCG.

A entrada então é uma serie de sequencias

ATAGA ATC GAT

E devemos retornar a árvore trie formado por ela, falando as arestas e qual a letra que ela liga na árvore, sengo que vamos anotar as arestas na ordem que são escritas.

1 2 A 2 3 T 3 4 A 4 5 G 5 6 A 3 7 C 1 8 G 8 9 A 9 10 T

Então a gente vai ser por exemplo, em python, apos um parse, uma lista de strings que são as cadeias

['ATAGA', 'ATC', 'GAT']

E a gente vai ter que construir a arvore trie. Quando eu estudei isso com minha querida professora de estrutura de dados, a gente implementava tudo na linguagem C, usando ponteiros, que de certa forma facilitava as coisas. Isso é possível de ser feito em Python também, usando classes, já que tudo é um ponteiro para objetos, assim como Java, eu acho, pelo menos aqui tem um exemplo legal pelo Robert Sedgewick, que aparece no curso dele de algorítimos do coursera.

Mas lendo um pouco eu vi essa implementação bem legal usando lista e dicionários. Então ao invés de usar ponteiros, a gente usa dicionários.

Bem primeiro a gente cria uma classe em python, para segurar essa estrutura

1
2
3
4
class Trie(object):
    def __init__(self):
        self.contagem = 1
        self.raiz = [self.contagem, {}]

Veja que atributo contagem é para manter a contagem dos nós que estão sendo criados, a gente começa em 1 que é a raiz, dai vamos incrementar ao inserir um nó, e como não temos nada, começamos com um dicionario vazio.
Em python, um dicionário vai ser uma estrutura que a gente cria com essa chaves, que vai ter uma chave e um valor para essa chave, separado por :.
Por exemplo:

{'A':adenina, 'T':timina, 'C':citosina, 'G': guanina}

A chave tem que ser um identificador, como uma letra ou número, pra gente vai ser as letras, e o valor associado com a chave vai ser…

Outra chave, ou o próximo nível da árvore. Legal né.
Poderíamos só usar o dicionario ser não tivéssemos que manter o registro da entrada de um novo nó na árvore, mas como temos que fazer isso, vamos usar uma lista, com o primeiro elemento da lista sendo o valor da ordem, e o segundo o dicionário.

1
2
3
4
5
6
7
    def inserir(self, sequencia):
        node = self.raiz
        for base in sequencia:
            if base not in node[1]:
                self.contagem = self.contagem +1
                node[1][base] = [self.contagem, {}]
            node = node[1][base]

Então agora vamos ver como vamos inserir algo.
Para inserir, criamos a função inserir, que recebe a sequência que desejamos inserir.
Primeiro pegamos a lista que é a raiz, veja que essa lista tem dois elementos, no elemento 0 temos a contagem que não vamos mexer, e no elemento 1 o dicionário.
Dai para cada base na sequência, usando o loop do for, vamos ver se a base é uma chave na lista, veja que usamos o not in node[1], que é o dicionário, então ele ta comparando essa base com o conjunto de chaves, se ele não encontrar, ele aumenta 1 a contagem, e cria uma nova lista, que é um nível abaixo da na árvore, e la adiciona uma lista com a contagem e um dicionário vazio, mas ele adiciona isso no dicionário da raiz, usando a base como chave, e pega esse nova lista feita como raiz para a próxima base.
Se base já estiver no dicionário, a gente não entra nesse if e simplesmente vamos para o dicionário que tem como chave a base em questão

Então para a primeira letra, no inicio a raiz vai ser

[1,{}]

E vai entrar um A, então isso vai ser adicionado ali no dicionario

[1,{A:[2,{}]}]

E assim sucessivamente.

Adicionando todas as sequências de exemplos a gente fica com o seguinte:

[1, {'A': [2, {'T': [3, {'A': [4, {'G': [5, {'A': [6, {}]}]}], 'C': [7, {}]}]}], 'G': [8, {'A': [9, {'T': [10, {}]}]}]}]

E dai para frente, é só formatar a saída, que basicamente é como imprimir a árvore em pré ordem, assim como é feito em árvores primarias. Imprimindo a aresta, e qual a base que gerou ela.

1 2 A 2 3 T 3 4 A 4 5 G 5 6 A 3 7 C 1 8 G 8 9 A 9 10 T

As vezes é meio difícil visualizar essa árvore, mas se a gente olhar o dicionário com uma indentação, da para começar a se ligar no que esta acontecendo.

[1, {'A': [2, {'T': [3, {'A': [4, {'G': [5, {'A': [6, {}] } ] } ], 'C': [7, {} ] } ] } ], 'G': [8, {'A': [9, {'T': [10, {} ]} ]} ]} ]

E para que nunca ouviu falar de arvore TRIE, basicamente a gente está escrevendo nos dados a seguinte árvore.

01

Fiz essa representação usando o pacote igraph, so para ajudar a visualizar, segue o código em R.

1
2
3
4
5
6
library(igraph)
exemplo<-graph.formula( 1--2,2--3,3--4,4--5,5--6,3--7,1--8,8--9,9--10)
V(exemplo)$name<-c("Raiz","A","T","A","G","A","C","G","A","T")
local<-structure(c(215, 184, 187, 190, 191, 195, 148, 242, 243, 246, 
236, 191, 145, 91, 43, 0, 97, 194, 146, 94), .Dim = c(10L, 2L))
plot(exemplo,layout=local)

Lembrando que para criar o local, primeiro eu so plotei o grafo usando o tkplot, depois que eu arrumei os vértices, eu peguei as coordenadas com a função get.coordinates(), dei um dput na saída dela para ter algo que eu pudesse escrever num objeto para colocar aqui. Assim eu pude omitir a parte do tkplot e ainda deixar meu grafo reproduzível.

Bem é isso ai, o script vai estar la no repositório recologia, ele so insere, porque não precisamos remover nem modificar nada para esse exercício, acho que ele vai ser mais de introdução para comparação de sequencias, que não não cheguei em exercícios de alinhamento no rosalind. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
class Trie(object):
    def __init__(self):
        self.contagem = 1
        self.raiz = [self.contagem, {}]
 
    def inserir(self, sequencia):
        node = self.raiz
        for base in sequencia:
            if base not in node[1]:
                self.contagem = self.contagem +1
                node[1][base] = [self.contagem, {}]
            node = node[1][base]
 
    def imprime(self):
        print self.raiz
 
def formata_recursivo(node):
    for base, node2 in node[1].iteritems():
        print node[0], node2[0], base
        formata_recursivo(node2)
 
if __name__ == '__main__':
    entrada = open("trie.in").readlines()
    dados=[]
 
    for linha in entrada:
        dados.append(linha.strip())
 
    print dados
 
    trie_sequencias = Trie()
 
    for sequencia in dados:
        trie_sequencias.inserir(sequencia)
 
    trie_sequencias.imprime()
 
    formata_recursivo(trie_sequencias.raiz)

Banco de Dados – Usando MySQL do R

Manga-Guide-to-Databases

Um banco de dados é uma coleção organizada de dados. Os dados são organizados em modelos que tentam representar a realidade de uma forma que o processamento(adição, remoção e consulta de dados) se torna mais eficiente, além de uma facilidade de guardar.

Para manter esses dados, a gente pode usar sistemas específicos para exercer todas essas funções da melhor forma possível, muitos deles são de uso geral, e existe a possibilidade de migrar entre eles sem começar tudo do zero. Acho que um bem comum para sistemas não tão grandes é o MySQL.

Mas primeiro, porque diabos estamos preocupados com isso? Porque dados em biologia tendem a desaparecer muito rápido, além de serem normalmente desorganizados. Nesse artigo aqui da revista CELL os autores falam como a cada ano da publicação de um artigo a chance de os dados sumirem aumenta em 17%. Na maioria dos editais de pesquisa, ou projetos, não existe uma preocupação com isso dentro da biologia. Alguma áreas, como em biologia molecular, isso ja está mais organizado, ao forçar todos os autores a publicarem seus dados numa base publica, o genbank. Mas em ecologia, as coisas não são bem assim.
Temos algumas revistas, como a Ecology da ESA que publica base de dados, mas eu acho que são casos a parte.
Mas além disso, eu sempre achei que uma maior preocupação com como os dados são guardados, pode ser muito útil mesmo localmente, ja que o ecossistema da biologia são normalmente laboratórios.
Muitos laboratórios, das mais diversas áreas na biologia, é comum você ver muitas pessoas coletando dados muitos similares(no sentido de todo mundo trabalhar no mesmo ambiente por exemplo), mas guardando em planilhas pessoais, cada um guardando de um jeito. Tudo bem que não existe nada de errado com isso, não é meu objetivo reclamar de como outras pessoas fazem seu trabalho. Mas existe algumas possibilidades que talvez sejam perdidas por isso.

Por exemplo, se cada pessoa guarda os dados do seu jeito, temos muitas dificuldades as vezes de conseguir unir esses dados depois, para fazer coisas maiores.
E continuando nesse pensamento, muitas pessoas vem e vão dos laboratórios, fazendo monografias, dissertações e teses, se todos os dados são guardados num banco de dados bem pensado e organizadinho, posteriormente o laboratório ganha uma base de dados, que pode por si só responder muitas perguntas, e além disso pode ajudar a começar a responder outras, já que dados ao longo do tempo são bem complexos de se conseguir.

Mas antes de se preocupar em montar um banco de dados, vamos ver como seria olhar para esses dados do R, so para ter um gostinho da idéia. Para tal vamos ver um pequeno exemplo usando o pacote RMySQL, que serve para acessar banco de dados administrasdos com o MySQL.

Primeiro, SQL é uma linguagem estruturada para “falar” com o banco de dados, que serve para criar eles, acessar os dados que você quer, alterar, enfim todas as atividades, e ela é diferente do programa que administra o banco de dados(mas todos tentam implementar SQL do mesmo jeito), mais ou menos isso.

O MySQL é um programa gerenciador de banco de dados que implementa essa linguagem, e podemos montar um servidor com vários banco de dados, cada banco de dados com seu próprio schema.

Mas vamos ver como funciona para acessar o banco de dados da UCSC de bioinformática.

Basicamente vão ser três atividades, a gente tem que se conectar ao banco de dados, pra isso a gente tem que ter o endereço dele, e um usuário, que pode ou não ter senha, e o gerenciador pode regular o que um usuário pode fazer. Nesse caso tudo esta explicado nessa pagina aqui para esse banco de dados em questão.

> ucscDb <- dbConnect(MySQL(),user="genome",host="genome-mysql.cse.ucsc.edu") > resultado <- dbGetQuery(ucscDb,"show databases;") > head(resultado) Database 1 information_schema 2 ailMel1 3 allMis1 4 anoCar1 5 anoCar2 6 anoGam1 > dbDisconnect(ucscDb) [1] TRUE

Aqui o primeiro comando serve para conseguir acesso ao banco de dados (usando a internet), veja que fornecemos um endereço, o usuário publico do banco de dados e o primeiro argumento é mais complicado, mas é o driver para fazer a conexão. Isso tudo nos retorna um objeto do R que é a conexão. Com ela agora a gente pode “fazer perguntas” para o banco de dados e recuperar os dados. Isso é o que fizemos com o dbGetQuery, a gente perguntou quais bancos de dados tem la, com o “show databases;”, e depois olhamos a saída, que é uma lista de bases que podemos acessar.

Certo, mas onde estão os dados, estão dentro dessas bases de dados, podemos escolher uma e acessar, por exemplo, se pegarmos a base de dados “hg19”, a gente acessa ele assim.

> hg19 <- dbConnect(MySQL(),user="genome", db="hg19",host="genome-mysql.cse.ucsc.edu") allTables <- dbListTables(hg19) > length(allTables) > [1] 11006 > allTables[1:5] [1] "HInv" "HInvGeneMrna" "acembly" "acemblyClass" "acemblyPep"

Como a gente desconectou la em cima, a gente conecta novamente, e podemos listar as tabelas que contém os dados. É um bocado de tabelas, com informações específicas. Podemos tentar olhar quais os atributos, os campos, colunas de uma tabela dessas assim:

> dbListFields(hg19,"affyU133Plus2") [1] "bin" "matches" "misMatches" "repMatches" "nCount" "qNumInsert" [7] "qBaseInsert" "tNumInsert" "tBaseInsert" "strand" "qName" "qSize" [13] "qStart" "qEnd" "tName" "tSize" "tStart" "tEnd" [19] "blockCount" "blockSizes" "qStarts" "tStarts"

Agora o mais legal é que o SQL é uma linguagem de busca, então podemos mais que olhar uma tabela, podemos pegar os dados de forma organizada, e até fazer buscas usando funções do SQL, por exemplo, podemos mandar um query para falar qual o tamanho de uma tabela, o número de linhas da seguinte forma.

> dbGetQuery(hg19, "select count(*) from affyU133Plus2") count(*) 1 58463

Bem aqui a gente ve que essa tabela tem realmente muitas linhas, isso implica que leva um tempo para baixar ela, mas temos um comando para baixar ela inteirinha como um dataframe no R, e tendo ela na memória no R, é so fazer o que queremos.

> affyData <- dbReadTable(hg19, "affyU133Plus2") > head(affyData) bin matches misMatches repMatches nCount qNumInsert qBaseInsert tNumInsert tBaseInsert strand 1 585 530 4 0 23 3 41 3 898 - 2 585 3355 17 0 109 9 67 9 11621 - 3 585 4156 14 0 83 16 18 2 93 - 4 585 4667 9 0 68 21 42 3 5743 - 5 585 5180 14 0 167 10 38 1 29 - 6 585 468 5 0 14 0 0 0 0 - qName qSize qStart qEnd tName tSize tStart tEnd blockCount 1 225995_x_at 637 5 603 chr1 249250621 14361 15816 5 2 225035_x_at 3635 0 3548 chr1 249250621 14381 29483 17 3 226340_x_at 4318 3 4274 chr1 249250621 14399 18745 18 4 1557034_s_at 4834 48 4834 chr1 249250621 14406 24893 23 5 231811_at 5399 0 5399 chr1 249250621 19688 25078 11 6 236841_at 487 0 487 chr1 249250621 27542 28029 1 . . . tStarts 1 14361,14454,14599,14968,15795, 2 14381,14454,14969,15075,15240,15543,15903,16104,16853,17054,17232,17492,17914,17988,18267,24736,29320,

Eu suprimi um pouco dos resultados, mas aqui a gente tem a tabela inteirinha num dataframe no R. Pronto para usar como quisermos. Mas o SQL vai muito além disso, ele permite comandos para selecionarmos exatamente o que queremos dos dados, além do que essas tabelas são montadas de forma estruturada, de forma que conseguimos recuperar diretamente os dados estruturadinhos.

Mas mais um exemplo colocando uma cláusula para pegar todas as colunas da tabela affyU133Plus2 onde a coluna misMatches estão entre 1 e 3.

> query <- dbSendQuery(hg19, "select * from affyU133Plus2 where misMatches between 1 and 3") > affyMis <- fetch(query) > affyMis[1:10,1:3] bin matches misMatches 1 585 723 3 2 586 740 2 3 73 986 3 4 587 741 1 5 587 985 1 6 587 938 1 7 588 784 3 8 588 1397 3 9 589 483 1 10 589 938 1

E daqui podemos analisar os dados

> quantile(affyMis$misMatches) 0% 25% 50% 75% 100% 1 1 2 2 3 Lembrando que como temos uma conexão com o banco de dados, sempre temos que fechar ela depois. > dbClearResult(query) [1] TRUE > dbDisconnect(hg19) [1] TRUE

Mas imaginem que legal, primeiro que teríamos uma base de dados disponível no laboratório que trabalhos direto, pessoas vem e vão, mas os dados sempre vão estar la. Provavelmente mais seguros, porque é mais fácil fazer uma cópia de segurança de um banco de dados que de mil planilha espalhadas nas mãos de muitas pessoas. Outra coisa é que para montar um banco de dados você precisa de um schema, uma estrutura de dados, isso pode nos forçar a pensar melhor sobre o que fazemos, sobre como coletamos dados, quais dados coletamos e porque coletamos, e pessoas que estiverem começando algum trabalho, podem se beneficiar sobre as caneladas passadas que foram corrigidas. Além do que assim, cedo ou tarde é fácil ter uma base de dados grande, se ela for alimentada por todas as pessoas que passam pelo laboratório. Além de que, quem passa pelo laboratório tem respaldo de dados para escrever projetos, o que facilitaria justificar quantidades, número de amostras ou pensar a quantidade de dados que vamos achar, ou se algum projeto é viável.

Outra coisa muito legal, é tornar mais fácil tudo ser reproduzível, só a gente guardar um script do R com os querys do banco de dados e analises realizadas (só controlar datas da busca, já que a base pode crescer sempre) e qualquer analise fica na mão.

Bem eu acho uma idéia bem legal, já que ninguém perde, porque todo mundo tem que coletar seus dados para concluir seja la o que for que estiver fazendo, mas a ciência ganha se tudo ficar guardadinho, e puder ser re-analisado conjugando dados dos vários trabalhos feitos anteriormente, e quem sabendo abrindo a possibilidade de responder coisas mais legais ainda.

Eu vi esse exemplo no curso de R do coursera chamado Getting and Cleaning Data, que mostra muitas outras formas de dados, e é legal para ver como a preocupação com esse tipo de coisa vem aumentando, com a tecnologia fazendo a gente cada vez conseguir mais dados, é importante conseguir guardar tudo direito, de forma estruturada para não se afogar em informação, mas banco de dados é uma área bem complexa, mas muito legal, e acho que não distante para começar empreitadas para guardar dados nos diversos laboratórios que cada um fica, mesmo disciplinas de pós-graduação no campo tendem a gerar muitos dados que ficariam muito legais se guardados num banco de dados bonitinho. Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

install.packages("RMySQL")
 
library(RMySQL)
 
ucscDb <- dbConnect(MySQL(),user="genome",host="genome-mysql.cse.ucsc.edu")
resultado <- dbGetQuery(ucscDb,"show databases;")
head(resultado)
dbDisconnect(ucscDb)
 
 
#Connecting to hg19 and listing tables
hg19 <- dbConnect(MySQL(),user="genome", db="hg19",host="genome-mysql.cse.ucsc.edu")
allTables <- dbListTables(hg19)
length(allTables)
allTables[1:5]
 
#Get dimensions of a specific table
dbListFields(hg19,"affyU133Plus2")
dbGetQuery(hg19, "select count(*) from affyU133Plus2")
 
#Read from the table
affyData <- dbReadTable(hg19, "affyU133Plus2")
head(affyData)
 
#Select a specific subset
query <- dbSendQuery(hg19, "select * from affyU133Plus2 where misMatches between 1 and 3")
affyMis <- fetch(query)
affyMis[1:10,1:3]
quantile(affyMis$misMatches)
 
dbClearResult(query)
dbDisconnect(hg19)

Teoria neutra na evolução e na ecologia

A teoria neutra, na evolução, começou em 1960 com Motto Kimura, basicamente diz que a frequência de alguns alelos podem variar por acaso devido a esses alelos terem um sucesso reprodutivo, um fitness neutro, ambos tem o mesmo valor de fitness.

Bastante simples, mas com implicações fortes. Lembrando que alelos são as formas de expressão de um gene.
Se a gente pensar na natureza e na espécie modelo mais comum em genética, a mosquinha Drosophila melanogaster que tiver, uma mutação, por exemplo, asas curtas.

08
08

Uma mosca, na natureza, sem a capacidade de voar está ferrada, ela vai morrer de fome. Ela vai ter um fitness muito menor que a mosca com asas normais e a frequência desse gene vai tender a sumir, isso é conhecido como seleção natural.

Agora talvez, outra mutações não tenham impacto direto no fitness. Estou supondo aqui, a titulo de exemplo, ja que eu não entendo de Drosophilas, mas a cor dos olhos, se não houver seleção natural qualquer, as fêmeas não verem moscas de olhos vermelhos como mais bonitas que as de olhos laranjas, o que aconteceria?

08
08

Bem essa mutação poderia ser considerada neutra, ou seja ela num melhora nem piora o fitness, mas isso não implica em nenhum momento que se você acompanhar essa população, a frequência desse alelo sera constante. Isso porque temos uma população, na natureza sempre finita, a a frequência é sempre o número de indivíduos que exibe uma característica, o alelo em questão dividido pelo tamanho total da população.

Nesse caso que temos 2 alelos em um gene (situação fictícia, eu estou assumindo que essa mutação e neutra e não existe mais variantes, sei que tem, mas é só para o exemplo ser fácil de entender), os alelos são, olhos vermelhos e olhos laranjas.

Mas esse exemplo é só para visualizar a situação, a ideia é que se esse gene é neutro, existe um processo onde se você começa cada alelo presente em 50% da população, 50 indivíduos, vamos supor, a qualquer momento que alguém usar um mata moscas aleatoriamente e mate, por cagada, 9 moscas de olhos vermelhos e 1 mosca da olhos laranjas, sobrara 41 moscas de olhos vermelhos a 49 de olhos laranjas, e se cada uma tem sei la, 1 filho, na próxima geração, veremos 45,5% de moscas de olhos vermelhos e 54.4% de moscas de olhos laranjas.

O problema é que quando eu bato o mata moscas, eu tento mata moscas, eu não fico olhando qual mosca eu vou matar, e pior que isso, nenhum tipo de moscas é mais rápido, ou mais ágil para fugir nesse exemplo, o que gera uma situação onde essa mudança na frequência dos alelos é completamente imprevisível, pra que direção ela vai.

Essa é um processo conhecido como deriva genética, ou genetic Drift, um efeito da teoria neutra.

Existem várias formas de simular a deriva genética, lembrando que a herança pode ser variável mesmo na gente. Lembre-se que mitocôndrias a gente so ganha da nossa mãe enquanto para genes nucleares alelo vem do nosso pai e outro alelo vem da nossa mãe, e tem crossing-over. Se começarmos a pensar muito nisso a gente não vai pra frente. Mas vamos usar aqui só pra visualizar o modelo do Wright–Fisher, que considera gerações discretas, indivíduos diploides, e cada nova geração tem um tamanho fixo, basicamente cada nova geração a gente pega um número aleatória de indivíduos da frequência da geração anterior.

01

Veja que cada linha é uma simulação, no eixo x temos as gerações, e no eixo x estamos representando a frequência de um alelo, lembrando que as duas frequências vão dar 1, então se um alelo tem 100% o outro tem 0, e vice versa.

Veja que em alguns casos, laranja e verde por exemplo, esse alelo domina, e o outro some, em outras simulações, ele some e só fica o outro alelo na população.

Apesar de imprevisível, a gente pode dizer que a deriva genética acaba sempre mudando as frequências de alelos, ainda mais, quase sempre sumiu com algum dos alelos, ela tirou alguém da jogada, e assim diminuiu a diversidade genética, já que ter uma população com 2 alelos é mais diversa que uma população com 1 alelo, seja ele qual for.

Certo, mas ai vendo tudo isso, um cara chamado Stephen P. Hubbell pensou se isso não poderia acontecer em outra escala, na escala de espécies. Será que existe uma teoria neutra para riqueza de espécies? Então ele pensou na teoria neutra unificada da biodiversidade e biogeografia (the unified neutral theory of biodiversity and biogeography), mas vamos dizer teoria neutra que é menor, ninguém tem saco de escrever esse nome por muito tempo.
A sua importância é que ela possibilita mais que muitos modelos de comunidade, ela faz predições claras em muitos níveis de organização, tanto de processos evolucionários como ecológicos.

Comunidades e genes são computacionalmente e matematicamente relacionadas, geralmente usando modelos iguais aos usados na genética, logo as lições que aprendemos sobre a deriva genética se aplicam a modelos de comunidade. A dinâmica ecológica numa escala local é conhecida de como deriva ecológica, e populações mudam ao acaso, devido a processos estocásticos como o nosso mata-moscas.

Hubbell propôs a teoria neutra como um modelo nulo para a dinâmica de indivíduos, para explicar a abundância de árvores tropicais ao na costa rica, barro colorado. No pacote vegan, se você abrir o pacote e digitar data(BCI), você pode ver parte dos dados que o Hubbell trabalhava, as abundâncias na floresta da reserva do Barro Colorado na costa rica, alias, como curiosidade, se você digitar data(), sem nenhum argumento, o R lista todos os data-set de exemplos disponíveis, de todos os pacotes abertos.

Certo, mas agora o controverso é que o Hubbell propôs, diferente do que comumente associamos a modelos nulos (que são como as coisas acontecem ao acaso normalmente, para teste de hipoteses), que é realmente assim que as coisas acontecem na dinâmica de comunidades.

Mas isso é controverso? Porque?

Porque estamos constantemente avaliando fatores que influenciam na riqueza e diversidade de espécies, e vem um cara e fala que tudo isso num é assim, que quem define quais espécies estão num local é o acaso. Realmente uma afirmação que vai contra muitos trabalhos que estamos acostumados a ver.

Mas se a gente não for muito extremista (e começar a achar que tudo tem que ser ou 100% estocásticos ou 100% determinísticos, sem acordos), a palavra chave para a teoria neutra do Hubbell é metacomunidade, uma coleção de comunidades similares conectadas por dispersão.

A metacomunidade é povoada inteiramente por indivíduos que são funcionalmente idênticos, segundo a teoria neutra. A teoria neutra fica então uma teoria de dinâmica de indivíduos, modelando nascimentos, morte, migração e mutação desses indivíduos, sendo a riqueza de espécies uma medida derivada desses parametros. Assumimos que dentro de uma guilda, como a de árvores do final de sucessão em floresta tropical, as espécies são essencialmente neutras, todo mundo tem um fitness igual, equivalente. Isso significa que probabilidades de nascimento, morte, mutação e migração são próximas ou equivalentes para todos os indivíduos. Assim, mudanças nos tamanhos das populações ocorrem devido a “random walks”, eventos estocásticos. Como a gente viu nas frequências gênicas. A gente pode simular praticamente igual fizemos ao gene, mas vamos usar aqui o pacote untb, que já tem essa simulação pronta.

02

Mas calma, isso não implica necessariamente na ausência de competição ou outros efeitos indiretos mediados pela dependência de densidade. A competição é vista como difusa e igual entre indivíduos.
Efeitos dependentes da densidade aparecem através de formas especificas no número total de indivíduos da comunidade, ou como características dos indivíduos relacionado as probabilidades de natalidade, mortalidade e especiação.

Um paradoxo gerado na teoria neutra de Hubbell é que na ausência de migração ou mutação, a diversidade declina para zero, riqueza para um, uma monodominância. Uma caminha aleatória (random walk), devido a equivalência de fitness, vai eventualmente resultar na perda de todas as espécies, exceto uma. Entretanto a perda de diversidade em uma comunidade local é esperada de forma muito, muito devagar, e é contra-balanceada por imigração e especiação. Assim, espécies não aumentam deterministicamente quando raras, isso faz o conceito de coexistência diferente do conceito que a gente vê nos modelos de Lotka-Volterra. Coexistência no caso da teoria neutra é um efeito estocástico balanceado pelo aparecimento de novas especies localmente. Lembrando que estamos pensando dentro de uma guilda, um nicho, não estamos comparando gramas com árvores aqui, ou onças com aranhas.

Mas se todas as comunidades fazem um random walk para a monodominâncias, como a diversidade é mantida em qualquer região em particular? Dois fatores mantém espécies em qualquer comunidade local. Primeiro, imigração para a comunidade local de indivíduos da metacomunidade. Mesmo com cada comunidade sofrendo um random walk para a monodominância, cada comunidade pode vir a ser dominada por uma espécies diferente, formando o pool de espécies dessa metacomunidade, já que todas as comunidades apresentam fitness diferentes. Assim comunidades separadas são preditas de se tornarem dominadas por espécies diferentes, e essas diferenças entre espécies locais ajudam a manter a diversidade na paisagem ocupada pela metacomunidade.

Um segundo fator, numa escala de tempo muito maior, é a especiação que dentro da metacomunidade mantém a biodiversidade. Mutação e consequentemente, especiação prove uma forma de introdução de diversidade na metacomunidade. Random walks em direção a extinção em comunidades que são tao lentas que são contrabalanceadas por especiação.

E agora temos uma pequena ideia do que se trada a teoria neutra da ecologia, e como ela veio da teoria neutra da genética. As simulações foi mais pra fazer alguma figurinha que realmente explicar algo, mas post sem figurinha é chato, e ainda as figuras foram feitas no ggplot2, que faz figuras bem bonitas, mas que até hoje eu não aprendi a usar direito. o script está no repositório do github, além de aqui em baixo, e precisamos definitivamente de muitos posts ainda pra entender a teoria do Hubbell, até onde ela vai, e como estimar parâmetros em metacomunidades para confrontar dados com teoria nesse caso aqui, mas paciência é uma virtude, sempre.

Referencia:
M Henry H Stevens 2011 A primer of ecology with R. Springer

library(untb)
library(ggplot2)
library(reshape)
library(RColorBrewer)
set.seed(123)
 
###################################################
### Deriva Genética
###################################################
#Original - http://statisticalrecipes.blogspot.com.br/2012/02/simulating-genetic-drift.html
 
#Parametros da Simulação de Deriva Genetica
N = 20           # Número de individuos diploides
N.chrom = 2*N    # Número de cromossomos
N.gen = 100      # Número de gerações
N.sim = 10       # Número de simulações
p = 0.2
q = 1-p
 
# Simulation
saida <- matrix(0,nrow=N.gen,ncol=N.sim,dimnames = list(paste("Geraçao",1:N.gen),paste("S",1:N.sim)))
saida[1,] = rep(N.chrom*p,N.sim) # Inicializando o numero de alelos na primeira geração
for(j in 1:N.sim){
  for(i in 2:N.gen){
    saida[i,j]<-rbinom(1,N.chrom,prob=saida[i-1,j]/N.chrom)
    }
  }
#Transformando em frequencia
saida<-data.frame(saida/N.chrom)
 
# Modificando os dados para a estrutura que da para plotar
sim_data <- melt(saida)
ggplot(sim_data, aes(x = rep(c(1:100), N.sim), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva genética") +
    xlab("Geração") + ylab("Frequencia alélica") +
    ylim(0,1) +
    labs(colour = "Simulações")
#ggsave("01.jpg")
 
###################################################
### Deriva Ecologica
###################################################
saida <- untb(rep(1:10, 90), prob=0, gens=5000, D=9, keep=TRUE)
 
sim_data <- melt(data.frame(species.table(saida)))
 
ggplot(sim_data, aes(x = rep(c(1:5000), 10), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva populacional") +
    xlab("Geração") + ylab("Abundância da espécie") +
    labs(colour = "Simulações")
#ggsave("02.jpg")

Crescimento Populacional e a Seleção Natural

A teoria da evolução é talvez a teoria mais importante para a biologia. Ela unificou muitos ramos científicos atribuídos a biologia que antes dela eram praticamente áreas distintas, como fisiologia, taxonomia, comportamento, etc.

Temos a frase celebre do nosso colega Theodosius Dobzhansky.

Nada em biologia faz sentido se não for à luz da evolução

Mas será que a gente realmente entende como a evolução funciona?

Quero dizer, provavelmente, todo mundo tem uma noção básica. O que a gente aprende na escola é tranquilo, mas o problema é que muitas vezes nos lembramos do exemplo básico e fazemos generalizações muito amplas, ou pior ainda, generalizações curtas e não conseguimos ver o processo como ele realmente ocorre.

Na verdade a teoria da evolução é muito mais do que vamos ver aqui, mas vamos olhar como funciona o negócio de variação no número de descendentes produzidos. Quem não lembra do exemplo da mariposa e da revolução industrial. A Biston betularia. Essa ai é aquela mariposa que existia de dois tipo, ou dois fenótipos, como preferir.

Tinhamos a versão, ou fenótipo mais claro.
claro

E a versão, ou fenótipo, mais escuro.
escuro

Ai antes da poluição, a mais clarinha era comum, e a escura era muito rara, com o acumulo de fuligem nos troncos das arvores, a clarinha começou a não ser muito boa para se camuflar, a escura, da cor da fuligem, começou a se dar bem e o cenário se inverteu, com a clarinhas passando a ser a mais rara e a escura sendo a mais comum.

Daqui pra frente vamos chamar de fenótipos sempre. Se a gente fosse la no campo e conta-se todas as mariposas existentes, poderíamos dizer que existe uma quantidade N de mariposas, essa quantidade sendo independente do fenótipo, que é a cor para as mariposas.

E elas vão vivendo, se reproduzindo a uma certa taxa, e bem isso nos faz lembrar da equação de crescimento populacional.

N_t = N_0 \cdot r^t

Lembre-se que estamos falando do crescimento discreto ai em cima, não confunda com a formula N_t = N_0 \cdot e^{r \cdot t} , a gente já discutiu sobre crescimento discreto versus crescimento contínuo aqui, de uma olhada para relembrar e ver a diferença.

Mas continuando, então temos nosso crescimento populacional, vamos ficar com ele por enquanto para facilitar o raciocinio, se começar a meter aquele número e de Euler vai ficar tudo mais difícil.

N_t = N_0 \cdot r^t

So que essa equação é valida para toda a população dessa determinada espécie, todos os indivíduos de uma região que compõem uma população de uma certa espécie, nossa mariposinha por exemplo.

Vamos pensar nos fenótipos claro e escuro como fenótipo A e fenótipo a, cada um deles então vai ter uma frequência na população, que chamaremos de f_A e f_a e cada um deles vai ter sua taxa de crescimento R, r_A e r_a, que também é conhecido como o fitness absoluto (Absolute fitness), que a princípio vamos considerar como igual a taxa de crescimento da população como um todo, escrevendo isso de forma mais clara, vamos assumir que r_A = r_a = R.

Primeiro, vamos pensar que o número inicial de indivíduos A seria igual a sua frequência no total inicial da população.

N_{A,0} = N_0 \cdot f_A

Certo, se a população inicial é N_0, sabemos que uma fração dela é de mariposas claras, a fração é f_A, por isso se multiplicar o total pela fração N_0 \cdot f_A da N_{A,0}

O mesmo vai ser valido para a população de a,

N_{a,0} = N_0 \cdot f_a

Lembrando que f_A + f_a = 1, tem que ser o total da população.

Vamos em frente, A e a são a população total dessa espécie.

N_0 = N_0 \cdot f_A + N_0 \cdot f_a

Sem muita mágica até aqui.
Certo, agora como o crescimento vai funcionar aqui?
Muito simples, para a próxima geração (vamos pensar apenas na diferença de uma geração para a próxima) temos a quantidade inicial vezes a taxa de crescimento, so que por cada um dos fenótipos.

Para A temos

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a

Mas como r_A = r_a = r, o fitness dos dois são iguais, então poderíamos simplesmente dizer.

N_{A,1} = N_0 \cdot f_A \cdot r

e

N_{a,1} = N_0 \cdot f_a \cdot r

Então podemos ver que a população total no tempo 1 N_1 vai ser

N_1 = N_0 \cdot f_A \cdot r + N_0 \cdot f_a \cdot r

Veja que o termo N_0, que é a população total inicial e o termo rocorrem nas duas multiplicações, então podemos simplificar isso.

N_1 = N_0 \cdot r \cdot (f_A + f_a)

Lembre-se que f_A + f_a = 1 já que são frequência, e esses dois A e a tem que ser o todo.

Então apos uma geração a frequência de A será a quantidade de indivíduos A pelo total de indivíduos da população, a frequência f_A será o seguinte.

{f_{A,1}} = {{N_0 \cdot f_{A,0} \cdot r}\over{N_0 \cdot r}}

Como é tudo uma multiplicação podemos cortar os N_0 e r e assim ficamos somente com.

f_{A,1}= f_{A,0}

A frequência na próxima geração não mudou nada.

Mas é isso que esperamos não? Porque nesse exemplo r_A = r_a = r, ou seja os dois fenótipos eram igualmente bons, mas e se r_A \ne r_a, um deles pode ser melhor que o outro, sei la r_A > r_a ou r_A < r_a, mas ambos os casos r_A \ne r_a, como fica esse mesmo esquema?

A principio, a próxima geração vai ter a mesma formula do crescimento.

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a,

Ou ainda, generalizando para o tempo t, temos

N_{A,t} = N_0 \cdot f_A \cdot {r_A}^t

E para a temos

N_{a,t} = N_0 \cdot f_a \cdot {r_a}^t

Ou seja a população N_t será

N_{t} = N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t

Agora vamos fazer aquela conta denovo, após t gerações, a frequência f_A vai ser a quantidade de individuos A sobre o total N

{f_{A,t}} = {N_0 \cdot f_A \cdot {r_A}^t \over{N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t}}

Agora nos podemos avaliar a frequência de a depois de t gerações assim como fizemos antes, mas podemos simplificar um pouco essa equação. Podemos começar tirando esse N_0 em cima e em baixo na divisão.

A frequência f_A vai ser o total dela sobre o tamanho total da população.
{f_{A,t}} = {f_A \cdot {r_A}^t \over{f_A \cdot {r_A}^t + f_a \cdot {r_a}^t}}

Agora podemos dividir em cima e em baixo por r_A e ficamos com o seguinte

{f_{A,t}} = {{f_A}\over{{f_A} + {f_a \cdot {{r_a}^t\over{{r_A}^t}}}}}

E finalmente temos a nossa equação para acompanhar como a frequência gênica muda ao longo das gerações.
Se a gente pensar um pouco, quando a razão {r_a}^t\over{{r_A}^t} é 1, ambos os fenótipos são igualmente bons, essa razão da um, e teremos  {f_A}\over{f_A + f_a \cdot 1}, como  f_A + f_a = 1, teremos  f_A dividido por 1 e nada nunca muda. Agora qualquer desvio na razão {r_a}^t\over{{r_A}^t} para mais ou para menos que 1, teremos mudanças.

Mas de qualquer forma, vamos jogar alguns números e fazer uma simulação no R.

Precisamos estabelecer algumas condições iniciais, e pensar por quantas gerações queremos simular.

# ************************************************** # Condições iniciais # ************************************************** N0 <- 1000 # População inicial rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a" fA <- 0.3 # Frequência do alelo "A" max_gen <- 20 # Número de gerações a simular

Veja que baseado naquelas quantidades, o resto dos ingredientes são derivados

# ************************************************** # Calculando variáveis derivadas # ************************************************** fa <- 1.0 - fA # Frequência do alelo "a" NA_0 <- N0 * fA # População inicial do alelo "A" Na_0 <- N0 * fa # População inicial do alelo "a"

Mas vamos la, criamos uma função para com aquelas condições iniciais e usando as equações que estabelecemos acima para ver o que acontece.

evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}

E testamos nossa funçãozinha.

> evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10) Nt NA_t Na_t fA_t fa_t [1,] 1000.000 300.0000 700.000 0.3 0.7 [2,] 1440.000 432.0000 1008.000 0.3 0.7 [3,] 1728.000 518.4000 1209.600 0.3 0.7 [4,] 2073.600 622.0800 1451.520 0.3 0.7 [5,] 2488.320 746.4960 1741.824 0.3 0.7 [6,] 2985.984 895.7952 2090.189 0.3 0.7 [7,] 3583.181 1074.9542 2508.227 0.3 0.7 [8,] 4299.817 1289.9451 3009.872 0.3 0.7 [9,] 5159.780 1547.9341 3611.846 0.3 0.7 [10,] 6191.736 1857.5209 4334.215 0.3 0.7 [11,] 7430.084 2229.0251 5201.059 0.3 0.7

Temos como entradas as quantidades iniciais, que repetimos na primeira linha, a geração zero, e temos na saída o tamanho da população total como primeira coluna, como ela vai mudando sendo cada geração uma linha, como cada fenótipo muda, bem como suas frequências ao longo das gerações, mas como olhar para uma matriz de números é meio chato, vamos representar esses dados em dois gráficos, um de como está a população e outro de como estão as frequências.

01

Certo, bem simples não, a população está crescendo exponencialmente, e as frequências genicas estão constantes, como previmos, porque as taxas de crescimento iniciais que colocamos para essa simulação são iguais.

rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a"

Agora lembre-se que nem tudo na vida são flores, se tivermos taxas de crescimento abaixo de 1, a população estará diminuindo.

Vamos pegar esse caso aqui.

#N0 = 1000 rate_A = 0.7 rate_a = 0.7 fA = 0.3

02

Mas ainda temos as frequências dos fenótipos iguais, aqui ta todo mundo se ferrando igual. Agora quando a gente diz que está havendo evolução é quando as frequência começam a mudar, quando aparece um novo mutante que derrepente se da melhor, se alguém tem a taxa de crescimento maior que sua contraparte, vemos o seguinte.

#N0 = 1000 rate_A = 2.0 rate_a = 1.5 fA = 0.02

03

Veja que aqui, nenhum dos fenótipos está necessáriamente se ferrando, só temos um fenótipo que é melhor que o outro, mas ninguém está diminuindo em quantidade ao longo do tempo, a população está crescendo como um todo, todos os fenótipos estão crescendo, apenas um mais rápido que o outro, mas isso faz com que a frequência dos fenótipos mude. Quando eu digo que as vezes a evolução pode não ser tão simples, porque como no exemplo das mariposas, a gente não pensa nesse caso assim, a gente sempre pensa que se alguém é pior, ele tem que estar se ferrando, como no próximo exemplo

#N0 = 1000 rate_A = 1.2 rate_a = 0.9 fA = 0.02

04

Essa é a situação que a gente tem na cabeça para as mariposas, pelo menos a primeira que vinha na minha, se a gente pensar na linha azul como as mariposas brancas e vermelho as mariposas escuras, a gente quer que alguém seja muito comum e alguém raro, ai esse cara mais comum, as mariposas brancas, começam a se dar mal por causa da revolução industrial, as mariposas escuras, que antes eram raras e ruim, começam a se dar bem, e temos uma inversão de quem é mais raro e quem é mais comum, mas vemos que a população decaiu um pouco e depois voltou a crescer. Mas como vimos, a evolução pode acontecer mesmo quando todo mundo está se dando bem, com a população crescendo, e além disso ainda, podemos pensar em outra situação.

#N0 = 10000 rate_A = 0.8 rate_a = 0.6 fA = 0.02

05

E é isso ai, ambos os fenótipos estão se ferrando, ambos estão diminuindo sua frequência. Cada vez que você olha a mata, você ve menos e menos daquela espécia que você observa, mas isso não quer dizer que a evolução não possa estar agindo ali. Veja que aqui não é quem se ferra e quem se da bem, aqui a situação é quem se ferra menos, quem pode perdurar aos tempos difíceis, talvez essa espécie seja extinta, mas talvez não, e quando os tempos mudarem, ela pode voltar a crescer, mas as frequências dos fenótipos estão alteradas já na próxima largada para o crescimento populacional.

Mas, não é a toa que nos humanos somos péssimos geradores de números aleatórios, pelo menos eu, sempre tendia a pensar em evolução, no caso de uma espécie que tem um fenótipo se tornando mais raro, porque a população dele está diminuindo, e outro se tornando mais comum porque a população dele está aumentando, mas esteja a população estacionada, aumentando ou diminuindo, as frequências de alelos, fenótipos, snp, como queria pensar, podem estar mudando, e todos esses cenários tem que estar na nossas cabeças.

So que a primeira coisa que vem na minha cabeça é um fenótipo se dando bem e outro se dando mal. Mas existem muitos mais possibilidades. Então é preciso se policiar para não tomar conclusões precipitadas, e conduzir uma análise de dados que de possibilidade de todos os casos serem avaliados corretamente.

Alias, para teste em sistemas mais complexos, sempre esses teste envolvem simulações com parâmetros aleatórios, porque muitas vezes o que podemos pensar não contempla todas as possibilidades, pensa em quantas árvores podemos desenhar com n espécies, eu nunca pensei que existissem tantas possibilidades assim, e é fácil se perder nas contas. Mas ao tentar algumas simulações, passar tudo para um modelo e ir tentando várias possibilidades de parâmetros, podemos talvez ver aquilo que nunca tinha passado pela nossa cabeça.

Bem é isso ai, as equações la do começo do post eu peguei de um curso do coursera, o Computational Molecular Evolution, ótimo curso alias, uma excelente seleção de tópicos para a ementa e conceitos bem complicados são explicados de maneira bem simples, pena que a maior parte pratica é usando o paup, mas agora é ir adaptando tudo pro R com as ferramentas que tivermos a mãos. Scripts no repositório recologia e aqui em baixo, e até o próximo post.

# **************************************************
# Condições iniciais
# **************************************************
 
N0      <- 1000 # População inicial
rate_A  <- 1.2  # Taxa de crescimento do alelo "A"
rate_a  <- 1.2  # Taxa de crescimento do alelo "a"
fA      <- 0.3  # Frequência do alelo "A"
max_gen <- 20   # Número de gerações a simular
 
# **************************************************
# Calculando variáveis derivadas
# **************************************************
 
fa   <- 1.0 - fA  # Frequência do alelo "a"
NA_0 <- N0 * fA   # População inicial do alelo "A"
Na_0 <- N0 * fa   # Polulação inicial do alelo "a"
 
# **************************************************
# Simulação
# **************************************************
 
evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}
 
evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10)
 
saida<-evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 50     rate_A = 1.2  rate_a = 1.2  fA = 0.3
saida<-evosimu(N0=50,rate_A=1.2,rate_a=1.2,fA=0.3,max_gen=20)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 0.7  rate_a = 0.7  fA = 0.3
saida<-evosimu(N0=100,rate_A=0.7,rate_a=0.7,fA=0.3,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topright",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 2.0  rate_a = 1.5  fA = 0.02
saida<-evosimu(N0=1000,rate_A=2.0,rate_a=1.5,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 1.2  rate_a = 0.9  fA = 0.02
saida<-evosimu(N0=1000,rate_A=1.2,rate_a=0.9,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 10000  rate_A = 0.8  rate_a = 0.6  fA = 0.02
saida<-evosimu(N0=10000,rate_A=0.8,rate_a=0.6,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")