Alinhamento de Seqüências Genéticas

1. Alinhamento de Seqüências Genéticas Fevereiro-2011 Daniel Guariz Pinheiro, PhD. Laboratório de Genética Molecular e Bioinformática Departamento de Genética Faculdade de Medicina de Ribeirão Preto Universidade de São Paulo

2. Introdução Seqüências Genéticas

3. Seqüenciamento de DNA Leroy Hood Roche/454 FLX ABI SOLiD ION Torrent • Seqüenciador • semi-automático 2002 2005 2007 2010 1986 2006 2008 2010 Gilbert & Sanger 1986 AppliedBiosystems 1977 • Seqüenciador • automático • comercial Illumina/Solexa GenomeAnalyzer Helicos HeliScope Pacific Biosciences • Métodos para o • seqüenciamento de DNA

4. Nova Geração de Seqüenciadores de DNA ABI 3730xl Roche/454 FLX Illumina/SolexaGA ABI SOLiD Adapted from Richard Wilson, School of Medicine, Washington University, “Sequencing the Cancer Genome”

5. 2008

6. SequenceReadArchive “(…) In mid-September 2010, the SRA contained >500 billion reads consisting of 60 trillion base pairs available for download (…) Almost 80% of the sequencing data are derived from the Illumina GA platform. The SOLiD™ and Roche/454 platforms account for 15% and 5% of submitted base pairs, respectively.(…)” “We’re growing by about 1 Tb/month.” NCBI’s staff scientist Martin Shumway (Leinonen R et. al., 2011)

7. Formato Fasta sequence.fa sequence.qual >SEQUENCE_1 cagtcagcatactcagtcagtcatgcatgctga gtcacttgcatgacgtcatgactgcatgactgc >SEQUENCE_1 1 9 7 15 20 21 16 26 31 37 38 ... 31 13 23 29 31 33 35 30 29 34 ... Extensões: .fa, .fasta, .fna

8. Qualidade O que queremos dizer com qualidade ? >SEQUENCE_1 1 9 7 15 20 21 16 26 31 37 38 ... 31 13 23 29 31 33 35 30 29 34 ...

9. Formato fastq Formato fastq Qualidadecodificadacomo um únicocaracterdatabela ASCII. sequence.fastq @SOLEXA01:1:1:27:1992#0/1 AGTACAAGAGACAGACATTCTTTTTTTTGACACAAG +SOLEXA01:1:1:27:1992#0/1 \FFFMXPYDDHJSUMVUJLPSNFRXZEDLNLHKHIT Extensões: .fastq

10. Introdução ATRIBUINDO Significado às seqüências

11. Há uma referência? • Reseqüenciamento • Existem seqüências produzidas a partir de um genoma/transcriptoma da mesma espécie da amostra ou de uma espécie relacionada que podem ser usadas como referências. Alinhamento com a referência. • Seqüenciamento de novo • Não há seqüências que podem ser usadas como referências. Este tipo de seqüenciamento exigirá uma montagem (assembly) das seqüências, utilizando apenas os dados obtidos desse seqüenciamento. Alinhamento entre as seqüencias geradas, que permitirá a obtenção de um consenso.

12. Seqüenciamento em pares • Seqüenciamento em pares • mate-pair • paired-ends (Korbelet al. , 2007) 36 bp 36 bp 36 bp 36 bp Referência: Referência: ~ 1928 bp a 4928 bp ~ 128 bp a ~428 bp >SOLEXA01:1:1:27:1992#0/1 >SOLEXA02:1:1:11:1992#0/1 >SOLEXA02:1:1:11:1992#0/2 >SOLEXA01:1:1:27:1992#0/2 paired-ends mate-pair

13. Alinhamento de Seqüências EmBioinformática, alinhamento de seqüências é uma forma de dispor as seqüências de DNA, RNA, ouproteínasparaidentificarregiões de similaridadequepodem ser conseqüência de relacionamentosfuncionais, estruturaisourelaçõesevolutivas entre elas.

14. Significado Biológico do Alinhamento de Seqüências • Definição de 3 termos importantes: • identidade: refere-se à fração de aminoácidosounucleotídeosidênticos entre pares de seqüênciasapós um alinhamentodessasseqüências; • similaridade: refere-se à fração de aminoácidosounucleotídeossimilares (com propriedadesfísico-químicassemelhantes – aminoácidosconservados) entre pares de seqüênciasapós um alinhamentodessasseqüências; • homologia: representa uma relação evolutiva entre as seqüências;

