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?
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?
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.
Use um animal de cada grupo em Taxonomy
Common Tree do NCBI