Analisando miRNAs com mirDeep2

Olá pessoal! Neste tutorial mostrarei como podemos identificar miRNAs a partir de dados de sequenciamento de nova geração usando a ferramenta mirDeep2.

O que são miRNAs?

MicroRNAs (miRNAs) são moléculas de RNA não-codificante (ncRNAs) com tamanho geralmente entre 21 e 25 pares de base (bps) que atuam no processo de regulação pós-transcricional da expressão gênica de organismos eucariotos. Os miRNAs interagem com o RNA mensageiro (mRNA) e reduzem a sua expressão através de diferentes mecanismos, como indução da clivagem (feita pela enzima Dicer), e redução da eficiência da tradução.

Grande parte dos miRNAs conhecidos são derivados de regiões gênicas do genomas, apesar de nem sempre estarem na porção codificante (coding DNA sequence, ou apenas CDS) composta pelos exons. Além disso, nem sempre o gene alvo do miRNAs é o mesmo do qual este é derivado, e em alguns casos a sequência final do miRNA não é um reflexo direito da porção genômica que o originou, visto que estas moléculas podem passar por diferentes modificações enzimáticas durante a sua biogênese. Nestes casos, essa sub-classe de miRNAs é denominada “isomiR”.

De modo geral, duas vias principais podem dar origem aos miRNAs. A primeira envolve o processamento feito pela enzima Drosha, uma ribonuclease que processa transcritos primários de miRNAs denominados pri-miRNAs, cuja clivagem origina os precursores de miRNAs (pre-miRNAs). Alternativamente, o pre-miRNAs podem ser derivados também do processo de splicing de alguns mRNAs realizado pelo spliceossoma. Neste segundo caso, os miRNAs originados a partir de introns são denominados “mirtrons”.

Biogênese de miRNAs

A mesma abordagem utilizada para o sequenciamento de mRNA, baseada no uso da transcriptase reversa, pode ser aplicada também para miRNAs, com apenas alguns ajustes no protocolo, como o uso adaptadores específicos, etapas de circularização (dependendo do protocolo) e seleção por tamanho. A partir dos miRNAs identificados nessas análises é possível se determinar genes-alvo, vias de regulação e variações na expressão destes RNAs de acordo com o contexto do organismo em análise.

Neste tutorial, faremos a identificação de miRNAs de modo a identificar moléculas já conhecidas e conservadas em outros organismos, bem como novos miRNAs ainda não descritos.

Estrutura do projeto

Nosso projeto deverá conter a seguinte estrutura:

project/
    alignments/
    bin/
    mirbase/
    reads/
    scripts/
        fastq_to_fasta.py
        index_reference.sh
        map_reads.sh
        annotate_mirnas.sh
    reference/
    reference_indexed/
    Makefile
    samples

Da mesma forma que nos tutoriais anteriores, neste utilizaremos um arquivo Makefile para centralizar e organizar os scripts que devem ser executados em cada etapa. Neste projeto, o arquivo terá a seguinte estrutura:

index_reference:
    bash scripts/index_reference.sh
convert_fastq_to_fasta:
    bash scripts/convert_fastq_to_fasta.sh
map_reads:
    bash scripts/map_reads.sh
annotate_mirnas:
    bash scripts/annotate_mirnas.sh

Instalando dependências

Primeiro vamos instalar alguns programas para fazer o download dos datasets e auxiliar nas etapas de conversão de arquivos. Os datasets com dados de sequenciamento serão derivados do NCBI Short Read Archive (SRA) e para isso precisamos instalar o pacote sra-toolkit. Já o processamento dos arquivos será realizado por um script em Python que usa a biblioteca BioPython. Para instalar estes pacotes, execute os seguintes comandos no terminal:

$ sudo apt-get install sra-toolkit python3-pip
$ sudo pip3 install biopython

Agora, vamos instalar as dependências para rodar o mirDeep2. Para isso, entre na pasta bin/ e execute os seguintes comandos:

