Análise de RNA Seq com genoma de referência - lmigueel/Bioinformatica GitHub Wiki

1. Sobre

O sequenciamento de alto rendimento de mRNA (RNA-Seq) tornou-se o método padrão para medir e comparar os níveis de expressão gênica em uma ampla variedade de espécies e condições. Os experimentos de RNA-Seq geram conjuntos de dados muito grandes e complexos que exigem software rápido, preciso e flexível para reduzir os dados lidos brutos a resultados compreensíveis. HISAT2 e StringTie são ferramentas para análises de experimentos de RNA-Seq. Juntos, eles permitem que os cientistas alinhem os reads a um genoma e reconstroem os transcritos, incluindo novas variantes (sequências não contidas na referência) e calculem a abundância desses transcritos em cada amostra. As análises de expressão diferencial podem ser realizadas pelo Sleuth ou DESeq2, além de outras ferramentas. Acesse os manuais desta WIKI para as análises posteriores.

2. Instalação

A maneira mais simples de instalar o HISAT2 e StringTie é através do conda:

 conda install -c bioconda hisat2

 conda install -c bioconda stringtie

O manual de cada um destes softwares podem ser encontrados aqui: HISAT2 e StringTie.

Você pode acessar o help dos programas via conda através da flag --help.

Os códigos fonte para uso interno, caso o conda não esteja funcionando, podem ser encontrados aqui: StringTie e HISAT2.

Parte 1. HISAT2

Primeira parte para a reconstrução dos transcritos a partir de um genoma de referência é criar um índice eficiente do genoma. Assim, as buscas dos reads no MAPEAMENTO fica mais precisa e eficiente.

Você precisa de dois arquivos:

ARQUIVO FASTA DO GENOMA
GFF DO GENOMA

Para fins apenas didáticos, supomos que o nome dos arquivos são genome.fasta e genome.gff.

Primeiro passo é utilizar os scripts Python para extrair informações dos sítios de splicing e exons do arquivo de anotação de genes (GFF). Estes arquivos já serão instalados pelo conda, e estão na pasta do source code (sem o prefixo hisat2_).

Depois utilizamos o hisat2-build para criar o índice chamado reference_genome.

hisat2_extract_splice_sites.py genome.gff > genome.ss
hisat2_extract_exons.py genome.gff > genome.exon
hisat2-build --ss genome.ss --exon genome.exon genome.fasta reference_genome

Agora podemos executar o mapeamento. Para fins didáticos nós teremos duas condições (tratado e controle) com três réplicas biológicas paired-end cada uma. No fim, teremos doze arquivos:

controle_rep1_1.fastq
controle_rep1_2.fastq

controle_rep2_1.fastq
controle_rep2_2.fastq

controle_rep3_1.fastq
controle_rep3_2.fastq

tratado_rep1_1.fastq
tratado_rep1_2.fastq

tratado_rep2_1.fastq
tratado_rep2_2.fastq

tratado_rep3_1.fastq
tratado_rep3_2.fastq

Executamos o mapeamento, com 5 processadores (-p 5) e com a saída de um arquivo .sam (-S). A opção --dta é super importante, pois ela gera um relatório feito sob medida para os programas de reconstrução, incluindo StringTie. Com esta opção, o HISAT2 requer comprimentos de âncora mais longos para a descoberta de novos transcritos. Isso leva a menos alinhamentos com short-anchors, o que ajuda os montadores a melhorar significativamente o uso de computação e memória. Caso você tenha arquivos single-end use a flag -U.

CUIDADO, antes de rodar todos de uma vez verifique o total de processadores disponíveis no servidor.

