# 2_alinhamento.sh # Criado em 15/08/2014 # Modificado em 27/11/2014 # Autor: Leandro Lima <llima@ime.usp.br> # Baixando o programa para alinhar as sequências (bwa) # A versão mais nova pode ser encontrada aqui: # http://sourceforge.net/projects/bio-bwa/files/ cd $mydir/tools # Quando fiz esse tutorial, a versão mais nova era a 0.7.10 wget http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.10.tar.bz2 tar -xvjf bwa-0.7.10.tar.bz2 cd bwa-0.7.10 make # para compilar o programa na sua máquina # Para ficar mais fácil de encontrarmos o bwa, vamos criar uma variável # para acessar a localização dele na sua máquina. Se você souber o que # significa adicionar o programa ao PATH, fique à vontade para fazer isso. BWA=$mydir/tools/bwa-0.7.10/bwa # O 'bwa' precisa que a referência tenha um índice (index), para que as # sequência sejam localizadas mais facilmente. A boa notícia é que o # próprio bwa é usado para criar esse índice. # Comando para criar o arquivo de índices do genoma. # Isso pode levar alguns minutos. $BWA index $hg19 # Instalando o 'samtools', para visualização de arquivos sam/bam # http://sourceforge.net/projects/samtools/files/latest/download cd $mydir/tools wget http://downloads.sourceforge.net/project/samtools/samtools/1.0/samtools-1.0.tar.bz2 tar -xvjf samtools-1.0.tar.bz2 cd samtools-1.0 make SAMTOOLS=$mydir/tools/samtools-1.0/samtools # Agora vamos para a parte do alinhamento. cd $mydir/results mkdir 2_aln cd 2_aln # ALINHAMENTO # Isso também vai levar bastante tempo, dependendo do número de sequências # a serem mapeadas. Basicamente é a parte mais demorada do pipeline. # O 'bwa' gera um arquivo de saída que contem informações sobre todas as # sequências usadas no alinhamento (como: onde foram mapeadas, se foram # mapeadas perfeitamente ou se houve inclusão, remoção ou troca de bases, # e até mesmo se a sequência não foi mapeada em nenhum lugar). Esse arquivo # (SAM, de "Sequence Alignment/Map") é bem grande, já que guarda bastante # informações sobre as sequências. No entanto, ele tem uma versão compactada, # ou binária (BAM), que pode ser gerada usando-se o programa samtools. # As 4 próximas linhas, que estão comentadas, equivalem exatamente às próximas # linhas que vêm neste script. Vamos usar o 'bwa' para mapear as sequências, # e depois pegar a saída do mapeamento (arquivo SAM) e transformar num # arquivo BAM. A diferença é que os comandos não comentados fazem isso sem # precisar guardar o arquivo SAM (grande), economizando espaço em disco. # Isso é muito bom, ainda contando com o fato de que o arquivo SAM pode ser # gerado usando-se o arquivo BAM através do programa 'samtools'. # Saiba mais sobre o 'bwa': # http://bio-bwa.sourceforge.net/bwa.shtml # Saiba mais sobre o 'samtools': # http://samtools.sourceforge.net/samtools.shtml # Você não precisa rodar as 4 linhas abaixo! # $BWA mem -t 4 $hg19 $SINGLE > aln_se.sam # $BWA mem -t 4 $hg19 $R1 $R2 > aln_pe.sam # $SAMTOOLS view -Sb aln_se.sam > aln_se.bam # $SAMTOOLS view -Sb aln_pe.sam > aln_pe.bam # Como explicado anteriormente, os comandos abaixo irão mapear # as sequências usando o 'bwa' (com 4 threads, ou seja, seu # mapeamento será executada em 4 processos paralelamente), # depois a saída desse alinhamento é passada como entrada # para o programa 'samtools' com a opção 'view', que irá pegar # um arquivo SAM (argumento 'S') e gerar um arquivo binário # (argumento 'b'). $BWA mem -t 4 $hg19 $SINGLE | $SAMTOOLS view - -Sb > aln_se.bam $BWA mem -t 4 $hg19 $R1 $R2 | $SAMTOOLS view - -Sb > aln_pe.bam # Depois, iremos juntar as sequências pareadas (método "paired-end") # com as sequências não-pareadas (método "single-end") usando # o programa 'samtools' com a opção 'merge'. $SAMTOOLS merge aln_se_pe.bam aln_se.bam aln_pe.bam # Agora, vamos visualizar algumas informaçõe sobre os arquivos # usando o programa 'samtools' com a opção 'view'. # Primeiramente, vamos entender como funciona um arquivo SAM. # Um arquivo SAM (ou BAM) contem algumas colunas com diversas # informações sobre cada sequência. # Coluna 1: o nome da sequência (que originalmente estava no # arquivo "fastq"). # Coluna 2: contem um número inteiro que é a soma de diversos # valores menores que apresentam características da sequência # e seu alinhamento. # Coluna 3: nome da sequência que serviu de referência para o # alinhamento (nesse caso, um cromossomo). # Coluna 4: a posição, na referência (cromossomo), em que o # alinhamento começou. # Coluna 5: qualidade do alinhamento/mapeamento # Saiba mais sobre as especificações dos arquivos SAM/BAM: # http://samtools.sourceforge.net/samtools.shtml#5 # Quant. de sequências pareadas no arquivo de single-end $SAMTOOLS view -f 0x1 aln_se.bam | wc -l # Quant. de sequências pareadas no arquivo de paired-end $SAMTOOLS view -f 0x1 aln_pe.bam | wc -l # Quant. de sequências pareadas no arquivo com todas as sequências (merge) $SAMTOOLS view -f 0x1 aln_se_pe.bam | wc -l # Quant. de sequências mapeadas $SAMTOOLS view -F 0x4 aln_se.bam | wc -l # Quant. de sequências não-mapeadas $SAMTOOLS view -f 0x4 aln_se.bam | wc -l # Obs: a soma do número de sequências mapeadas com o número de sequências # não-mapeadas sempre deve ser igual ao total de sequências $SAMTOOLS view aln_se.bam | wc -l $SAMTOOLS view -F 0x4 aln_pe.bam | wc -l $SAMTOOLS view -F 0x4 aln_se_pe.bam | wc -l # Observe que o número de seqs. mapeadas do arquivo com todas as seqs. # deve ser a soma das seqs. single-end (189558) com as seqs. paired-end (14797410) # 189558 + 14797410 = 14986968 # Depois de ver essas estatísticas, podemos remover os arquivos # separados, para economizar espaço em disco. rm aln_se.bam aln_pe.bam
Acesse aqui o pipeline completo: http://www.estudarcomputacao.com/2014/09/pipeline-de-exoma-completo-parte-0.html
