Alinhamentos

Alinhamentos podem ser locais ou globais.

Alinhamento local
Alinhamento local é usado quando queremos identificar, numa coleção de sequências, aquelas que apresentam um alinhamento legal com a nossa sequência pesquisada (Query)
O algorítimo computacional vai utilizar fragmentos da sua sequência (de tamanho = word size, W) e alinhar com as sequências do banco de dados
- ele vai descartar todas sequencias que não coincidirem em pelo menos uma das janelas da sua.
Vamos para a página do BLAST e verificar o word size entrando em "nucleotide blast" e abrindo "Algorithm parameters".
- qual o word size para o BLAST regular (blastn) e para a versão "gulosa" megablast?
- abrindo a aba blastp, qual o word size utilizado para descarte como default?

Matriz de pontuação
Escolhidas as sequencias que compartilham uma janela com a sua, o programa vai alinhar a vizinhança observando se a pontuação é "boa" ou "não compensa mais continuar"
Para pontuar, o programa usa uma matriz
- para comparação de nucleotídeos é usada identidade mas a pontuação é ajustada
- para comparação de proteínas é usada similaridade. A mais usada é BLOSUM62 feita a partir de domínios considerados homólogos mas não muito idênticos para não cair em uma matriz de identidade (máximo de 62% de identidade entre os domínios)

Gaps
Ao estender o alinhamento, o programa tira pontos quando abre um gap e tira pntos quando estende um gap
Quando o escore cai e não se recupera, não compensa mais caminhar, e o alinhamento pára (por isso, local)
O alinhador local não quer chegar ao alinhamento final, ele só quer identificar sequências com alguma coisa (de tamanho significativo) em comum
O fundamento teórico é que a função gênica está quase sempre confinada em domínios contínuos de uma proteína
Quando alinhamos mRNA (cDNA) com genoma e há introns por exemplo, o programa gera múltiplos HSPs

Estatística
E-value é o número de alinhamentos esperados com escore igual ou melhor que o obtido, dado que a base de dados tem um tamanho grande. É uma barreira estatística, caso o valor de e-value seja menor que o limiar estabelecido,as sequencias não se parecem por acaso, mas por serem homólogas.

Tipos de BLAST e quando usar
Vamos alinhas estas sequencias com a base de dados de proteínas denominada "nr" utilizando para isso BLASTp
Anotar uma sequencia envolve propagar a identificação das mais similares, tente fazer isso.
Alguns programas reconhecem trechos que as localizam subcelularmente como Pstort, determine sua possível sublocalização (veja no final do resultado as probabilidades)
BLASTn é a busca onde a query e database são de nucleotídeos. Tente identificar esta sequencia com BLASTn limitando a pesquisa a Homo sapiens: Human genomic + transcript
Compare o resultado com uma busca com BLASTx limitada ao Organism: Homo sapiens. BLASTx traduz a sequencia nucleotídica nas seis possíveis fases de leitura e compara as seis proteínas teóricas com a base de dados de proteínas.
Veja a tabela do código genético: códon é degenerado
Vamos repetir o que foi feito neste trabalho, utilizar tBLASTx com a mesma sequencia e alinhar contra o transcriptoma humano na forma de EST, escolha database EST e limite para Organism: Homo sapiens.
Faltou tBLASTn, quando a query é protéica, mas procura-se na base nucleotídica. Pesquise nas ESTs humanas como acima utilizando a sequencia da telomerase de Euplotes. Para conferir, abra a sequencia do best hit e use em um BLASTx. O que obteve?

Busca de homólogos remotos com PSI-BLAST
PSI-BLAST usa a BLOSUM62 só na primeira busca, é como um BLASTp
Em seguida, faz uma matriz com as sequências que capturou, se numa posição ocorrem D, E, K, R nas proporções 35%, 20%, 15%, 30%, ele dará peso alto para esses resíduos, nessa posição.
Inicie a busca com a sequencia da esfingomielinase da aranha. Mude o Expect threshold para 1e-5 e max seq para 20000 ePSI-BLAST Threshold para 1e-5.
Veja o Taxonomy reports dessa iteração. Algum fungo, alguma bactéria? Run PSI-Blast iteration 2, 3, etc... Corynebacterium tem PLD?

Vamos a Ribeirão Preto fazer BLAST como um bioinformata
ssh cebi15@143.107.223.182
ls
mkdir eusoujacu (faça um diretório pra vc ou use o da outra aula)
cd eusoujacu
ls /home/treinamento/blast_aula/FASTAS
less /home/treinamento/blast_aula/FASTAS/myc.seq - é a sequencia de um fator de transcrição c-myc
less /home/treinamento/blast_aula/454data/tumor.seq - quantas sequencias tem ai?
cat less /home/treinamento/blast_aula/454data/tumor.seq | grep ">" (pare com control C)
cat less /home/treinamento/blast_aula/454data/tumor.seq | grep ">" -c (pergunte ao Dedel se é muita sequencia)
blastall (enter)
blastall -p blastn -i /home/treinamento/blast_aula/FASTAS/myc.seq -d /home/treinamento/blast_aula/454data/tumor.seq -a 20 -F F -e 1e-10 -o saida_normal
blastall -p blastn -i /home/treinamento/blast_aula/FASTAS/myc.seq -d /home/treinamento/blast_aula/454data/tumor.seq -a 20 -F F -e 1e-10 -o saida_tabulada -m 9
blastall -p blastn -i /home/treinamento/blast_aula/FASTAS/myc.seq -d /home/treinamento/blast_aula/454data/tumor.seq -a 20 -F F -e 1e-10 -o saida_tabulada -m 8
Quantos transcritos de c-myc tem na amostra de tumor?
Repita agora com p53.seq. É expressa no tumor?

Anotação com KO
Vamos então verificar se todos os 25 mil genes humanos são expressos no tumor?
less /home/treinamento/blast_aula/CDS/h.sapiens.nuc - veja que são CDS (ATG até TAA, TAG ou TGA)
Contando os CDS: cat /home/treinamento/blast_aula/CDS/h.sapiens.nuc | grep ">" -c
KO será a query (-i) e tumor a database, mas vamos usar a versão gulosa megablast
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

Interpretando o resultado com comandos de linux
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”

Alinhamento global
Serve para quando vc quer comparar totalmente as sequencias mais importantes do seu estudo
E é base para filogenia. Use a ferramenta TaxOnTree com os identificadores desse arquivo e depois desse. Abra com FigTree. No linux baixe o java. Abra outro terminal com control t, entre em Downloads e use java -jar figtree para abrir o programa em ambiente gŕafico. O primeiro arquivo tem sequencias da primeira enzima da bissíntese de Glicina, não essencial
O segundo arquivo tem sequencias de uma enzima da síntese de Valina, Isoleucina e Leucina, essenciais, embora esta enzima permaneça nos animais, a via se perdeu
Teste ambos arquivos com MultiAlign para entender como fica um alinhamento múltiplo global.