15. Identificação das seqüências • Reseqüenciamento • Alinhamento: Conjunto de Seqüências X Seqüências Referências (Ex.: Genoma) >seq1 gcagtcagtcacacatgtca... >seq2 cgcgcatgcgcgtactctat... >seq3 tcgagcatcatcagtcgtca... >seq4 tatgctttatagcgagtcat... ..... >chrX atcacacatgtcacatggtcag ggcatcagtcagtcagtcatgc gcgcgcatgcgcgtactctatc tcatgcgtcagtcatgcatgcg agcagtcatgcatgcatcgcac tgcatcatacgtcatgcatgaa ..... • Objetivos: • - Eliminar as sequência sem hit • - Eliminar as sequência com hits múltiplos (ambiguous) • - Guardar as sequência com hit único (unambiguous)

16. Montagem de seqüências • Seqüenciamento de novo • Alinhamento: Conjunto de Seqüências X Conjunto de Seqüências Consensus : Seq A Seq B Seq C SeqD Seq E Seq F Seq G ACAGTACGACAGTACGACCAGTACGATAGCAGTACGATACGACCGA TCCAGTACGATAGCAGTACGATCAG GCACAGTACGACCAGTACGATACAGGAAC CAGGTACGATACGACGGACGGGG ACAGTACGACAGTACGAAAC GTACGACCAGTACGATACACT AACGACAGTACGAAACGGG TATAGGTACGATACGACGGAC

17. Introdução algoritmos para alinhamento de seqüências

