mySQL

Análise de transcriptoma usando a base de dados Kegg Orthology

BLAST usando 25 mil CDS humanas como query versus 500 mil transcritos de tumor de mama como database:

megablast -i h.sapiens.nuc -d tumor.seq -D 3 -F F -a 20 -p 97 -s 80 -o megakegg

-i        = input
-d        = database
-D 3    = saida tabulada
-F F   = filtro de baixa complexidade desligado (F = False)
-a 20    = use 20 processadores
-p 97    = mínima identidade 97% (sequenciamento pode ter até 3% de erro)
-s 80    = mínimo escore 80
-o       = nome do arquivo de saída

Interprete o resultado com comandos de shell

cat megakegg | awk ‘{print $1}’ | sort | uniq -c | sort -k 1,1 -n -r > resultado

awk = nenhuma condição (aspas) e ação de imprimir coluna 1 (chaves)
sort = ordena as queries (até já estão ordenadas, só pra garantir)
uniq = mostra cada query uma única vez
uniq -c = idem, mas mostra quantas de cada uma haviam
sort -k 1,1 = ordena pela coluna 1, que é o número de ocorrências das queries
sort -n = entende valores como números ao invés de palavras
sort -r = reverso, ordena do maior para o menor

Selecione algum hit e procure a identificação no arquivo h.sapiens.nuc

$cat h.sapiens,nuc | grep “hsa:1234”

Há como fazer de forma global usando banco de dados (MySQL)

Para isso é necessário preparar tabelas a serem utilizadas. Bora aprender um pouco de MySQL?

1) ls /home/treinamento/mysql_aula e veja quantos arquivos

2) crie um diretorio mysql_aula:

    mkdir mysql_aula

3) e entre nele com cd mysql_aula

