Pipeline de exoma completo - parte 0 (Bioinformática)

Atualização com as partes já publicadas.

Parte 0: referências e conceitos básicos
http://www.estudarcomputacao.com/2014/09/pipeline-de-exoma-completo-parte-0.html

Parte 1: controle de qualidade das sequências
http://www.estudarcomputacao.com/2014/11/pipeline-de-exoma-completo-parte-1.html

Parte 2: alinhamento/mapeamento das sequências
http://www.estudarcomputacao.com/2014/12/pipeline-de-exoma-completo-parte-2.html

----------------------------------------------------------

# Parte 0: referências e conceitos básicos
# Criado em 15/08/2014
# Modificado em 30/09/2014
# Autor: Leandro Lima 

# Oi pessoal. Tudo bem?
# Meu nome é Leandro de Araújo Lima. No momento (fim de 2014) estou concluindo meu
# doutorado em Bioinformática na Universidade de São Paulo (USP), e estou fazendo
# estágio sanduíche no Children's Hospital of Philadelphia (http://www.chop.edu),
# mais especificamente no Center for Applied Genomics (http://caglab.org).

# Você pode me encontrar através do e-mail llima@ime.usp.br ou através da página
# http://www.ime.usp.br/~llima, que tem mais informações sobre mim.

# Este pipeline foi feito para quem não tem experiência com Bioinformática, e
# mesmo para quem não tem experiência com o uso do terminal do Linux (ou Mac,
# por exemplo). Se você tem experiência com o terminal e/ou Bioinformática, mas
# não tem experiência com análise de exomas, este pipeline também é pra você.

# Antes de começar, você precisa ter algumas noções de Biologia Molecular.
# Se você não sabe ou não lembra dois conceitos de DNA, cromossomo e genoma
# (humanos), recomendo que você dê uma relembrada. E, é claro, exoma!
# http://pt.wikipedia.org/wiki/Ácido_desoxirribonucleico
# http://pt.wikipedia.org/wiki/Sequência_de_DNA
# http://pt.wikipedia.org/wiki/Cromossomo
# http://pt.wikipedia.org/wiki/Genoma
# http://pt.wikipedia.org/wiki/Genoma_humano
# http://pt.wikipedia.org/wiki/Exoma

# Não continue sem saber os conceitos acima. Dê pelo menos uma olhada nesses
# artigo da Wikipedia. Alguns deles, inclusive, são bem curtos.

# Recomendo também que você leia alguns post no blog da Genomika.
# Os posts estão em português, e muito bem escritos.
# http://www.genomika.com.br/blog/

# Continuando...

# Todos os comandos a seguir deverão ser executados em um terminal, seja este
# no Linux, BSD ou Mac (por exemplo). Esse código que você está vendo, incluindo
# este texto nas linhas que começam com '#' (jogo-da-velha), são um código de
# uma linguagem própria para ser interpretada e usada no terminal. Esses programas,
# que se comunicam com o Sistema Operacional, são chamados de "shell". O que
# você vai executar são comandos do "bash", que é um tipo específico de shell.
# Saiba mais sobre o "bash": http://pt.wikipedia.org/wiki/Bash (artigo muito legal!)

# Os comandos poderão/deverão ser executados copiando-os e colando-os no terminal.
# Todas as linhas que NÃO começam com o "#" são comandos executáveis. Já as linhas
# que começam com o "#" são comentários, e serão descartadas no terminal.


# OK, agora vamos começar!

# Abra o terminal do Linux e crie um diretório de trabalho. Você pode fazer isso
# usando o comando "mkdir ". Certifique-se de ter uma
# máquina com pelo menos 4GB de memória e 100GB de espaço em disco.

# Caso você não tenha tanto espaço em disco, entre em contato, por favor!
# Posso disponibilizar arquivos menores.

# O "mydir" vai ser o diretório de trabalho. Você pode substituir $PWD por outro
# diretório qualquer ou simplesmente copiar a linha abaixo para tornar o diretório
# em que você está o diretório de trabalho. "mydir" vai passar a ser uma variável
# que vai ser referenciada usando o caracter $ na frente (ou seja, $mydir)

