AULA PRÁTICA BIOLOGIA MOLECULAR COMPUTACIONAL​
: EXPLORANDO DADOS NGS BMC Monitores: Lucas Silva e Julio Sosa E:​
​
[email protected], ​
[email protected] 17 setembro de 2015 Professor : João Carlos Setubal Objetivos Ferramentas Comandos Básicos Protocolo Utilizar as ferramentas básicas de alinhamento e montagem de transcritoma para obter os genes diferencialmente expressos entre duas amostras. ●
●
●
●
●
●
●
●
fastqc tophat2 cufflinks cuffmerge cuffdiff trimmomatic homer2 linux Para trabalhar com nossos dados, vamos precisar saber alguns comandos básicos do Linux. Podem procurar mais informação no site: http://wiki.ubuntu­br.org/ComandosBasicos 1) Vamos comecar pelo básico, certifique­se de que a pasta atual é : /home/biomolcomp/ Para isso digite o comando $ pwd 2) Crie um diretório contendo o nome da dupla no seguinte formato: $ mkdir aluno1_aluno2 3) Entre no diretório criado $ cd aluno1_aluno2 4) Crie um link simbólico dos arquivos /home/biomolcomp/data/adrenal_1.fastq /home/biomolcomp//data/adrenal_2.fastq /home/biomolcomp/data/brain_1.fastq /home/biomolcomp/data/brain_2.fastq O comando para a criação do link simbólico deve ser: $ ln ­s /home/biomolcomp/data/adrenal_1.fastq . $ ln ­s /home/biomolcomp//data/adrenal_2.fastq . $ ln ­s /home/biomolcomp/data/brain_1.fastq . $ ln ­s /home/biomolcomp/data/brain_2.fastq . 5) Agora é hora ver na prática como um arquivo fasta é. Utlize: $ less adrenal_1.fastq ​
(para sair digite a letra Q) Faça o mesmo para os outros arquivos. $ less adrenal_2.fastq​
(para sair digite a letra Q) (Nota: mais informação do formato FASTQ em https://en.wikipedia.org/wiki/FASTQ_format​
) 6) Seria interessante conseguir saber o numero de sequências totais nos arquivos, não é? Temos algo para isso. O comando ​
wc​
nome_do_arquivo (word count) $ wc adrenal_2.fastq O problema aqui é que o comando conta o número de linhas totais. Mas podemos utilizar uma união com o comando ​
grep​
(para procurar apenas o cabeçalho das reads. E só depois fazer a contagem.) Vamos lá! $ grep ‘@ERR’ adrenal_1.fastq | wc 7) Vamos fazer outros testes agora. Abra os arquivos das reads. Copie as 3 sequências aleatórias para um bloco de notas do windows. As sequencias devem ter o seguinte formato >seq_1 GCTACGATCGATCGATCGATCGATCGATCG >seq_2 TCGATCGATCGATCGATCAGCTAGCTAGCAT >seq_3 CGATCGATCGATCGATCGATCGATCGATCGA Agora colamos isso no ​
blastn​
(http://blast.ncbi.nlm.nih.gov/) e checamos em quais genes essas sequências são alinhadas. Note como as reads as vezes não se alinham completamente. Possivelmente isso deve surgir devido aos erros inerentes ao sequenciamento. Filtragem 8) Vamos tentar melhorar isso filtrando as sequências de entrada. Digite o comando ​
fastqc​
no terminal para checar como a qualidade do sequenciamento está: $ fastqc Agora abra o arquivo ​
fastq​
do seu diretório utilizando o programa fastqc (Nota: Mais informações do programa FastQC no link: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/​
) 9) No próximo passo vamos filtrar os arquivos fasta. Essa etapa é importante para a diminuição dos erros gerados durante o sequenciamento. Possuir reads mais puras gera menor tempo de execução do alinhamento, bem como o reaproveitamento de reads que seriam descartadas na etapa de alinhamento por possuir quantidade elevada de mismatchs contra o genoma. A ferramenta utilizada nessa etapa será o trimmomatic. O seguinte comando faz: *Remoção adaptadores *Remoção bases do início com baixa qualidade ou Ns *Remoção bases do fim com baixa qualidade ou Ns *Percorre o read com uma janela de 4, removendo quando a qualidade média por base é menor do que 15 *Descarta reads com cumprimento menor do que 20 bases Antes vamos criar um link simbólico do arquivo contendo os adaptadores. $ ln ­s /home/biomolcomp/data/TruSeq3­PE­2.fa TruSeq3­PE­2.fa $ java ­jar /home/biomolcomp/Trimmomatic­0.33/trimmomatic­0.33.jar PE ­threads 1 adrenal_1.fastq adrenal_2.fastq adrenal_1_paired.fastq.gz adrenal_1_unpaired.fastq.gz adrenal_2_paired.fastq.gz adrenal_2_unpaired.fastq.gz ILLUMINACLIP:TruSeq3­PE­2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:20 (Nota: mais informação do programa Trimmomatic em http://www.usadellab.org/cms/?page=trimmomatic​
) Repita o passo com os dados do cerebro Descompacte os arquivos que foram considerados filtrados e que mantiveram os pares após a filtragem. $ gzip ­d *_paired.fastq.gz ( ' *_paired.fastq.g ' faz com que todos os arquivos que iniciem com qualquer cadeia de caracteres e terminam com ​
_paired.fastq.gz​
sejam incluídos na ​
descompactação​
) Agora utilize novamente a ferramenta ​
fastqc​
para checar a qualidade dos dados. Mapeamento 10) Para realizar o mapeamento, primeiro temos que criar links simbólicos do genoma de referência e de nosso arquivo .gtf (​
lembre­se de estar no diretório aluno1_aluno2) $ ln ­s ../data/Ensembl_GRCh37/genome.* . $ ln ­s ../data/gencode19_lncrnas_custom.gtf . É possível olhar o conteúdo da pasta com ‘ls ­l’ 11) Agora vamos rodar TopHat2 utilizando os arquivos gerados até aqui. $ tophat ­p 1 ­G gencode19_lncrnas_custom.gtf ­o thout_adrenal genome adrenal_1_paired.fastq adrenal_2_paired.fastq Faremos o mesmo com a amostra de cérebro (Nota: informação formato .bam em http://genome.sph.umich.edu/wiki/SAM_Format​
) Montagem 12) A montagem dos transcritos pode ser feita com a ferramenta Cufflinks $ cufflinks ­p 4 ­o clout_adrenal thout_adrenal/accepted_hits.bam Repita o passo com a amostra de cérebro Expressão diferencial 13) Quais genes estão diferencialmente expressos? para responder a pergunta, primeiro vamos nos entrar na pasta aluno1_aluno2 e criar um arquivo de anotação de referência com o cuffmerge $ nano assemblies.txt nano é um editor de texto, e esse comando abre o editor criando o arquivo assembies.txt. Dentro desse arquivo vamos escrever: ./clout_adrenal/transcripts.gtf ./clout_brain/transcripts.gtf (para sair pressionar Ctrl + x) Salvar pressionando Y $ cuffmerge ­g gencode19_lncrnas_custom.gtf ­s genome.fa ­p 1 assemblies.txt Agora sim, vamos rodar o Cuffdiff $ cuffdiff ­o diff_out ­b genome.fa ­p 1 –L A,B ­u merged_asm/merged.gtf ./thout_adrenal/accepted_hits.bam ./thout_brain/accepted_hits.bam Explorando os resultados 14) Dentro da pasta diff_out encontramos vários resultados interessantes. Com​
ls podemos checar alguns desses arquivos gerados. $ ls Vamos manter o foco no arquivo gene_exp.diff . $ more gene_exp.diff Vamos agora selecionar apenas aqueles genes dados como diferencialmente expressos pelos testes estatísticos do cuffdiff. Se quisermos saber a quantidade de genes que caem nesse grupo, basta apenas digitar. $ grep 'yes$' gene_exp.diff |wc Vamos imprimir na tela agora apenas o nome dos genes dados como diferencialmente expressos. Incluiremos também a informação do fold­change relativo das duas amostras. $ grep 'yes$' gene_exp.diff |cut ­f3,10 (o comando após o pipe ' | ' faz com que apenas as colunas 3, 10 sejam incluídas) Nesse momento salvaremos 2 arquivos diferentes. Um com a informação do nome do gene. E um com a informação do nome do gene e do fold. $ grep 'yes$' gene_exp.diff |cut ­f3 > genes_names_signif.txt $ grep 'yes$' gene_exp.diff |cut ­f3,10 > gene_fold_signif.txt (o ' > ' redireciona toda a saída que seria mostrado na tela para um arquivo de nome escolhido. Nesse caso gene_fold_signif.txt) Podemos utilizar os comando feitos até aqui para checar este arquivo. $ less gene_fold_signif.txt Nesse momento já poderíamos definir genes interessantes que possuem nível de expressão diferentes nas diferentes amostras. Mas tentaremos enriquecer isso de uma forma mais completa. Para isso utilizaremos um comando do pacote homer. $ findMotifs.pl genes_names_signif.txt human homer_results/ ­start ­400 ­end 100 ­len 8,10 Os parametros ­start e ­end definem o intervalo em volta do inicio dos genes, onde motivos de ligação fatores de transcrição serão procurados. O parametro len define o tamanho dos motivos buscados 8,9,10. $ cd homer_results $ firefox homer_results/homerResults.html $ firefox homer_results/knownResults.html Referências Differential gene and transcript expression analysis of RNA­seq experiments with TopHat and Cufflinks. Trapnell C1​
, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, ​
Rinn JL, Pachter L. Trimmomatic: A flexible trimmer for Illumina Sequence Data Anthony M. Bolger, ​
​
Marc Lohse and Bjoern Usadel Simple Combinations of Lineage­Determining Transcription Factors Prime cis­Regulatory Elements Required for Macrophage and B Cell Identities Heinz S, Benner C, Spann N, Bertolino E et al