hisat2 --dta -p 5 -x reference_genome -1 controle_rep1_1.fastq -2 controle_rep1_2.fastq -S controle_rep1.sam
hisat2 --dta -p 5 -x reference_genome -1 controle_rep2_1.fastq -2 controle_rep2_2.fastq -S controle_rep2.sam
hisat2 --dta -p 5 -x reference_genome -1 controle_rep3_1.fastq -2 controle_rep3_2.fastq -S controle_rep3.sam
hisat2 --dta -p 5 -x reference_genome -1 tratado_rep1_1.fastq -2 tratado_rep1_2.fastq -S tratado_rep1.sam
hisat2 --dta -p 5 -x reference_genome -1 tratado_rep2_1.fastq -2 tratado_rep2_2.fastq -S tratado_rep2.sam
hisat2 --dta -p 5 -x reference_genome -1 tratado_rep3_1.fastq -2 tratado_rep3_2.fastq -S tratado_rep3.sam

Antes de executarmos o StringTie, necessitamos converter os arquivos para .bam. Caso deseje verificar algumas métricas essenciais para análises de .bam, acesse AQUI. Utilizaremos o samtools.

samtools sort -o controle_rep1.bam controle_rep1.sam
samtools sort -o controle_rep2.bam controle_rep2.sam
samtools sort -o controle_rep3.bam controle_rep3.sam
samtools sort -o tratado_rep1.bam tratado_rep1.sam
samtools sort -o tratado_rep2.bam tratado_rep2.sam
samtools sort -o tratado_rep3.bam tratado_rep3.sam

Com os arquivos .bam em mãos, vou executar a reconstrução dos transcritos.

Parte 2 - StringTie

Antes de começar necessitamos criar um arquivo contendo o nome de todos os arquivos que serão utilizados na função merge, que junta todas as coordenadas dos genes encontrados em uma única referência através dos overlaps observados.

echo -e 'controle_rep1_stringtie.gtf\ncontrole_rep2_stringtie.gtf\ncontrole_rep3_stringtie.gtf\ntratado_rep1_stringtie.gtf\ntratado_rep2_stringtie.gtf\ntratado_rep3_stringtie.gtf ' > mergelist.txt

O arquivo mergelist.txt está da forma:

controle_rep1_stringtie.gtf
controle_rep2_stringtie.gtf
controle_rep3_stringtie.gtf
tratado_rep1_stringtie.gtf
tratado_rep2_stringtie.gtf
tratado_rep3_stringtie.gtf

Agora podemos partir para a análise do StringTie e gerar os arquivos de coordenadas (GTF) de cada biblioteca. O nome do arquivo deve ser o mesmo no arquivo mergelist.txt.

stringtie controle_rep1.bam -p 5 -G genome.gtf -o controle_rep1_stringtie.gtf 
stringtie controle_rep2.bam -p 5 -G genome.gtf -o controle_rep2_stringtie.gtf 
stringtie controle_rep3.bam -p 5 -G genome.gtf -o controle_rep3_stringtie.gtf 
stringtie tratado_rep1.bam -p 5 -G genome.gtf -o tratado_rep1_stringtie.gtf 
stringtie tratado_rep2.bam -p 5 -G genome.gtf -o tratado_rep2_stringtie.gtf 
stringtie tratado_rep3.bam -p 5 -G genome.gtf -o tratado_rep3_stringtie.gtf

Agora podemos executar o merge dos arquivos .GTF e criar uma única referência.

stringtie --merge -G genome.gtf -o stringtie_merged.gtf mergelist.txt

Agora podemos extrair as sequências correspondentes.

Parte 3 - Extraindo as sequências

Agora que temos o arquivo final de coordenadas genômicas (stringtie_merged.gtf), podemos extrair as sequências fasta do arquivo através do programa GFFREAD. Vamos a instalação.

Conda:

 conda install -c bioconda gffread

Para extrair as sequências (arquivo all_transcripts.fasta), basta realizar:

gffread stringtie_merged.gff -g genome.fasta -w all_transcripts.fasta

PRONTO!

Já possuo o transcriptoma de referência

Caso você já possua o transcriptoma de referência (caso for um arquivo de cDNA), basta partir para as análises de quantificações (ACESSE AQUI) e de expressão diferencial (ACESSE AQUI).

Citação

Pertea, Mihaela, et al. "Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown." Nature protocols 11.9 (2016): 1650.