mydir=$PWD

cd $mydir # 'cd' muda o diretório de trabalho para $mydir

# REFERÊNCIAS
mkdir reference # 'mkdir' cria um diretório
cd reference 

# Baixando o genoma humano (versão de montagem: GRCh37, para o 1000 genomas)
# O comando 'wget' serve para baixar um arquivo.
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/2.8/b37/human_g1k_v37.fasta.gz
gunzip human_g1k_v37.fasta.gz # gunzip para descompactar o arquivo .gz

# Saiba mais sobre as versões usadas no 1000 genomas aqui:
# http://www.1000genomes.org/faq/which-reference-assembly-do-you-use

# esse comando vai criar uma variável que representa o caminho para o
# genoma e pode ser acessada pelo terminal
hg19=$mydir/reference/human_g1k_v37.fasta

# Entenda mais sobre as versões do genoma no link abaixo
# http://en.wikipedia.org/wiki/Reference_genome

# Observe que o arquivo tem a extensão 'fasta'
# Em poucas palavras, um arquivo fasta é um arquivo texto com sequências.
# Nesse arquivo do genoma humano, cada sequência é um dos cromossomos
# humanos, e mais algumas sequências de regiões ainda não montadas
# completamente, ou de DNA mitocondrial.
# Você pode ler mais sobre as especificações do arquivo fasta nos links abaixo.
# http://rosalind.info/glossary/fasta-format/
# http://en.wikipedia.org/wiki/FASTA_format

# Para visualizar um arquivo de texto, você pode usar o comando "less". Por exemplo:

less $hg19

# Pressione "q" para sair

# Você vai observar que o arquivo contem, além das sequências, os nomes (rótulos)
# da cada sequência. As linhas com esses nomes têm um sinal de maior (">") logo
# no início. Você pode usar o comando "grep" para ver quais os nomes das sequências
# que estão no arquivo do hg19. Por exemplo:

grep '^>' $hg19
# o sinal circunflexo diz para o grep que você só quer pegar o sinal ">" que esteja
# no início da linha.

# AMOSTRAS
# Agora vamos criar um diretório e baixar os arquivos das sequências que deverão ser mapeadas
cd .. # vai para o diretório pai
mkdir data
cd data

# O "1000 genomas" é um projeto que tem o objetivo de encontrar as variantes mais
# frequentes em diferentes populações. Foram escolhidas mais de 1000 pessoas com
# origens diferentes (afro-americanos, europeus, asiáticos, etc.) e o genoma dessas
# pessoas foi utilizado para encontrar as frequências das variações em cada população,
# além de deixar públicos os dados brutos das sequências.
# Saiba mais sobre o projeto "1000 genomas": http://www.1000genomes.org/

# O sequenciamento das amostras do 1000 genomas foi feito usando o método pareado
# ("paired-end"). Nesse método, a máquina que faz o sequenciamento faz a leitura
# das sequências de duas formas, que é usando cada uma das pontas da sequência
# gerada como início (ou seja, uma leitura do começo para o final, e outra do
# final para o começo). Essa leituras são colocadas em arquivos diferentes,
# mas têm o mesmo nome, só mudando a indicação de se é a sequência 1 (do
# início) ou a 2 (do fim). Depois de baixar as amostras, vamos ver isso.
# Também há um arquivo separado com as sequências cujo par não foi encontrado.
# Saiba mais sobre o sequnciamento pareado ("paired-end"):
# http://en.wikipedia.org/wiki/Paired-end_tag
# http://www.cureffi.org/2012/12/19/forward-and-reverse-reads-in-paired-end-sequencing/
# http://technology.illumina.com/technology/next-generation-sequencing/paired-end-sequencing_assay.html