$ sudo apt-get install bowtie
$ wget -O http://eddylab.org/software/squid/squid.tar.gz
$ tar -xvf squid.tar.gz
$ cd squid-1.9g
$ bash configure
$ make
$ sudo make install
$ cd ..
$ git clone https://github.com/eb00/randfold
$ cd randfold/src
$ bash configure
$ make
$ sudo cp randfold /usr/local/bin/
$ cd ..
$ git clone https://github.com/rajewsky-lab/mirdeep2
$ cd mirdeep2
$ perl install.pl
$ cd ../..

Obtendo o dataset

Neste tutorial vamos analisar dados de sequenciamento de miRNAs de uma amostra de soro equino. Os dados de sequenciamento serão obtidos do SRA, e para fazer o download vamos utilizar a ferramenta fastq-dump, distribuida com o pacote sra-toolkit que instalamos. Para baixar os dados, edite o arquivo arquivo samples e adicione o seguinte código de acesso:

ERR1699030

Agora, pela linha de comando, entre na pasta reads/ e execute o seguinte comando:

$ fastq-dump ERR1699030

Ao executar este comando, o um arquivo ERR1699030.fastq será gerado. Note que neste caso estamos realizando o download dos dados manualmente, mas caso tenhamos várias amostras para analisar, podemos executar esta etapa com um script simples, download_data.sh, como o mostrado abaixo.

cd reads/
while read sample; do
    fastq-dump $sample
done < ../samples

Agora, precisamos baixar a sequência do genoma do Equus caballus em formato FASTA, disponível no Ensembl (link). Salve o arquivo na pasta reference/ , descompacte e renomeie o arquivo FASTA para reference_raw.fa. Agora, pelo terminal, entre na pasta reference/ e execute o seguinte comando:

$ cat raw_reference.fa | cut -d' ' -f1 > reference.fa

Este comando servirá para remover do header de cada sequência dos genomas as descrições, mantendo apenas o ID. Na prática, isso evita algumas incompatibilidades na forma como o Bowtie e o mirDeep2 lidam com a estrutura do arquivo FASTA.

Por fim, precisamos baixar alguns dados de miRNAs já conhecidos, tanto do nosso organismo de interesse (Equus caballus) quanto de um organismo próximo (que neste caso será o Homo sapiens). Para isso, vamos utilizar alguns dados derivados do mirBase, que podem ser obtidos nos links abaixo:

• E. caballus miRNAs;
• H. sapiens miRNAs;

Agora, descompacte ambos arquivos na pasta mirbase.

Pre-processamento das leituras

Agora precisamos converter as leituras do formato FASTQ obtidas do SRA para o formato FASTA aceito pelo mirDeep2. Para isso, é preciso garantir que não tenhamos sequências repetidas (redundância), e ajustar o header para um formato específico do programa que segue a estrutura:

>{id da amostra}_{id da sequencia}_x{ocorrência da sequência}
{sequencia}

Exemplo:

>001_1_x30
GCATTTGTGGTTCAGTGGTGGAATTCTCGC
>001_2_x1
TTTAGTTGGTTTTCGGAACTGAGC
>001_3_x2
GAGCAGTCAGGATGGCCGAGCGGTC
>001_4_x1
AAGCTCGGTCTGAGGCCAAACAG
>001_5_x1
AAGAGTGGAGAGGATCTTTGTC

Para fazer isso, vamos utilizar ofastq_to_fasta.py, executado no Makefile pela regra fastq_to_fasta. Este script deve conter o seguinte código:

from Bio import SeqIO
def clear_sequence(_sequence):
    sequence = ''
    _sequence = _sequence.upper()
    for base in _sequence:
        if base in 'ATCGN':
            sequence += base
        else:
            sequence += 'N'
    return sequence
