Quantificação de reads de dados de RNA Seq - lmigueel/Bioinformatica GitHub Wiki

1. Sobre

Kalisto é um programa para quantificar abundâncias dos transcritos de dados de RNA-Seq, ou mais geralmente de reads de NGS. Ele é baseado na nova ideia de pseudoalinhamento para determinar rapidamente a compatibilidade de reads, sem a necessidade de alinhamento. Em benchmarks com dados RNA-Seq padrão, kallisto pode quantificar 30 milhões de leituras humanas em menos de 3 minutos em um computador desktop Mac usando apenas as sequências de leitura e um índice de transcriptoma que leva menos de 10 minutos para ser construído. O pseudoalinhamento de reads preserva as principais informações necessárias para a quantificação e, portanto, o kallisto não é apenas rápido, mas também tão preciso quanto as ferramentas de quantificação existentes.

O manual do kallisto pode ser acessado AQUI.

2. Instalação

Via conda:

conda install kallisto

Para o download interno acesse: https://pachterlab.github.io/kallisto/download

3. Uso sobre os dados de RNA-Seq

As principais etapas do kallisto são duas:

  1. CRIAR UM ÍNDICE DO FASTA (kallisto index)
  2. QUANTIFICAR (kallisto quant)

Vamos começar pelo index. Para fins didáticos, vamos supor que você possui um arquivo chamado _all_transcripts.fasta e uma análise por tratado e controle, ambos paired-end, gerando 12 arquivos .fastq. Essa ideia é a mesma na seção Análise de RNA-Seq com genoma de referência.

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

Gerando o index chamado index_transcritos:

kallisto index -i index_transcritos all_transcripts.fasta

Agora vamos quantificar cada biblioteca .fastq. Importante lembrar de executar bootstraps (-b), alterar o número de processadores (-t) e adicionar o nome da pasta output (-o).

kallisto quant -i index_transcritos -o controle_rep1 -b 1000 -t 7 controle_rep1_1.fastq controle_rep1_2.fastq
kallisto quant -i index_transcritos -o controle_rep2 -b 1000 -t 7 controle_rep2_1.fastq controle_rep2_2.fastq
kallisto quant -i index_transcritos -o controle_rep3 -b 1000 -t 7 controle_rep3_1.fastq controle_rep3_2.fastq
kallisto quant -i index_transcritos -o tratado_rep1 -b 1000 -t 7 tratado_rep1_1.fastq tratado_rep1_2.fastq
kallisto quant -i index_transcritos -o tratado_rep2 -b 1000 -t 7 tratado_rep2_1.fastq tratado_rep2_2.fastq
kallisto quant -i index_transcritos -o tratado_rep3 -b 1000 -t 7 tratado_rep3_1.fastq tratado_rep3_2.fastq

IMPORTANTE: Caso você tenha arquivos single-end você deve utilizar as flags --single -l X -s 0.001 arquivo.fastq, onde o valor de X é o tamanho do seu fragmento (pode ser o do read). O desvio padrão pode ser calculado, mas se o tamanho pela anállise de qualidade está consistente, coloque 0.001 para precauções.

Note que criaremos as pastas do controle e tratado, com seus respectivos nomes. Dentro de cada pasta existem três essenciais arquivos:

a) abundance.h5: arquivo contendo informações dos bootstraps e é utilizado pelo Sleuth, por exemplo. b) abundance.tsv: possui as quantificações de cada gene do arquivo all_transcripts.fasta c) run_info.json: possui as informações do alinhamento

4. Montando a tabela final

A tabela final pode ser realizada por você em um script in-house, mas ela pode ser gerada através do script abundance_estimates_to_matrix.pl que está no pacote TRINITY dentro da pasta /util/. Acesse aqui: script. Acesse o pacote Trinity aqui.

Executar este script é bem simples, mas precisamos ter instalado o edgeR, que deverá ser executado no R:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("edgeR")

Caso não tenha o R instalado, execute:

 conda install -c r r

ou execute o download do R no Linux AQUI.

Para executar o script abundance_estimates_to_matrix.pl precisamos que todos os arquivos .tsv tenham seu nome modificado de forma eficiente. Para isso, dentro da pasta onde estão todas as pastas quantificadas, eu costumo copiar os arquivos .tsv com seus nomes correspondentes. Como, por exemplo:

cp controle_rep1/abundance.tsv controle_rep1.tsv
cp controle_rep2/abundance.tsv controle_rep2.tsv
cp controle_rep3/abundance.tsv controle_rep3.tsv
cp tratado_rep1/abundance.tsv tratado_rep1.tsv
cp tratado_rep2/abundance.tsv tratado_rep2.tsv
cp tratado_rep3/abundance.tsv tratado_rep3.tsv

Agora sim. Um detalhe importante é o arquivo gene-transcrito. Caso você tenha, o kallisto criará dois arquivos: uma tabela para os transcritos e outra agrupada para os genes. O bom que este script já libera a tabela de counts, TPM e TMM (requer edgeR instalado). Como não temos este arquivo tabular, vamos adicionar o argumento 'none'.

perl ~/trinity/util/abundance_estimates_to_matrix.pl --est_method kallisto --gene_trans_map none controle_rep1.tsv controle_rep2.tsv controle_rep3.tsv tratado_rep1.tsv tratado_rep2.tsv tratado_rep3.tsv

No final teremos os arquivos:

*isoform.counts.matrix: matriz de counts (pode ser utilizada no [DESeq2](https://github.com/lmigueel/Bioinformatica/wiki/DESeq2)

*isoform.TPM.not_cross_norm: matriz de TPM 

*isoform.TMM.EXPR.matrix: matriz normalizada por TMM (pode ser utilizada no em [PCA](https://github.com/lmigueel/Bioinformatica/wiki/PCA--(Principal-Component-Analysis)-em-dados-de-RNA-Seq)

*TMM_info.txt: estatística dos factores de normalização pela profundidade do sequenciamento.

*runTMM.R: script em R para realizar a normalização em TMM

A partir deste ponto, pode ser executadas análises de expressão diferencial.

Citação

Nicolas L Bray, Harold Pimentel, Páll Melsted and Lior Pachter, Near-optimal probabilistic RNA-seq quantification, Nature Biotechnology 34, 525–527 (2016), doi:10.1038/nbt.3519

Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011 May 15;29(7):644-52. doi: 10.1038/nbt.1883. PubMed PMID: 21572440.