Análise de expressão gênica diferencial com RNA-Seq (reference-guided)
Olá pessoal! Neste tutorial mostrarei como podemos realizar a análise de expressão diferencial a partir de dados de RNA-Seq. Este tipo de análise tem se tornada extremamente popular com o barateamento do sequenciamento de DNA promovido pelas plataformas de nova geração (NGS), mas ainda temos pouco material disponível em português que aborde o seu passo a passo.
O que é RNA-Seq
O sequenciamento de RNA não é uma coisa nova. Na verdade, antes mesmo de existirem técnicas para análise de sequenciamento de DNA, como os métodos descritos por Sanger e Maxam & Gilbert (ambos de 1977), sequências de RNAs não-codificantes (como tRNAs) há haviam sido publicadas, sendo obtidas por degradação química ou enzimática.
Entre os anos 80 e começo dos anos 2000, diferentes abordagem foram desenvolvidas para a análise do transcriptoma (conjunto de transcritos de um organismo em uma condição) através do sequenciamento de Sanger, como EST (expressed sequence tags) e SAGE (serial analysis of gene expression). Esses métodos, no entanto, era limitados pelo baixo throughput das plataformas de Sanger, que, por “rodada”, permitiam no máximo que 96 amostras fossem sequências por vez. Por conta disso, neste período os microarranjos (DNA microarray) se tornaram uma alternativa útil para o estudo em larga escala da expressão gênica, mas estes apresentavam uma outra limitação: apenas regiões conhecidas poderiam ser analisadas.
Com o advento do sequenciamento de nova geração (NGS), o throughput gerado por cada rodada de sequenciamento passou de algumas centenas de de kilobases (Kbs) para vários gigabases (Gbs). Apesar de ter sido inicialmente direcionado para o estudo de sequências genômicas, o NGS foi rapidamente adotado também por muitas outras “ciências ômicas”, dentre elas a transcriptômica. Para fins de diferenciação, adotou-se o termo RNA-Seq para se referir à análise do transcriptoma através destas novas tecnologias, apesar dos primeiros trabalhos produzidos com esta abordagem ainda terem usado o nome EST.
A análise do transcriptoma com RNA-Seq pode ser realizada com ou sem uma sequência genômica de referência. No primeiro caso chamamos esta análise de reference-guided, e esta abordagem geralmente é empregada quando estamos analisando genomas microbianos ou eucarióticos “modelo” (ex: Homo sapiens, Mus musculus). Já a análise sem referência (de novo) é geralmente realizada para organismos eucarióticos não-modelo.
No caso da análise reference-guided, os resultados do sequenciamento, as leituras (reads), são alinhadas contra o genoma, e a contagem de leituras mapeadas em cada gene é usada para pedir a sua expressão. Já no caso de novo, é necessário se fazer uma montagem dos transcritos, para que depois seja possível realizar a sua quantificação.
Como a expressão gênica varia de acordo com as condições na qual o organismo se encontra, é possível se utilizar a análise de RNA-Seq para mensurar o efeito de diferentes tratamentos na expressão de diferentes genes simultaneamente. Este tipo de análise, denominada expressão gênica diferencial, permite se entender como o perfil de expressão de um determinado organismo é alterado a ser submetido a uma determinada condição. Entretanto, diferente de métodos clássicos de biologia molecular, como os baseados em PCR, neste caso não precisamos escolher um gene normalizador.
Após identificarmos genes com expressão aumentada (up-regulated) ou diminuída (down-regulated), podemos mapear processos biológicos que estão sendo modulados em resposta às mudanças na qual o organismo de encontra. Este tipo de abordagem é extremamente útil para se entender diferentes processos biológicos, como a dinâmica de expressão gênica em células cancerosas submetidas a certos medicamentos, ou em plantas submetidas à uma determinada condição de estresse, por exemplo.
O que faremos?
De modo a entender as etapas da análise e dados de RNA-Seq, vamos realizar uma análise reference-guided da expressão gênica de Leptospira interrogans e comparar os perfis de expressão desta bactéria in vitro e in vivo (durante a infecção). Os dados usados neste projeto são derivados deste artigo (Caimano et al, 2014), e são públicos :D
Instalando dependências
Para este projetos, precisamos de um ambiente Ubuntu Linux (≥16.04), com R ≥ 3.6. Para que possamos realizar as análises, vamos instalar os seguintes programas e bibliotecas:
$ apt install rna-star samtools python-htseq \
libxml2-dev python3-pip libssl-dev \
sra-toolkit
$ pip3 install pandasAlém disso, também precisamos instalar alguns pacotes do R. Para isso, crie um script install_packages.R na pasta scripts/ e execute com o comando Rscript scripts/install_packages.R. Este script deve conter o seguinte código:
#!/usr/bin/env Rscriptif (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager",
repos='http://cran.us.r-project.org')
BiocManager::install(c("DESeq","DESeq2","edgeR", "HTSFilter", "limma"), dependencies=TRUE)
Estrutura do projeto
Neste projeto, usaremos a seguinte estrutura de arquivos e diretórios:
project/
alignments/
differential_expression/
read_count/
reads/
reference/
reference.fasta
reference.gff3
reference_indexed/
scripts/
install_packages.R
download_data.sh
index_reference.sh
count_reads_on_genes.sh
merge_count_tables.py
differential_expression.R
Makefile
samples
samples_groupsNote que com exceção das pastas scripts/ e reference/, todas as demais estão vazias. Isso ocorre pois todos os demais arquivos serão obtidos via scripts, desde o download dos dados de sequenciamento até a análise de expressão diferencial :D
Os arquivos do genoma de referência (reference.fasta e reference.gff3) estão disponíveis neste link, basta baixar de descompactar na pasta reference/.
Por padrão eu costumo organizar as etapas na ordem que elas devem ser executadas utilizando arquivos do tipo Makefile. Neste caso, o conteúdo do arquivo fica assim:
download_data:
bash scripts/download_data.shindex_reference:
bash scripts/index_reference.shmap_reads:
bash scripts/map_reads.sh
bash scripts/index_sort_reads.shcount_reads:
bash scripts/count_reads_on_genes.sh
python3 scripts/merge_count_tables.pydifferential_expression_analysis:
mkdir -p differential_expression; \
cd differential_expression; \
Rscript ../scripts/differential_expression.R
Obtendo o dataset
Os dataset serão obtidos a partir do banco de dados SRA (short read archive) do NCBI, e para isso precisamos utilizar a ferramenta fastq-dump, que é parte do pacote sra-toolkit. Estes datasets representa 2 condições (“tratamentos”) diferentes às quais a bactéria foi submetida (in vitro e in vivo), cada um com 3 replicatas, sendo cada um indexado SRA com um código de acesso diferentes. Estes códigos de acesso devem ficar no arquivo samples, sendo eles:
SRR1071264
SRR1071263
SRR1071262
SRR1071261
SRR1071260
SRR1071259Além disso, precisamos também ter um arquivo contendo as informações de quais amostras correspondem a cada tratamento. Nestas informações devem ficar no arquivo samples_groups:
sample,group
SRR1071264,CONTROL
SRR1071263,CONTROL
SRR1071262,CONTROL
SRR1071261,TREATMENT
SRR1071260,TREATMENT
SRR1071259,TREATMENTPara utilizar o fastq-dump precisamos informar como argumento um código de acesso de SRA. Ex:
$ fastq-dump SRR1071264Este comando fará o download da amostra de código de SRR1071264 e gerará um arquivo com o nome SRR1071264.fastq. Caso queiramos baixar cada um dos datasets disponíveis no arquivo samples podemos repetir este mesmo comando para cada uma das amostras …. ou fazer um script para repeti-lo para cada código de acesso (e é exatamente isso que faremos 😜). O script download_data.sh , localizado na pasta scripts/, deve conter o seguinte código:
mkdir -p reads/
cd reads/while read sample; do
fastq-dump $sample
done < ../samples
Ao final, a pasta reads/ conterá um arquivo .fastq para cada código de aceso no arquivo samples. Agora, basta digitar o comando make download_data no diretório raiz para baixar os dados.
Alinhamento das leituras
Apesar precisamos alinhar as sequências contra o genoma de referência, mas antes precisamos indexar o mesmo. Este processo é necessário pois facilita a execução do algoritmo de busca. Diferente de programas de alinhamento clássicos, como o BLAST, os algoritmos utilizados para alinhamento de sequências de leitura de NGS geralmente utilizam estruturas de dados compactadas, como as produzidas a partir de método de Burrows–Wheeler (BWT). Este processo torna a busca mais rápida, mas por outro lado exige um pré-processamento (indexação) e reduz a “sensibilidade” do alinhamento.
Neste tutorial utilizarmos a ferramenta STAR, um programa de alinhamento de leituras específico para RNA-Seq. Diferente de programas de mapeamento genéricos, como o Bowtie e o BWA, o STAR utiliza um algoritmo splice-aware, que considera as funções exon-intron durante o processo de alinhamento. Alternativamente, é possível também utilizar o TopHat para esta finalidade.
Para realizarmos a indexação, podemos o seguinte comando no diretório do projeto:
mkdir -p reference_indexed
STAR \
--runMode genomeGenerate \
--genomeFastaFiles reference/reference.fasta \
--genomeDir reference_indexedEste comando também pode ser embutido em um script, index_reference.sh , na pasta scripts/ , para ser executado com o comando make index_reference. A principal vantagem de manter os comando em scripts e organiza-los com o Makefile é garantir a padronização e replicação do protocolo de análise.
Agora que temos a sequência indexada, podemos fazer o alinhamento. Este processo deve ser realizado para cada arquivo de leituras, e para não precisarmos repetir este passo vários vezes vamos criar outro script, map_read.sh. Este script deve conter o seguinte código:
while read sample; do
mkdir -p alignments/$sample
STAR \
--runThreadN 4 \
--genomeDir reference_indexed \
--readFilesIn reads/$sample.fastq \
--outFileNamePrefix alignments/$sample/outdone < samples
Após o alinhamento, precisamos realizar uma etapa de ordenamento de indexação do mesmo. O processo de ordenamento, realizado pelo script index_sort_reads.sh, é importante pois o possibilita facilita o processo de contagem das leituras mapeadas em cada região, pois caso contrário, o algoritmo teria de percorrer praticamente todo o arquivo de alinhamento para encontrar leituras mapeadas sobre cada gene.
cd alignments
while read sample; do
cd $sample
samtools view -Sb -F4 outAligned.out.sam > unsorted.bam
samtools sort -o sorted.bam unsorted.bam
samtools index sorted.bam
cd ..
done < ../samplesContando as leituras alinhadas em cada gene
Agora que as sequências estão alinhadas, vamos realizar a contagem de quantas leituras mapearam em cada gene. Para isso, utilizaremos a ferramenta HTSeq. Este programa recebe como input o arquivo do alinhamento e uma anotação do genoma de interesse no forma GFF3. Além disso, precisamos informar qual tipo de feature (ex: CDS, gene, exon, mRNA) queremos analisar, bem como o atributo que será usado como key. Da mesma forma que no alinhamento, precisamos repetir esta etapa para cada amostra, e para isso criamos o script count_reads_on_genes.sh.
FEATURE_TYPE=CDS
GENE_KEY=locus_tag
cd read_count
while read sample; do
htseq-count \
-t $FEATURE_TYPE \
-i $GENE_KEY \
-f bam \
../alignments/$sample/sorted.bam \
../reference/reference.gff3 \
> $sample.raw.count
cat $sample.raw.count > $sample.count
rm $sample.raw.count
done < ../samplesAgora, precisamos concatenar as contagens para as diferentes amostras em um mesmo arquivo, de modo que possamos usar como entrar para a análise de expressão. Para isso, usamos o seguinte scripts em Python (existem várias formas de fazer isso, inclusive com o próprio R, mas eu, particularmente preciso manipular dados tabulares usando o Pandas).
import pandas as pdsamples = [l.strip("\n") for l in open("samples").readlines()]
samples_df = [pd.read_csv('read_count/%s.count'%sample, sep="\t", header=None, names=['locus_tag', sample]) for sample in samples]
merged_df = samples_df[0][['locus_tag']]for s,sample in enumerate(samples):
merged_df[sample] = samples_df[s][sample]merged_df.drop(merged_df.tail(5).index,inplace=True)
merged_df.to_csv("read_count/merged.csv", index=False)
A execução destes dois scripts é realizada pela regra count_reads do Makefile. Ao fim, teremos um arquivo merged.csv na pasta read_count/ , além dos arquivos específicos de cada amostra. Agora, só falta a análise da expressão diferencial :D
Análise da expressão diferencial
Agora chegamos à análise principal. O último script que será executado é o differential_expression.R, que é chamado o comando Rscript scripts/differential_expression.R pela regra differential_expression_analysis do Makefile. O código está apresentado abaixo.
#!/usr/bin/env Rscriptlibrary(edgeR)
library(limma)
library(RColorBrewer)
library(mixOmics)
library(HTSFilter)rawCountTable <- read.table('../read_count/merged.csv', header=TRUE, row.names=1, sep=',')
sampleInfo <- read.table("../samples_groups", header=TRUE, row.names=1, sep=',')dgeFull <- DGEList(rawCountTable, group=sampleInfo$group)
dgeFull$sampleInfo <- sampleInfo
dgeFull$samples$condition <- relevel(sampleInfo$group, ref = "CONTROL")pseudoCounts <- log2(dgeFull$counts+1)boxplot(pseudoCounts, col="gray")dge <- calcNormFactors(dgeFull)
dge <- estimateCommonDisp(dge)
dge <- estimateTagwiseDisp(dge)dgeTest <- exactTest(dge)
filtData <- HTSFilter(dge)$filteredData
dgeTestFilt <- exactTest(filtData)
hist(dgeTest$table[,"PValue"], breaks=50)
hist(dgeTestFilt$table[,"PValue"], breaks=50)resFilt <- topTags(dgeTestFilt, n=nrow(dgeTest$table))sigDownReg <- resFilt$table[resFilt$table$FDR<0.05,]
sigDownReg <- sigDownReg[order(sigDownReg$logFC),]
sigUpReg <- sigDownReg[order(sigDownReg$logFC, decreasing=TRUE),]write.csv(sigDownReg, file="sigDownReg.csv")
write.csv(sigUpReg, file="sigUpReg.csv")plotSmear(dgeTestFilt,
de.tags = rownames(resFilt$table)[which(resFilt$table$FDR<0.01)])volcanoData <- cbind(resFilt$table$logFC, -log10(resFilt$table$FDR))
colnames(volcanoData) <- c("logFC", "negLogPval")plot(volcanoData, pch=19)
O que este script faz é ser as contagens de leituras do arquivo merged.csv e as informações de grupos (tratamentos) do arquivo sample_groups. Depois disso, para cada amostras são calculados fatores de normalização ( calcNormFactors(...)), de modo que seja possível comparar as quantificações nos diferentes datasets. Já as funções estimateCommonDisp(...) e estimateTagwiseDip(...) analisam a distribuição dos dados nos diferentes genes, seja usando verossimilhança (likehood) e distribuição binomial negativa, seja por estatística Bayesiana. A partir da distribuição aferida, é feita então a análise gene a gene entre os diferentes grupos pela função exactTest(...), e os valores são filtrados pelo pacote HTSFilter. Por fim, os genes com expressão diferencial são então extraídos com a função topTags(...).
Ao final do processo, os genes com expressão diferencial são salvos nos arquivos sigDownReg.csv e sigUpReg.csv . Ambos arquivos contêm os mesmos genes diferencialmente expressos, mas no primeiros a ordem é dos down-regulated par aos up-regulated, e no segundo na ordem oposta. Além disso, também temos um arquivo PDF gerado pelo R com os gráficos da análise. Dentre os gráficos gerados, um dos mais pertinentes para a análise de RNA-Seq diferencial é o o Volcano Plot (mostrado abaixo).
Ainda poderíamos realizar uma análise de enriquecimento funcional para avaliar os processos biológicos que estão mais ativados em cada condição, mas deixarei esta para o próximo capítulo 😄.