for s, sample in enumerate(open('samples')):
    sample = sample.strip('\n')
    records = SeqIO.parse(
        open("reads/%s.fastq"%sample),
        'fastq')
    
    sequence_uniq = set([])
    sequence_count = {}
    for record in records:
        sequence = str(record.seq)
        sequence_uniq.add(str(record.seq))
        if sequence not in sequence_count:
            sequence_count[sequence] = 0
        sequence_count[sequence] += 1
    
    with open('reads/%s.fasta'%sample, 'w') as writer:
        for seq, sequence in enumerate(sequence_uniq):
            if len(sequence) < 20:
                continue
            sequence = clear_sequence(sequence)
            writer.write(">%s_%s_x%s\n%s\n"%(
                str(s+1).zfill(3),
                seq+1, 
                sequence_count[sequence],
                sequence))

Na pasta do projeto, execute o seguinte comando:

$ make fastq_to_fasta

Alinhamento das leituras

Da mesma forma que na análise de RNA-Seq diferencial, antes de realizar o alinhamento da leitura precisaremos indexar a sequência de referência. Para isso, vamos utilizar o seguinte código no nosso script index_reference.sh, localizado na pasta scripts/.

mkdir -p reference_indexed
cd reference
bowtie-build --large-index  \
    --bmax 16777216 \
    --dcv 256 \
    --threads 4 \
    reference.fa \
    ../reference_indexed/reference

obs: Os argumentos --bmax e --dcv controlam alguns parâmetros que impactam diretamente o uso de memória durante o processo de indexação, e foram otimizados para funcionar na minha máquina (com 16 Gb de RAM). Dependendo do caso, pode ser necessário ajustar estes valores, abaixando o valor de --bmax para evitar erros do tipo segmentation default (core dump).

O script map_reads.sh, executará o processo de alinhamento propriamente dito. Este alinhamento é realizado por um script do próprio mirDeep2, mapper.pl, que utilizar o Bowtie como base e converte o seu resultado para um formato mais adequado para as análises posteriores. A automatização deste processo é feita pelo script map_reads.sh, executado pela regra map_reads do Makefile.

while read sample; do
    mkdir -p alignments/$sample
    cd alignments/$sample
    rm $sample.aligned.aff
    perl ../../bin/mirdeep2/bin/mapper.pl 
         ../../reads/$sample.fasta \
         -p ../../reference_indexed/reference \
         -c -i -j -l 18 -m \
         -t $sample.aligned.aff
    cd ../../
done < samples

Dentro da pasta do projeto, execute os seguintes comandos no terminal:

$ make index_reference
$ make map_reads

Identificação de miRNAs

Agora, vamos realizar o processo de identificação do miRNAs a partir do alinhamento. Para isso, utilizamos o script principal do mirDeep2, mirDeep2.pl, que recebe o resultado do alinhamento em formato .aff gerado script mapper.pl e os mirRNAs conhecidos em formato FASTA. Este processo será realizado pelo script annotate_mirnas.sh , executado pela regra annotate_mirnas do Makefile.

while read sample; do
    cd alignments/$sample
    perl ../../bin/mirdeep2/bin/miRDeep2.pl \
         ../../reads/$sample.fasta \
         ../../reference/reference.fa \
         $sample.aligned.aff \
         ../../mirbase/mirbase_eca.fasta \
         ../../mirbase/mirbase_hsa.fasta \
         none
    cd ../../
done < samples

Execute este script com o comando:

$ make annotate_mirnas

Obs: Este etapa pode demorar algumas horas para ser finalizada.

Ao final, na pasta alignments/ teremos uma subpasta com o nome de cada amostra (neste caso, ERR1699030), que conterá os arquivos com o relatório da anotação em formato BED, HTML e CSV. Um exemplo do relatório em HTML está apresentado abaixo.

Exemplo de resultado da identificação de miRNAs gerada pelo mirDeep2.

Agora que temos os miRNAs identificados, é possível identificar genes alvo com uso de ferramentas TargetScan, bem como realizar a quantificação com ferramentas do próprio mirDeep2 para analisar a expressão diferencial estas moléculas em diferentes condições, de modo similar com a análise que realizados no tutorial de RNA-Seq.