# 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 sequencias (bwa)# A versao mais nova pode ser encontrada aqui:# http://sourceforge.net/projects/bio-bwa/files/cd $mydir/tools# Quando fiz esse tutorial, a versao mais nova era a 0.7.10wget http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.10.tar.bz2tar -xvjf bwa-0.7.10.tar.bz2cd bwa-0.7.10make # para compilar o programa na sua maquina# Para ficar mais facil de encontrarmos o bwa, vamos criar uma variavel# para acessar a localizacao dele na sua maquina. Se voce souber o que# significa adicionar o programa ao PATH, fique a vontade para fazer isso.BWA=$mydir/tools/bwa-0.7.10/bwa# O 'bwa' precisa que a referencia tenha um indice (index), para que as# sequencia sejam localizadas mais facilmente. A boa noticia e que o# proprio bwa e usado para criar esse indice.# Comando para criar o arquivo de indices do genoma.# Isso pode levar alguns minutos.$BWA index $hg19# Instalando o 'samtools', para visualizacao de arquivos sam/bam# http://sourceforge.net/projects/samtools/files/latest/downloadcd $mydir/toolswget http://downloads.sourceforge.net/project/samtools/samtools/1.0/samtools-1.0.tar.bz2tar -xvjf samtools-1.0.tar.bz2cd samtools-1.0makeSAMTOOLS=$mydir/tools/samtools-1.0/samtools# Agora vamos para a parte do alinhamento.cd $mydir/resultsmkdir 2_alncd 2_aln# ALINHAMENTO# Isso tambem vai levar bastante tempo, dependendo do numero de sequencias# a serem mapeadas. Basicamente e a parte mais demorada do pipeline.# O 'bwa' gera um arquivo de saida que contem informacoes sobre todas as# sequencias usadas no alinhamento (como: onde foram mapeadas, se foram# mapeadas perfeitamente ou se houve inclusao, remocao ou troca de bases,# e ate mesmo se a sequencia nao foi mapeada em nenhum lugar). Esse arquivo# (SAM, de "Sequence Alignment/Map") e bem grande, ja que guarda bastante# informacoes sobre as sequencias. No entanto, ele tem uma versao compactada,# ou binaria (BAM), que pode ser gerada usando-se o programa samtools.# As 4 proximas linhas, que estao comentadas, equivalem exatamente as proximas# linhas que vem neste script. Vamos usar o 'bwa' para mapear as sequencias,# e depois pegar a saida do mapeamento (arquivo SAM) e transformar num# arquivo BAM. A diferenca e que os comandos nao comentados fazem isso sem# precisar guardar o arquivo SAM (grande), economizando espaco em disco.# Isso e muito bom, ainda contando com o fato de que o arquivo SAM pode ser# gerado usando-se o arquivo BAM atraves 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# Voce nao 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 irao mapear# as sequencias usando o 'bwa' (com 4 threads, ou seja, seu# mapeamento sera executada em 4 processos paralelamente),# depois a saida desse alinhamento e passada como entrada# para o programa 'samtools' com a opcao 'view', que ira pegar# um arquivo SAM (argumento 'S') e gerar um arquivo binario# (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 sequencias pareadas (metodo "paired-end")# com as sequencias nao-pareadas (metodo "single-end") usando# o programa 'samtools' com a opcao 'merge'.$SAMTOOLS merge aln_se_pe.bam aln_se.bam aln_pe.bam# Agora, vamos visualizar algumas