18. Alinhamentos de Seqüências • Alinhamento Global (e.g.Algoritmo de Needleman-Wunsch) • As seqüências envolvidas devem ser alinhadas de um extremo ao outro. Adequado quando as seqüências possuem aproximadamente o mesmo tamanho. Seq X : C A T T A G C A G C C T | . ||||| Seq Y : - A G T A – - A G C - - • Alinhamento Local (e.g. Algoritmo de Smith–Waterman) • Procura-se alinhar apenas as regiões mais similares, independente da localização relativa de cada região. Seq X[4,10]: T A G C A G C | |||| Seq Y[3,7]: T A - - A G C • Alinhamentos (Global/Local) (DNA/Protein) • FASTA (http://fasta.bioch.virginia.edu/fasta_www2/fasta_list2.shtml) • EMBOSS Align (http://www.ebi.ac.uk/Tools/emboss/align/)

19. Problema • Transformar uma seqüência de caracteres em outra: • Operações: • inserção • deleção • substituição • Custo de operação: • Score de substituição • Penalidade para Gaps (inserção/deleção) • Qual é a quantidade de operações mínima ? • Como achar a séries de operações que vai garantir que usamos a quantidade de operações mínima ? Exemplo: Scores: Match: 2 Mismatch (S): -1 Gap(I): -2 Gap(D): -2 ACGT || G-GT Score (4-2-1): 1 2 matches: 4 1 gap: -2 1 mismatch: -1

20. Soluções • Método força bruta (busca exaustiva) • Praticamente inviável • Algoritmos de Programação Dinâmica • Smith-Waterman; Needleman-Wunsch; • SW é um algoritmo para achar o alinhamento mais provável com uma estrutura certa; • Por razões de tempo e espaço, não pode ser usado para alinhamento de sequências de larga escala; • Utilizações de aproximações (heurísticas); • Geralmente, quanto mais rápida for a aproximação, mais distante estará a resposta da solução “correta”;

21. Matriz de Programação Dinâmica GG A • > Score (-2-1): -3 • 1 gap: -2 • 1 mismatch: -1 • > Score(-1-2): -3 • 1 mismatch: -1 • 1 gap: -2 • > Score(-4-2): -6 • 2 gaps: -4 • 1 gap: -2 GG A GG A traceback D(i-1, j-1) + s(xi, yj) (diagonal -> match/mismatch) D(i -1, j) + g (acima -> gap acima) D(i, j -1) + g (esquerda -> gap esquerda) D(i, j) = max Exemplo: Scores: Match: 2 Mismatch (S): -1 Gap(I): -2 Gap(D): -2 ACGT || G-GT Score (4-2-1): 1 2 matches: 4 1 gap: -2 1 mismatch: -1

22. BLAST • Basic Local Alignment Search Tool • http://blast.ncbi.nlm.nih.gov/ • Heurística: dicionário de palavras E-value (S): número de diferentes alinhamentos com scores equivalentes ou melhores que S que são esperados ocorrer ao acaso em buscas em um banco de dados aleatório, do mesmo tamanho, com a mesma composição de bases; QUANTO MENOR... MELHOR!!! NÃO CONFUNDIR COM P-value (probabilidade)

23. BLAT • BLAT—TheBLAST-LikeAlignmentTool • http://genome.ucsc.edu/ • Estruturalmente diferente (BLAST) • Além de outros pontos, o Blat constrói um índice do banco de dado de seqüências (database) (k-mers) e faz as buscas na seqüência a qual se deseja consultar (query); • Blat é mais rápido, porém menos sensível; • Possui código especialmente para lidar com intros em alinhamentos RNA/DNA; • Comumente utilizado para localizar uma determinada seqüência no genoma ou determinar a estrutura de exons de um RNA; • Pode ser utilizado para alinhar seqüências de Roche/454;

24. Alinhamento de seqüências curtas • BLAST/BLAT são lentos demais para alinhar milhões de sequências (Illumina: 35bp-100bp/SOLiD: ) • Premissas: • Não precisamos de um alinhamento sofisticado como SW; • Não precisamos de estatísticas com e-value; • Normalmente, sabemos a quantidade de mismatches máximas que queremos;

25. Alinhamentos baseados em Hashingtable • Idéia dos algoritmos de alinhamentos baseados em hashingtables: read: agtcgtat | ||| || genoma: acggcacgaggaactcgaatctgacgcatgcagtacta Se admitirmos 2 mismatches entre a minha sequência e o genoma. • Se separados em 4 fragmentos, vão existir pelo menos 2 fragmentos sem mismatches !

26. seeds • 6 possibilidades de seeds, com no mínimo 2 fragmentos de match perfeito read: agtcgtat agtc---- 1 --tcgt-- 2 ----gtat 3 ag--gt-- 4 --tc--at 5 ag----at 6

27. Busca nas tabelas hash • São construídas 6 listas das palavras achadas nas leituras; • Para cada 6 possibilidades de palavra no genoma, procurar na lista determinada para ver se existe uma possibilidade de matching; • Como buscar sequências das palavras nas listas de palavras? • Algoritmo de hashing • Possui uma função capaz de transformar uma cadeia de caracteres (string) em valores (índices);

28. Alinhamento de Seqüências com hashing • Softwares • ELAND (Anthony. J. Cox, 2006, unpublished data), • MAQ (Li H et al., 2008) • SOAP (Li R et al., 2008) • Características: • Para detectar mais mismatch, precisamos de mais seeds: • Mais mismatch => mais tempo • Algoritmo mais sofisticado para o alinhamento vai requerer mais tempo: • Indels/gaps => mais tempo • Problemas com hashing: • Memória e tempo • Precisa de CPUs múltiplos com muita memória. • Necessidade de métodos menos “glutões”

29. Burrows–Wheeler • Algoritmo usado normalmente em softwares de compressão (.bzip2) • Em alinhadores de seqüências: • Bowtie (Langmead B et al., 2009) • BWA • BWA-SHORT (Li H. andDurbin R., 2009) • BWA-SW (Li H. andDurbin R., 2010)

30. Bowtie • http://bowtie-bio.sourceforge.net • Burrows-Wheeler; • Reduz a quantidade de memória e de tempo para alinhar sequências curtas; • Podem ser usadas seqüências Illumina e SOLiD • Deficiências: • Não tem garantia de retornar todos os hits com mismatches (exceto com opção --best) • Limite de 3 mismatches (demora mais) • Reads longos reduz a velocidade • Não tem indels

31. LF mapping

32. InexactMatching

33. Introdução Bowtie

34. BowtieIndexBuilder: bowtie-build Usage: bowtie-build [options]* <reference_in> <ebwt_outfile_base>reference_in comma-separated list of files with ref sequencesebwt_outfile_base write Ebwt data to files with this dir/basenameOptions: -f reference files are Fasta (default) -c reference sequences given on cmd line (as <seq_in>) -C/--color build a colorspace index -a/--noauto disable automatic -p/--bmax/--dcv memory-fitting -p/--packed use packed strings internally; slower, uses less mem -B build both letter- and colorspace indexes --bmax <int> max bucket sz for blockwise suffix-array builder --bmaxdivn <int> max bucket sz as divisor of ref len (default: 4) --dcv <int> diff-cover period for blockwise (default: 1024) --nodc disable diff-cover (algorithm becomes quadratic) -r/--noref don't build .3/.4.ebwt (packed reference) portion -3/--justref just build .3/.4.ebwt (packed reference) portion -o/--offrate <int> SA is sampled every 2^offRate BWT chars (default: 5) -t/--ftabchars <int> # of chars consumed in initial lookup (default: 10) --ntoa convert Ns in reference to As --seed <int> seed for random number generator -q/--quiet verbose output (for debugging) -h/--help print detailed description of tool and its options --usage print this usage message --version print version information and quit hg18.1.ebwt hg18.2.ebwt hg18.3.ebwt hg18.4.ebwt hg18.rev.1.ebwt hg18.rev.2.ebwt [/data/indexes]$ bowtie-build /data/hg18.fa hg18 $BOWTIE_INDEXES=“/data/indexes”

35. BowtieIndexInspector: bowtie-inspect Usage: bowtie-inspect [options]* <ebwt_base> <ebwt_base> ebwt filename minus trailing .1.ebwt/.2.ebwt By default, prints FASTA records of the indexed nucleotide sequences to standard out. With -n, just prints names. With -s, just prints a summary of the index parameters and sequences. With -e, preserves colors if applicable.Options: -a/--across <int> Number of characters across in FASTA output (default: 60) -n/--names Print reference sequence names only -s/--summary Print summary incl. ref names, lengths, index properties -e/--ebwt-ref Reconstruct reference from ebwt (slow, preserves colors) -v/--verbose Verbose output (for debugging) -h/--help print detailed description of tool and its options --help print this usage message [/data/indexes]$ bowtie-inspect -s hg18

36. BowtieAligner: bowtie Usage: bowtie [options]* <ebwt> {-1 <m1> -2 <m2> | --12 <r> | <s>} [<hit>] <m1> Comma-separated list of files containing upstream mates (or the sequences themselves, if -c is set) paired with mates in <m2> <m2> Comma-separated list of files containing downstream mates (or the sequences themselves if -c is set) paired with mates in <m1> <r> Comma-separated list of files containing Crossbow-style reads. Can be a mixture of paired and unpaired. Specify "-" for stdin. <s> Comma-separated list of files containing unpaired reads, or the sequences themselves, if -c is set. Specify "-" for stdin. <hit> File to write hits to (default: stdout) Input: -q query input files are FASTQ .fq/.fastq (default) -f query input files are (multi-)FASTA .fa/.mfa -r query input files are raw one-sequence-per-line -c query sequences given on cmd line (as <mates>, <singles>) -C reads and index are in colorspace -Q/--quals <file> QV file(s) corresponding to CSFASTA inputs; use with -f -C --Q1/--Q2 <file> same as -Q, but for mate files 1 and 2 respectively -s/--skip <int> skip the first <int> reads/pairs in the input -u/--qupto <int> stop after first <int> reads/pairs (excl. skipped reads) -5/--trim5 <int> trim <int> bases from 5' (left) end of reads -3/--trim3 <int> trim <int> bases from 3' (right) end of reads --phred33-quals input quals are Phred+33 (default) --phred64-quals input quals are Phred+64 (same as --solexa1.3-quals) --solexa-quals input quals are from GA Pipeline ver. < 1.3 --solexa1.3-quals input quals are from GA Pipeline ver. >= 1.3 --integer-quals qualities are given as space-separated integers (not ASCII)Alignment: -v <int> report end-to-end hits w/ <=v mismatches; ignore qualities or-n/--seedmms <int> max mismatches in seed (can be 0-3, default: -n 2) -e/--maqerr <int> max sum of mismatch quals across alignment for -n (def: 70) -l/--seedlen <int> seed length for -n (default: 28) --nomaqround disable Maq-like quality rounding for -n (nearest 10 <= 30) -I/--minins <int> minimum insert size for paired-end alignment (default: 0) -X/--maxins <int> maximum insert size for paired-end alignment (default: 250) --fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (default: --fr) --nofw/--norc do not align to forward/reverse-complement reference strand --maxbts <int> max # backtracks for -n 2/3 (default: 125, 800 for --best) --pairtries <int> max # attempts to find mate for anchor hit (default: 100) -y/--tryhard try hard to find valid alignments, at the expense of speed --chunkmbs <int> max megabytes of RAM for best-first search frames (def: 64) Reporting: -k <int> report up to <int> good alignments per read (default: 1) -a/--all report all alignments per read (much slower than low -k)-m <int> suppress all alignments if > <int> exist (def: no limit) -M <int> like -m, but reports 1 random hit (MAPQ=0); requires --best --best hits guaranteed best stratum; ties broken by quality --strata hits in sub-optimal strata aren't reported (requires --best)Output: -t/--time print wall-clock time taken by search phases -B/--offbase <int> leftmost ref offset = <int> in bowtie output (default: 0) --quiet print nothing but the alignments --refout write alignments to files refXXXXX.map, 1 map per reference --refidx refer to ref. seqs by 0-based index rather than name --al <fname> write aligned reads/pairs to file(s) <fname> --un <fname> write unaligned reads/pairs to file(s) <fname> --max <fname> write reads/pairs over -m limit to file(s) <fname> --suppress <cols> suppresses given columns (comma-delim'ed) in default output --fullref write entire ref name (default: only up to 1st space)Colorspace: --snpphred <int> Phred penalty for SNP when decoding colorspace (def: 30) or --snpfrac <dec> approx. fraction of SNP bases (e.g. 0.001); sets --snpphred --col-cseq print aligned colorspaceseqs as colors, not decoded bases --col-cqual print original colorspacequals, not decoded quals --col-keepends keep nucleotides at extreme ends of decoded alignmentSAM: -S/--sam write hits in SAM format --mapq <int> default mapping quality (MAPQ) to print for SAM alignments --sam-noheadsupppress header lines (starting with @) for SAM output --sam-nosqsupppress @SQ header lines for SAM output --sam-RG <text> add <text> (usually "lab=value") to @RG line of SAM headerPerformance: -o/--offrate <int> override offrate of index; must be >= index's offrate-p/--threads <int> number of alignment threads to launch (default: 1)--mm use memory-mapped I/O for index; many 'bowtie's can share --shmem use shared mem for index; many 'bowtie's can shareOther: --seed <int> seed for random number generator --verbose verbose output (for debugging) --version print version information and quit -h/--help print this usage message [/data]$ bowtie hg18 > -c "AGGAATTGCGGGAGGAAAATGGGTAGTTAGCTATTT,AGGGCCCATAGCAACAGATTTCTAGCCCCCTGAAGA" > --best --strata --tryhard -m 1

37. Conclusão considerações finais

38. Conclusão • Alinhamento global: Alinhamento de 2 sequências com mesmo tamanho: • Algoritmo de Needleman-Wunsch • Alinhamento local: Alinhamento de 2 seqüências, uma curta e a outra muito mas longa: • Algoritmo de Smith-Waterman • Encontram o alinhamento mais provável; • Lentos para alinhamentos contra o genoma inteiro; • Baseados em um modelo matemático, os outros, são baseados em heurísticas, sem prova formal de obtenção da solução ótima;

39. Conclusão • BLAST: Utiliza heurísticas (k-tuples); • Maior sensibilidade; • Possui estatísticas, o E-value além do Score; • Pode ser usado para Sanger (megablast), mas é muito lento com seqüências Roche/454; • BLAT: Utiliza heurísticas (semelhante ao BLAST - índice do banco de dados na memória) • Blat é mais rápido, porém menos sensível; • Lida melhor com intros em alinhamentos RNA/DNA, bom para determinar estrutura de exons de RNAs; • Pode ser utilizado para alinhar seqüências de Roche/454;

40. Conclusão • Next-GenerationSequenceAlignments • Primeiros programas: Hashing • Illumina e SOLiD; • ELAND (Anthony. J. Cox, 2006, unpublished data), • MAQ (Li H et al., 2008) • SOAP (Li R et al., 2008) • Requerem muita memória; • O nível de sensibilidade depende do programa e das opções; • A partir de 2009: Burrows-Wheeler • Illumina e SOLiD; • Bowtie (Langmead B et al., 2009) • BWA (Li H. andDurbin R., 2009) • Requerem menos memória e são mais rápidos;

41. Conclusão • Novas plataformas de seqüenciamento irão surgir exigindo novos programas de alinhamento; • Não há um programa perfeito para todas as situações; • É importante entender como os programas funcionam e como a configuração pode influenciar os resultados; • Heurística utilizada; • Argumentos; rapidez (tempo/memória) sensibilidade

42. Visualização • IGV (Genome Browser) • http://www.broadinstitute.org/software/igv/home • Formatos de arquivos: • BAM, BED, Birdsuite Files, CBS, CN, Cytoband, FASTA, GCT, genePred, GFF, GISTIC, HDF5, IGV, LOH, MAF, PSL, MUT, RES, SAM, SampleInformation, SEG, SNP, TAB, TDF, TrackLine, TypeLine, WIG

43. Daniel Guariz Pinheiro dgpinheiro@gmail.com

44. Fim Exercício

45. Bowtie • http://lgmb.fmrp.usp.br/~daniel/downloads/cvbioinfo2011/ • cvbioinfo2011_p1.fa • cvbioinfo2011_p2.fa