# Baixando dados das sequências de uma amostra do projeto 1000 genomas
# Sequências do início
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/HG01242/sequence_read/SRR100169_1.filt.fastq.gz
# Sequências do fim
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/HG01242/sequence_read/SRR100169_2.filt.fastq.gz
# Sequências sem par
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/HG01242/sequence_read/SRR100169.filt.fastq.gz

# Vamos criar variáveis para referenciar o caminho completo das sequências
R1=$mydir/data/SRR100169_1.filt.fastq.gz
R2=$mydir/data/SRR100169_2.filt.fastq.gz
SINGLE=$mydir/data/SRR100169.filt.fastq.gz

# Vamos visualizar um dos arquivos "fastq". Como ele está compactado,
# vamos usar o comando "zcat". Como não queremos que todo o conteúdo
# do arquivo seja jogado de uma vez só na tela, vamos usar o comando
# "less" para abrir o conteúdo do arquivo e olhá-lo com calma.
# Para pegar a saída de um programa (ex: zcat) e direcionar para a entrada
# de outro programa (ex: less), usamos o sinal da barra vertical ("|"),
# também conhecido como "pipe". Por exemplo:
zcat $R1 | less # Lembre-se de pressionar "q" para sair

# Observação: os dois comandos acima podem ser resumidos usando-se 
# o comando "zless" (ex: zless $R1).

# O comando "wc" ("word count") serve para contar palavras em um
# arquivo. Mas além de contar palavras, ele também tem uma opção
# para contar linhas (-l). Por exemplo:
wc -l $R1 # isso vai contar o número de linhas do arquivo $R1

# Arquivos "fastq" são uma extensão do arquivo "fasta",
# mas que possibilita incluir a qualidade das bases.
# Essa qualidade é importante para sabermos o quanto devemos confiar
# se cada base é mesmo o que foi lido no processo de sequenciamento.
# No arquivo "fastq", as sequências são mostradas a cada 4 linhas.
# Veja isso usando o comando "head", que mostra as primeiras linhas
# de um arquivo texto (use a opção -n para definir a quantidade,
# que por padrão é de 10 linhas). Nos exemplos abaixo, serão mostradas
# a 8 primeiras linhas (ou seja, 2 sequências de cada arquivo. Compare
# os nomes das sequências de R1 com as sequências de R2.
zcat $R1 | head -n8
zcat $R2 | head -n8

# A primeira linha, que inicia com "@", apresenta o nome da sequência
# (no arquivo "fasta" era o símbolo ">", lembra?).
# A segunda linha apresenta a sequência propriamente dita.
# A terceira linha contem somente um sinal de "+".
# A quarta linha contem as qualidades de cada letra representada na
# 2a linha (cada caractere representa um valor de qualidade da base).
# Entenda mais sobre o formato de arquivo de sequências "fastq"
# http://en.wikipedia.org/wiki/FASTQ_format
# http://maq.sourceforge.net/fastq.shtml
# http://rosalind.info/problems/tfsq/

# Exercícios
# 1 - Quantas sequências há em cada arquivo?
numero_de_linhas=`zcat $R1 | wc -l`
expr $numero_de_linhas / 4
# ou, em uma linha só...
zcat $R1 | echo $((`wc -l`/4))

# O número de linhas/sequências do arquivo R2 deve ser exatamente
# igual ao número de linhas do arquivo R1, pois para cada sequência
# do arquivo R1, há o par correspondente no arquivo R2.

# Infelizmente, não podemos fazer como nos arquivos "fasta",
# buscando pelas linhas que começam com "@" com o grep e 
# depois contá-las. Isso porque o caracter "@" também é
# usado como caracter de qualidade, então ele não só aparece
# no início das linhas com nomes (o mesmo ocorre com o "+").


# 2 - Como imprimir somente os nomes (labels) das sequências (reads)?
zcat $R1 | awk 'NR % 4 == 1' | less
# Observação: o programa "awk" está sendo usado para verificar uma
# condição, que é se a linha sendo lida é a primeira de uma série de 4
# (ou seja, se o número da linha dividido por 4 deixa resto 1)