4) copie o material para usar com cp /home/treinamento/mysql_aula/* . {olhe o ponto!}
    vc pode dar um less em alguns desses arquivos (saia apertando "q")

5) Crie um arquivo tabulado com apenas algumas colunas do blast:

    cat megakegg | awk -v OFS="\t" '($1 != "#") {print $1,$2,$3,$11,$12}' > megakegg_tab

    -v OFS = insere um tabulador entre as colunas selecionadas

    Mas pode ser feito com cut:

    cat megakegg | grep –v “#” | cut –f 1,2,3,11,12 > megakegg_tab

    cat = joga o conteúdo do arquivo na tela
    grep = seleciona linhas contendo algo, mas com -v seleciona linhas que não contêm algo
    cut = como o nome diz, corta, e -f aponta as colunas a pegar
    > = sem aspas, no terminal, faz salvar com o nome indicado após o sinal

6) Importante: estando nessa mesma pasta, onde estão os arquivos, acesse o mySQL com seu username e password (user11 senha11):

    mysql -u <username> -p <ENTER>

7) Informe a senha

8) Listar bases de dados disponíveis é a primeira coisa que se faz ao conectar!

    show databases;

9) Informe ao mysql qual banco de dados vai usar (banco de dados tem o mesmo nome do seu usuário: database11)

    use <database>;

10) Verifique que o banco ainda está vazio (sem tabelas):

    show tables;

    VEJA O PROCEDIMENTO

11) Crie uma tabela para armazenar os dados selecionados do BLAST:

create table result_blast (cds varchar(15), subject varchar(50), identity double(5,2), evalue varchar(10), score int, index cds_idx (cds));

12) Verifique a estrutura da tabela criada:

desc result_blast;

13) Carregue o resultado do blast para a tabela:

load data local infile 'megakegg_tab' into table result_blast;

14) Verifique a tabela criada, um select vale por um ls em mysql!

select * from result_blast limit 10;

- veja que o mesmo arquivo agora se encontra em colunas no banco de dados

15) a consulta anterior utilizando awk |sort|uniq já nos informou quais "hsa" (CDSs) deram mais hits e portanto estão mais presentes nessa amostra de tumor. Que tal utilizar o banco pra isso?

select cds, count(cds) from result_blast group by cds limit 10;

- mais simples não? dizemos que vc deu um count (a)grupando pelos cds

16) Crie uma tabela a partir dessa contagem de CDSs (basta acrescentar create table ao select):

create table hsa_count select cds, count(cds) as hits from result_blast group by cds;

- repare só que vc renomeou counts como "hits"

17) Verifique a nova tabela criada hsa_count:

select * from hsa_count limit 10;

- e veja quantos hsa estão na tabela com select count(*) from hsa_count;

18) Cansado de ver apenas símbolos "hsa", sem saber o que eles são? Crie uma tabela contendo descrições dos genes:

create table hsa_description (cds varchar(15), description varchar(150), index cds_idx (cds));

19) Verifique a tabela criada, sempre é bom:

select * from hsa_description limit 10;

20) Carregue as descrições na tabela recém criada:

load data local infile 'hsa_description' into table hsa_description;

21) Verifique os dados carregados, sempre bom também, vai que tá errado:

select * from hsa_description limit 10;

22) Agora altere a primeira tabela hsa_count para conter também a coluna de descrição do gene:

alter table hsa_count add column description varchar(150);

23) Visualize nova estrutura da tabela, tá bom.. só se quiser:

select * from hsa_count limit 10;

24) Atualize a tabela hsa_count com as descrições dos genes:

update hsa_count, hsa_description set hsa_count.description = hsa_description.description where hsa_count.cds = hsa_description.cds;

- tranquilo? vc igualou o campo description das duas tabelas (pois em uma ele estava vazio) onde os cds eram os mesmos nas duas tabelas.

25) Verifique a tabela agora com as descrições:

select * from hsa_count limit 10;

- ahhh, agora sim!

26) Qual o nome do carinha que da mais hit mesmo?

select * from hsa_count order by hits desc limit 10;

- vc deu um order by, de forma decrescente

Agora vamos relacionar nossos "hsa" a grupos KO ao qual pertencem? O que é KO mesmo?

Vamos simular o uso de uma database subindo o arquivo dump dela!
Imaginando que vc copiou para seu diretorio mysql_aula o arquivo databaseaula.sql (itens 2 a 4)... De fora do mysql fatemos isso:
mysql -u userX -psenhaX databaseX < databaseaula.sql

27) Crie uma tabela que relaciona KOs e CDS:

create table hsa_ko (cds varchar(15), ko varchar(11), hits bigint default 0, index cds_idx (cds), index ko_idx(ko));

desc hsa_ko;

28) Perdido? Não sabe quantas tabelas já criou?

show tables;

29) Carregue os dados na tabela recém criada:

load data local infile 'hsa_ko.list' into table hsa_ko;

30) Inclua a contagem de hits na tabela recém criada, nela só genes que tenham classificação KO vão entrar:

update hsa_ko, hsa_count set hsa_ko.hits = hsa_count.hits where hsa_ko.cds = hsa_count.cds;

- mesma sintaxe de antes, faz o campo hits pegar o valor da outra tabela sempre que o cds for o mesmo

31) Verifique o numero de pares "hsa" que têm KO (inclusive os CDS que não tiveram hits):

select count(*) from hsa_ko;

32) Vamos excluir pares hsa x ko cujo CDS não obtive hits e não nos interessam:

delete from hsa_ko where hits = 0;

33) Verifique o numero de pares hsa x ko restantes:

select count(*) from hsa_ko;

Muito menos né, claro, nem todos CDS que tem KO estão transcritos na amostra

34) A qual KO o danadinho que deu mais hit pertence? Vamos dar uma olhada tambem em um KO especifico, o K00799

select * from hsa_ko order by hits desc limit 10;
select * from hsa_ko where ko='K00799';

35) Agora crie mais uma tabela contendo o número de CDSs e o total de hits para cada ko.

create table ko_hits select ko, count(distinct cds) as total_cds, sum(hits) as total_hits from hsa_ko group by ko;

- vc criou uma tabela com um select, tranquilo, com atenção no ko, contou os cds que ele tem, somou seus hits como total hits, (a)grupando por ko

36) Verifique tabela criada:

select * from ko_hits limit 10;

37) Já sei, mais uma vez cansado de ver identificadores sem saber o que eles são né? Então crie maaaais uma tabela com descrições dos KO:

create table ko_description (ko varchar(11) primary key, description varchar(150));

38) Popule (em mysqelês se fala assim) essa tabela:

load data local infile 'ko_desc' into table ko_description;

39) Verifique os dados:

select * from ko_description limit 10;

40) Adicione coluna de descrição do ko na tabela ko_hits, tá lembrado como né?

alter table ko_hits add column ko_desc varchar(150);

41) Atualize a tabela ko_hits com as descrições:

update ko_hits, ko_description set ko_hits.ko_desc = ko_description.description where ko_hits.ko = ko_description.ko;

- como anteriormente, faz o update de um campo que estava NULL usando a outra tabela como referência, onde os ko são os mesmos, claro

42) Dá uma oiadinha:

select * from ko_hits limit 10;

43) Qual KO tem mais hits?

select * from ko_hits order by total_hits desc limit 10;

44) Pergunte ao Miguelito um KO legal. Efetue a busca LIKE na tabela de KOs:

select * from ko_hits where ko_desc like '%tumor%';

    ou

select * from ko_hits where ko_desc like '%catalase%';

- veja que % é usada como coringa, pode ter qualquer coisa antes, ou depois

45) Pra finalizar, crie uma tabela que relaciona KOs e vias metabólicas!!!

create table KOmap (path varchar(25), ko varchar(25), path_desc varchar(150));

46) Popule a tabela:

load data local infile 'KO2map' into table KOmap;

47) Agora é quebradera, já sabemos usar select, create, show, load, alter, order by, update, limit, where, like, count ... Tudo que um biólogo precisa saber certo? Errado, falta o JOIN!

    Misture dados de tabelas distintas para facilitar a observação dos resultados:
    Aproveita, e mistura vários comandos vai... select, join, order by e limit

select ko_hits.*, KOmap.path, KOmap.path_desc from ko_hits inner join KOmap on ko_hits.ko = KOmap.ko order by total_hits desc limit 10;

Observe a redundância, um KO pertence a diferentes vias metabólicas e elas podem ser visualizadas nessa consulta.

- vc pode adicionar algum "where" usando LIKE (a descrição contendo alguma coisa) para achar coisas nesse JOIN

select ko_hits.*, KOmap.path, KOmap.path_desc from ko_hits inner join KOmap on ko_hits.ko = KOmap.ko where ko_desc like '%tumor%' order by total_hits desc limit 10;

APRENDEU A PERGUNTAR BEM??? Além de mostrar como bioinformatas trabalham com bases de dados, o tutorial ensina quase todos os bons comandos de mySQL (ask well)

Kegg Pathway:

        Pathway tem as relações de vias bioquímicas.. E são muitas!

        Procure em Pathway informações sobre a via Prostate cancer (item 6.2)

        Que tal verificar sistema de reparo de mismatch? Quer os genes humanos? Simples, troque em cima por Homo sapiens e olhe para a direita.

        Ah quer Corynebacterium diphtheriae NTCC? Então selecione e olhe o lado esquerdo!
        Mas que tal Trichomonas vaginalis? Quem é essa cara do lado esquerdo?

        Encontre POU5F1 (Oct4) no Kegg, veja o pathway em que está, E desvende seu grupo de ortólogos (KO).
        Faça o mesmo para IFNG (Interferon Gamma). Verifique sua distribuição taxonômica e compare com TETHERIN ou BST2 examinando seu KO K06731 e estudando sua distribuicao taxonomica.
        Use um animal de cada grupo em Taxonomy Common Tree do NCBI