# Relatório Final
**Disciplina:** Bioinformática II: Análises de transcriptoma
**Nome:** Mirella Cunha Melaré
**Prof:** Daniel Guariz Pinheiro
### Análise transcriptômica
Análise transcriptômica é o estudo do transcriptoma -conjunto de transcritos de RNA expressos em circunstâncias ou tecidos específicos. Esse tipo de análise tem sido utilizado para avaliar diferenças de expressão em condições contrastantes, em tecidos, ou mesmo em períodos diferentes. Também tem sido utilizado para avaliar o perfil transcriptômico em diagnósticos de doenças, e descobertas de novos biomarcadores, ou até na medicina personalizada avaliando a resposta a drogas (1).
Uma das abordagens para análise de transcriptomas são os microarrays, que analisa um conjunto de sequências hibridizadas em uma superficie. Entretanto para avaliar a expressão gênica por essa técnica, é necessário a pre-determinação de genes que vão servir como "sonda" para hibridização. Atualmente, a tecnologia de sequenciamento de RNA tem sido mais utilizada, pois se obtêm um grande conjunto completo de dados (conjunto de transcritos) de uma só vez. Para isso é extraído o RNA da amostra desejada, e em seguida é produzido o cDNA a partir da mesma que é então submetido ao processo de sequenciamento de nova geração (NGS). Há diferentes plataformas de sequenciamento, mais dentre as mais atuais, todas realizam o sequenciamento de milhões de fragmentos de DNA em paralelo. Comparando a tecnologia de microarray, o RNA-Seq pode mensurar tanto genes super-expressos quanto genes de baixa abundância, necessitando de pouco material genético (1,2).
A bioinformática entre nesse contexto a fim de realizar análises de quantificação dos níveis de expressão gênica, identificação de splicing alternativos e as isoformas geradas devido a esse processo, novos transcritos, e genes colapsados. Para isso ela atua na montagem desses fragmentos seja por meio do mapeamento das reads contra um genoma de referência, disponível em bancos de dados públicos, seja montando de novo, ou seja sem a utilização de um genoma. Cada base é sequenciada várias vezes, quanto mais vezes as bases forem sequenciadas, maior a profundidade, maior a acurácia nos dados (1,2).
Para a prática da disciplina em questão, foi realizada uma análise transcriptômica a partir de dados simulados.
Como não foi possível utilizar dados publicos ou reais, devido a espaço no servidor, foi gerado dados a partir de genes escolhidos a partir do banco de dados púplico NCBI - National Center for Biotechnology Information.
### Simulação dos dados
A primeira etapa consistiu da escolha de um organismo para se trabalhar. Foi escolhido a espécie *Gallus gallus* (frango), pois é o organismo com o qual trabalho na pós graduação.
O genoma do *Gallus gallus* apresenta um tamanho de 1043Mb, uma média de conteúdo GC de cerca de 41% e aproximademente 49681 proteínas já identificadas (NCBI).
Foram escolhidos 5 genes da espécie escolhida para simulação dos dados de sequenciamento:
* Mtor: alvo mecanicista da rapamicina quinase
* EIF4EBP1: Fator de iniciação de tradução de eucariotos
* RPS6KB1: Proteina Ribossomal S6 quinase
* FLRT3: Proteina transmembrana rica em repetições de leucina
* LRRN2: Proteina neuronal rica em repetições de Leucina.
Para buscar esses genes diretamente no NCBI via linha de comando, foi necessário a instalação do pacote E-direct:
Os comentários estão a frente do "#" pois assim não são executados.
```
# Para baixar o arquivo de uma página da web, foi utilizado o seguinte comando na linguagem perl:
/bin/bash
perl -MNet::FTP -e \
'$ftp = new Net::FTP("ftp.ncbi.nlm.nih.gov", Passive => 1);
$ftp->login; $ftp->binary;
$ftp->get("/entrez/entrezdirect/edirect.tar.gz");'
# Para descompactar o pacote:
gunzip -c edirect.tar.gz | tar xf -
# Para remover o arquivo compactado e liberar espaço:
rm edirect.tar.gz
# Para adicionar o caminho no $PATH :
builtin exit
export PATH=${PATH}:$HOME/edirect >& /dev/null || setenv PATH "${PATH}:$HOME/edirect"
# Para a instalação do pacote edirect, é preciso executar o seguinte script:
./edirect/setup.sh
```
Para baixar os genes selecionados do NCBI, foi utilizado o comando `esearch` e esses já foram sendo adicionados ao arquivo transcriptoma.fa:
```
# O primeiro comando tem apenas um simbolo ">" enquanto que os seguintes tem ">>". Isso faz com que se evite sobreescrever, pois se não no final, o arquivo iria conter apenas um gene e não todos os 5.
esearch -db nucleotide -query XM_417614.6 | efetch \
-format fasta > transcriptoma.fa
esearch -db nucleotide -query XM_424384.6 | efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_025141968.1| efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_015283423.2 | efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_425820.6 | efetch \
-format fasta >> transcriptoma.fa
```
- Geração dos arquivos de abundância
Foram criados manualmente dois arquivos utilizando um editor de textos para representação de duas condições biológicas contrastantes. Os mesmos números de acesso utilizados para a geração do transcriptoma, foram colocados na coluna 1. Na segunda coluna foi dada uma proporção fictícia de cada transcrito (entre 0 e 1), levando em consideração que os genes podem estar em diferentes quantidades em cada condição biológica. Os arquivos gerados estão representados abaixo:
abundanceAsimu.txt:
XM_417614.6 0.50
XM_424384.6 0.05
XM_025141968.1 0.10
XM_015283423.2 0.20
XM_425820.6 0.15
abundanceBsimu.txt
XM_417614.6 0.05
XM_424384.6 0.05
XM_025141968.1 0.15
XM_015283423.2 0.50
XM_425820.6 0.25
- Arquivos fastq
Para as análises transcriptomicas é necessário que as reads estejam em formato fastq. Foram geradas 2 réplicas de cada condição. E para cada reṕlica e cada condição foram realizados os passos a seguir similarmente feitos com a condição A réplica 1:
-- Geração de sequências:
Foram geradas 250000 sequências de fragmentos aleatórios utilizando o transcriptoma gerado inicialmente e levando em consideração a abundância contida no arquivo de abundancia referente respectiva condição biológica. As mesmas possuem 300pb com desvio padrão de 30 pb.
```
generate_fragments.py -r transcriptoma.fa \
-a ./abundance_A.txt \ #arquivo de abundância
-o ./tmp.frags.r1 \ #arquivo gerado no processo
-t 25000 \ #quantidade de fragmentos a serem gerados
-i 300 \ # Tamanho em pares de bases dos segmentos gerados
-s 30 # Desvio padrão do tamanho dos fragmentos
```
-- Renomeação das sequências:
As sequências foram renomeadas para "SAMPLEA" juntamente com um número hexadecimal sequencial começando por “0000000001” e as bases foram posicionadas em uma única linha (até 1000 pb).
O comando SED tem como função substituições. Nesse caso, substituirá a descrição adicional.
```
cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \
-if FASTA \
-of FASTA \
-p SAMPLEA1 \
-w 1000 | \
sed 's/^>\(\S\+\).*/>\1/' \
> ./frags_A1.fa
```
-- Simulação:
Foi realizada uma simulação considerando leituras “paired” 2x “151” bp , com parâmetros configurados ([https://www.ebi.ac.uk/goldman-srv/simNGS/runfiles5/151cycPE/s_4_0099.runfile]) usando simNGS:
```
cat ./frags_A1.fa | simNGS -a \
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-p paired \
/usr/local/bioinfo/simNGS/data/s_4_0099.runfile \
-n 151 > ./SAMPLEA1.fastq 2> SAMPLEA1.err.txt
```
-- Geração das reads 1 e 2:
As reads geradas foram desintercaladas as em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq no diretório "raw"
```
#Para criar um diretório é utilizado o comando mkdir:
mkdir -p ./raw
#Separação das reads:
deinterleave_pairs SAMPLEA1.fastq \
-o ./raw/SAMPLEA1_R1.fastq \
./raw/SAMPLEA1_R2.fastq
```
```
# Remover arquivos desnecessários:
rm -f ./tmp.frags.r1.1.fasta ./frags_A1.fa ./SAMPLEA1.fastq ./SAMPLEA1.err.txt
#Verificar o número de sequências geradas
wc -l ./raw/SAMPLEA1_R1.fastq
wc -l ./raw/SAMPLEA1_R2.fastq
```
#### Quantidade de sequências geradas por read de cada condição:
100000/4 = 25000
- O comando acima deu o valor de 100000, mais cada sequencia ocorre no arquivo a cada 4 linhas, por isso foi dividido por 4.
#### Número de sequências submetidas ao processamento:
No total 200000 (25000x8) sequências foram geradas e enviadas para processamento.
### Obtenção das referências
Foram utilizadas as sequências fasta e o arquivo de anotação no formato gff do genoma de referência do organismo escolhido (*Gallus gallus*)
```
# Foi criado outro diretório denonimanado "ref"
mkdir ./refs
# cd é para entrar no diretório
cd ./refs
# Para baixar ambos arquivos foi utilizado o comando wget + o respectivo link:
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.gff.gz
# Para descompactar os arquivos:
gunzip *.gz
```
### Limpeza dos arquivos referência
Foi feito uma limpeza dos arquivos de referência baixados. Para o arquivo em formato FASTA foi feito uma limpeza do cabeçalho, excluindo tudo depois do espaço na linha de identificação, utilizando o script cleanfasta.sh:
```
# Todos os scripts devem conter essa primeira linha para funcionar:
#!/bin/bash
# É preciso passar o arquivo fasta na linha de comando para execução do script.
# Se não passar, vai aparecer na tela uma mensagem dizendo que está faltando o arquivo fasta.
# Se não existir o arquivo, aparecerá a mensagem dizendo que o arquivo não foi encontrado.
infile=$1
if [ ! ${infile} ]; then
echo "Missing input fasta file"
exit
fi
if [ ! -e ${infile} ]; then
echo "Not found input fasta file (${infile})"
exit
fi
# Eliminação das descrições necessárias na linha de indentificação. Após o espaço, tudo será eliminado.
sed 's/^>\(\S\+\).*/>\1/' ${infile}
```
- Permissão para execução:
Após ser ter feito o script, deve-se sempre dar permissão para executá-lo:
`chmod a+x cleanfasta.sh`
- Execução do script:
Para executar o script basta colocar o "./" na frente do seu nome mais os argumentos necessários para sua execução. Nesse caso é necessário ser repassado somente o arquivo no formato FASTA que será submetido ao processo de limpeza. Por último é colocado o nome do arquivo de saída (limpo):
`./cleanfasta.sh GCF_003025095.1_Triha_v1.0_genomic.fna > genome.fa`
- Ajustes do arquivo GFF:
Para os ajustes do arquivo de anotação do genoma de referência, foi executado o script fixNCBIgff.sh que se encontra no repositório Github.
Similarmente a execução do script anterior, é necessário passar como argumento a anotação genômica no formato gff e por último é colocado o arquivo de saída.
`fixNCBIgff.sh GCF_003025095.1_Triha_v1.0_genomic.gff genome.gff`
O script fixNCBIgff.sh possui dependências:
- fixGFFgenestruct.pl
- fixGFFdupID.pl
- sortGFF.pl
- Testes com GFF
Para a conversão de GFF versão 3 para um arquivo GTF:
`gffread genome.gff -g genome.fa -T -o genome.gtf`
Para obtenção de um arquivo com a sequência FASTA do trasncriptoma:
`gffread genome.gff -g genome.fa -w transcriptome.fa`
Para obtenção de um arquivo com a sequência FASTA do proteoma:
`gffread genome.gff -g genome.fa -y proteome.fa`
### Testando alinhadores
Alinhamento é o preciso de identificação de onde as reads são similares á sequência utilizada como referência.
Por exemplo, na figura abaixo pode-se observar o alinhamento das reads em um genoma de referência formando por último uma sequência consenso.
Foram testados 3 alinhadores: Bowtie2, TopHat e STAR.
- Bowtie2: O Bowtie 2 é uma ferramenta rápida e eficiente para alinhar as reads provenientes de sequenciamento contra um genoma de referência. O programa primeiramente indexa o genoma com um índice FM para facilitar o mapeamento. O alinhamento do que o programa realiza é por deafult to tipo end-to-end, ou seja, alinhamento global(3).
```
# Indexação do genoma de referência: É preciso passar o arquivo no formato fasta do genoma de referência que será indexado e o nome do índice. Como output tem-se vários arquivos com o nome genome seguido de números, similarmente a um índice de livros.
bowtie2-build -f genome.fa genome
```
```
# Alinhamento das reads contra o genoma de referência
# -x para se colocar o genoma de referência.
# No parâmetro -q deve-se colocar as reads com os números -1 ou -2 se referindo a reads 1 e reads 2 respectivamente.
# -S é tipo do arquivo que será o output, nesse caso será no formato SAM.
bowtie2 -x ./ref/genome -q -1 ./raw/SAMPLEA1_R1.fastq -2 ./raw/SAMPLEA1_R2.fastq -S bowtie2-test.sam
```
```
# Conversão de SAM para BAM
# BAM é um formato binário
samtools view -b bowtie2-test.sam > bowtie2-test.bam
# Para realizar a ordenação por coordenadas do arquivo:
samtools sort bowtie2-test.bam -o bowtie2-test-sorted.bam
# Para indexar o arquivo já ordenado:
samtools index bowtie2-test-sorted.bam
```
- TopHat: O software TopHat, utiliza o bowtie2 para fazer o alinhamento das sequências, e aquelas sequências que não foram mapeadas, ele tenta alinha utilizando informações de splicing alternativo, pois ele identifica potenciais exons (4).
Utilizando as informações de mapeamento pelo bowtie2, o TopHat constrói um banco de dados de possíveis junções de splicing e em seguida realiza o mapeamento das reads contra essas junções para confirmá-las (4).
```
Tophat2 --num-threads 2 \ #Número de processadores
--GTF ./refs/genome.gtf \ #Arquivo de anotação no formato gtf
--microexon-search \ #Buscar por alinhamentos incidentes para micro-exons
--coverage-search \ # Habilita a cobertura baseado na busca por junções
--read-edit-dist 2 \ #Alinhamento de reads que apresenta mais que 2 distancias serão descartadas
--b2-very-sensitive \ #Tipo de alinhamento:Muito sensitivo
./refs/genome \
./raw/SAMPLEA1_R1.fastq ./raw/SAMPLEA1_R2.fastq
```
- O alinhador STAR foi utilizado no script de análise transcriptômica e será discutido mais adiante.
**A partir desse ponto foi feito o script preprocess3.sh de processamento de dados, que envolve primeiro o controle de qualidade dos dados, seguido do processamento dos mesmos.**
O Script preprocess3 se encontra depositado no repositório github. Para acessá-lo clique [aqui](https://github.com/MMelare/Aula-de-Bioinform-tica-II-An-lise-de-Transcriptomas/blob/master/preprocess3.sh)
### Controle de qualidade
Um primeiro passo a se fazer com os dados provenientes do RNASeq é a checagem da qualidade das reads obtidas, pode-se observar se há algum problema com os dados, pois é importante saber antes de iniciar as análises.
Um programa que possui essa função e que foi utilizado para checar a qualidade dos dados, é o FASTQC.
O FASTQC fornece um rápido resumo de onde pode haver problemas nos dados. Esse resumo é dado também em uma interface gráfica, onde se tem graficos e tabelas com os resultados (5).
Foi realizado 2 controles de qualidade utilizando o FASTQC. Um antes do processamento e outro após o processamento, a fim de checar se os dados estão apropriados para prosseguir com as análises.
Como exemplo, segue-se uma das análises principais a serem consideradas realizadas pelo FASTQC: Qualidade da sequência por base.
- Pre-processamento:
- Pós-processamnto:
**Discussão:** Pode-se observar que antes do processamento haviam algumas bases que apresentavam um score baixo, ou seja de baixa qualidade atingindo o quadrante vermelho, score inferior a 20. Após o processamento essas bases de baixa qualidade foram removidas e as restantes possuem um bom score de qualidade acima de 20.
Um score de 20 significa um erro a cada 100 bases, enquanto que um score de 30 significa 1 erro a cada 1000 bases.
### Processamento de dados
Após a checagem da qualidade dos dados, os mesmos foram submetidos a etapa de processamento, isto é, a remoção de bases de baixa qualidade e trimagem de sequencias remanescentes de adaptadores.
Para a etapa processamento de dados 2 programas foram utilizados: Atropos e PrinSeq.
- Atropos: O Atropos é um software que realiza a trimagem de sequencias de adaptadores remanescentes do preparo das bibliotecas para sequenciamento. O programa apresenta dois algoritmos: O primeiro denominado "insert-match" é baseado na correspondecia das reads, assim as mesmas são alinhadas (Read1 contra a fita reversa da read2) e o que sobrar são os adaptadores que são então trimados. As sequências que não forem trimadas nessa primeira etapa são submetidas ao segundo algortimo denominado "adapter-match". Como o próprio nome diz, é baseado no alinhamento entre os adaptadores que são então trimados (6).
Depois da trimagem dos adaptadores as reads são submetidas ao processamento pelo programa PrinSeq.
- PrinSeq: O Prinseq é responsável, nesse caso, pela remoção de bases de baixa qualidade, qualquer sequência que ainda tenha permanecido de adaptadores e outros artefatos (7).
Após o processamento dos dados utilizando o PrinSeg, os dados foram novamente submetidos a checagem de controle de qualidade e em seguida ao alinhamento das reads contra o genoma de referência.
As etapas a seguir estão contidas no script rnaseq-ref.sh. Esse script primeiramente executa do script anterior preprocess3.sh, e depois realiza o alinhamento, a montagem, e análise de quantificação e expressão diferencial dos genes. Os comentários dos comandos se encontram em ambos os scripts.
**- Número de sequências que foram processadas:**
SAMPLEA1.atropos_final.prinseq_1.fastq: **24194**
SAMPLEA1.atropos_final.prinseq_2.fastq: **24194**
SAMPLEA1.atropos_final.prinseq_1_singletons.fastq: **293**
SAMPLEA1.atropos_final.prinseq_2_singletons.fastq: **0**
SAMPLEA2.atropos_final.prinseq_1.fastq: **24164**
SAMPLEA2.atropos_final.prinseq_2.fastq: **24164**
SAMPLEA2.atropos_final.prinseq_1_singletons.fastq: **283**
SAMPLEA2.atropos_final.prinseq_2_singletons.fastq: **0**
SAMPLEB1.atropos_final.prinseq_1.fastq: **24215**
SAMPLEB1.atropos_final.prinseq_2.fastq: **24215**
SAMPLEB1.atropos_final.prinseq_1_singletons.fastq: **290**
SAMPLEB1.atropos_final.prinseq_2_singletons.fastq: **0**
SAMPLEB2.atropos_final.prinseq_1.fastq: **24222**
SAMPLEB2.atropos_final.prinseq_2.fastq: **24222**
SAMPLEB2.atropos_final.prinseq_1_singletons.fastq: **277**
SAMPLEB2.atropos_final.prinseq_2_singletons.fastq: **0**
**Total de reads processadas: 194733**
A parte de alinhamento das reads contra um genoma de referência, montagem e análise de expressão gênica diferencial foi realizada utilizando o script rnaseq-ref.sh. O mesmo também tem em sua execução o script preprocess3.sh.
Para acessar o script rnaseq-ref.sh clique [aqui](https://github.com/MMelare/Aula-de-Bioinform-tica-II-An-lise-de-Transcriptomas/blob/master/rnaseq-ref.sh)
### STAR
O programa STAR - Spliced Transcripts Alignment to a Reference, realiza o alinhamento das reads contra o genoma de referência. O STAR difere de outros alinhadores desenvolvidos por utilizar como estratégia alinhamento contendo splicing alternativo (9).
Para cada read que o STAR alinha, ele busca pela sequência mais longa que corresponde a uma ou mais localizações no genoma de referência. Cada parte da read que é mapeada separadamente é chamado de seed. Em seguida, o STAR irá procurar por porções das reads que não foram mapeadas tentando identificar splicing alternativo e então juntar essas possíveis regiões e mapeá-las contra o genoma (9).
- Alinhamento A1:
24486 de fragmentos alinhados, desses:
24193 são paired-ends;
293 são single-end.
98.78% do total de fragmentos foram mapeados.
- Alinhamento A2:
24447 de fragmentos alinhados, desses:
24164 são paired-ends;
283 são single-end.
98.81% do total de fragmentos foram mapeados.
- Alinhamento B1:
24505 de fragmentos alinhados, desses:
24220 são paired-ends;
290 são single-end.
98.79% do total de fragmentos foram mapeados.
- Alinamento B2:
24497 de fragmentos alinhados, desses:
24193 são paired-ends;
277 são single-end.
98.86% do total de fragmentos foram mapeados.
**Total de fragmentos que foram mapeados: 97935**
### Montagem
Montagem do transcriptoma é o processo de idenficação dos transcritos e suas variantes que estão sendo expressas sob determinada condição ou em um determinado tecido ou periodo. O objetivo da montagem é reconstruir o conjunto completo de sequências de todos os transcritos. Um transcriptoma pode ser definido então, como um conjunto de sequencias contiguas que representam regiões transcritas (RNAs). A montagem de transcriptômica pode ser realizada de duas maneiras: baseando-se em uma referência ou não (montagem *de novo*), embora uma combinação de ambos também possa ser utilizado (10). O esquema abaixo representa uma fluxograma das etapas de montagem a serem realizadas a partir da obtenção das reads, provenientes do sequenciamento, quando se tem ou não um genoma para ser utilizado como referência. Ambas as abordagens serão discutidas a seguir.
### Montagem do transcriptoma baseando-se em uma referência:
Essa abordagem é utilizada quando se tem acesso a um genome sequenciado completo do organismo que está sendo estudando, ou seja, o transcriptoma é reconstruido por meio de um mapeamento contra sequencias previamente conhecidas. Quando há uma referência de boa qualidade, essa estratégia é recomendada, pois possui alta sensibilidade. A acurácia dessa abordagem depende de um alinhamento correto, sendo importante levar em consideração os splicings, assim como realizado pelo alinhador STAR. Nessa abordagem contaminantes não vem a ser um grande problema pois não se alinham contra o genoma de referência. Há dois montadores que são mais populares: Cufflinks e StringTie (10).
#### Cufflinks e StringTie
Cufflinks e StringTie são montadores de transcriptomas a partir de dados de RNA-Seq, utilizando um genoma de referência. Esses montadores têm a vantagem de inferir splicings alternativos, criando uma tabela somente com as junções de splicing. Durante a aula utilizamos ambos os programas para realizar a montagem do transcriptoma. Ambos estão no script rnaseq-ref.sh, onde foi comentado os parãmetros utilizados (11,12).
Cufflinks Também é o nome de um conjunto de ferramentas utilizadas para quantificar o nivel de expressão gênica, os quais foram utilizados nessa análise. Os seguintes passos foram reaizado utilizando a montagem do StringTie. Embora na aula ambas foram feitas, a análise apresentada foi refeita com outro organismo, e nesse caso foi utilizado somente a montagem do StringTie. As ferramentas estão descritas a seguir e comentadas no mesmo script.
**-Cuffcompare**
Após a montagem do transcritoma, foi realizada um comparação entre a mesma e o genoma de referência para avaliar a qualidade de montagem. Quando não se tem um genoma de referência, essa comparação pode ser feita utilizando organismos próximos filogeneticamente.
**-Cuffmerge**
Como output das montagens, foram gerados vários arquivos uma para réplica e condição. Essas montagens foram então unidas em um único arquivo no formato gtf, para posteriormente se realizar a análise de expressão gênica.
**-Cuffquant**
Essa ferramenta computa os perfis de expressão dos genes e transcritos, originando arquivos de abundancia no formato cxb.
**-Cuffdiff**
Foi realizada a comparação dos níveis de expressão dos genes e transcritos utilizando o cuffdiff. Essa ferramenta fornece várias informações, entre elas, quais genes estão superexpressos ou inibidos em determinada condição biológica. Além disso, relata também genes que sofreram diferencialmente splicing alternativo. Também mede as diferenças entre isoformas, regiões cds, entre outros.
**-Cuffnorm**
Por último foi realizada uma normalização dos níveis de expressão para que fiquem na mesma escala para futuras análises.
### Avaliação da expressão gênica diferencial
A tabela abaixo foi realizada buscando o ID do gene referente ao ID dado no arquivo de saída do programa cuffdiff. Em seguida foram adicionados as respectivas proporções dadas para cada gene dos arquivos de abundância gerados. A expressão diferencial, log2(fold change) é dada no arquivo de expressão diferencial gênica, saída do programa cuffdiff.
Podemos observar que os valores de log2 - fold change, estão, em sua maioria, de acordo com os valores propostos nos arquivos de abundancia. O gene mais diferencialmente expresso foi o XM_417614.6, sendo mais expresso na condição A. O gene XM_424384.6, embora tenha dado uma diferença de expressão significante, não está concordante com os valores do arquivo da abundancia, pode ter ocorrido algum erro durante o processo de simulação e análise dos dados. Os outros genes também estão de acordo com as proporções do arquivo de abundância.
### Informações adicionais aos resultados de expressão gênica diferencial
Para fundir dados ou informações de 2 arquivos diferentes, baseando-se em coluna comum (usando um IDentificador comum) entre os dois arquivos formatados como tabelas, utilizaremos os seguintes passos:
Queremos realizar a junção de 2 arquivos:
- gene_exp.diff: Saída do programa cuffdiff
- Uma tabela do NCBI que contem descrições dos genes do organismos selecionado (*Gallus gallus*).
Essa tabela do NCBI foi baixada diretamente via linha de comando, utilizando o ID taxonomico da espécie:
as ferramentas E-utilities do NCBI.
```
esearch -db gene -query "9031[Taxonomy ID]" \
| efetch -format tabular > ./ref/gene_result.txt
```
A partir dessa tabela selecionamos as colunas de interesse: GeneID, Symbol e Description, que foram salvos em um novo arquivo denominado gene_info.txt:
`cut -f 3,6,8 ./ref/gene_result.txt > ./ref/gene_info.txt`
O cufflinks/cuffmerge, durante o seu processamento, pode fundir genes considerados inicialmente distintos, que ficam unidos por vírgula. Para o script Merger.R que realizará a fusão dos arquivos funcionar, foi necessário utilizar a função abaixo para separar os identificadores que ficaram unidos:
```
splitteR.R --x="./output/cuffdiff/gene_exp.diff" \
--col.x="gene" \
--by.x="," \
--out="./output/cuffdiff/gene_exp.diff.splitted.txt"
```
Em seguida, foi realizada a fusão dos dados, dois 2 arquivos, o gene_info.txt e o gene_exp.diff.splitted.txt.
Foi utilizado o parâmetro –all.x para que todas as linhas do arquivo gene_exp.diff.splitted.txt estejam no resultado final, no arquivo a ser denominando gene_exp.diff.spplitted.merged.txt, mesmo que não tenha uma linha correspondente no arquivo gene_info.txt. Essa última situação pode acontecer quando o cufflinks cria um novo mapeamento gênico, que não corresponde a nenhum outro contido no GFF/GTF original, usado como guia para alinhamento e montagem do transcritoma.
Para essa fusão o script abaixo foi utilizado com a seguinte linha de comando:
```
MergeR.R --x="./output/cuffdiff/gene_exp.diff.spplitted.txt" \
--all.x \
--y="./refs/gene_info.txt" \
--by.x="gene" \
--by.y="Symbol" \
--out="./output/cuffdiff/gene_exp.diff.spplitted.merged.txt"
```
Abaixo está a imagem ta tabela resultante, onde temos o simbolo, os Genes e os outros resultados da tabela gerada no cuffdiff. Não está na imagem, mas mais a frente tem-se a descrição dos genes.
### Montagem *de novo*
Nessa abordagem a reconstrução do transcriptoma é feita sem utilizar um genoma de referência. Essa técnica pode ser mais complexa pois há reads que alinham em multiplos locais do genoma. Quanto menor a read, maior a chance dela se alinhar em multiplos sitios. Além disso, pode-se ter um grande número de isoformas o que também pode comprometer a acurácia da montagem. A complexidade dessa abordagem requer uma demanda maior do uso da bioinformática. O montador atual mais utilizado é o Trinity, o qual se baseia na construção de grafos de Bruijn, onde os kmers-derivados das reads são montados com base em sua sobreposição (k-1) (10).
A figura abaixo representa um esquema do grafo de Bruijn, onde "a" representa as reads, "b" são os kmers de 5 pb originados a partir das reads, "c" Grafo de Bruijn, onde leva-se em conta a sobreposição de k-1 bases para ir juntando as sequencias certas e contiguas, formando a sequencia consenso (d).
Outros montadores utilizam o algoritmo OLC - Overlap-layout-consensus, onde a sobreposição se da alinhando read a read, o que requer muito tempo e muitos recursos computacionais. Logo, grafos de Bruijn são preferidos para realizar montagem *de novo*. O programa Trinity, também reconstrói transcritos a partir de splicing alternativo e sequencias parálogas fazendo uma clusterização de contigs sobrepostos e gerando um grafo de Bruijn para cada cluster separadamente, para então montar a sequencia consenso (10).
Foi feita uma aula prática de montagem *de novo* onde fui utilizado o programa Trinity. Para essa motagem foi utilizado o script rnaseq-denovo.sh, também comentado.
Para acessar o script rnaseq-denovo.sh clique [aqui](https://github.com/MMelare/Aula-de-Bioinform-tica-II-An-lise-de-Transcriptomas/blob/master/rnaseq-denovo.sh)
Como input foram utilizadas as reads processadas no formato fastq. No output passado será criado um diretório de montagem.
Para execução do script foi utilizado o seguinte comando:
`./rnaseq-denovo.sh ./output/processed/prinseq ./output`
O script rnaseq-denovo.sh inicialmente realizou um processo inicial de renomear as reads para que os cabeçalhos contenham o sufixo /1 e /2, respectivamente para as reads R1 e R2. Isso foi feito tanto para as sequencias paired end, como para as singletons.
O Trinity também tem uma opção para executar a montagem utilizando como guia os alinhamentos das reads contra o genoma referência. Os alinhamentos utilizados são os obtidos com o programa STAR. Para executar este modo de montagem, foi criado o script rnaseq-ref-trinity.sh.
Esse script está disponível no github, para acessá-lo clique [aqui](https://github.com/MMelare/Aula-de-Bioinform-tica-II-An-lise-de-Transcriptomas/blob/master/rnaseq-ref-trinity.sh)
Agora, o input então são os arquivos de alinhamento no formato bam.
Na primeira etapa deste script antes de realizar a montagem é realizado uma junção de todos os arquivos de alinhamento em um só utilizando a ferramenta samtools. O script também está comentado os parametros utilizados.
Para a execução do script o seguinte comando foi utilizado:
`./rnaseq-ref-trinity.sh ./output/star_out_final ./output`
- Para avaliação das montagens foi criado um banco de dados denominado "Reftrans" que é resultado da indexação do transcriptoma gerado no inicio com os genes escolhidos:
```
makeblastdb -title "RefTrans" \ #Titulo do banco de dados
-dbtype nucl \ # Tipo do banco: nucleotideos
-out RefTrans \ #output
-in transcriptoma.fa \ #input
> makebleastdb.log.out.txt \
2> makeblastdb.log.err.txt
```
Após a execuçã do comando anterior os seguintes arquivos foram gerados:
- RefTrans.nhr
- RefTrans.nin
- RefTrans.nsq
Em seguida foi realizada uma comparação da montagem montada com o trinity e o banco de dados criado, utilizando a ferramenta blastn:
- Trinity de novo:
```
blastn -max_hsps 1 \ #Número maximo de HSPs - segmentos que apresentam um alto score de alinhamento.
-max_target_seqs 1 \ # número de sequências alinhadas que se deseja manter
-num_threads 8 \ #Número de processadores
-query ./output/trinity_assembled/Trinity.fasta \ #Sequência que se deseja comparar
-task blastn \ #Tarefa: Será realizado o blastn
-db ./RefTrans \ #Subject: banco de dados a ser utilizado
-out ./Trinity_x_RefTrans.txt \ # Nome do output
-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \ #Opção de visualização do alinhamento
> ./Trinity_x_RefTrans.log.out.txt \
2> ./Trinity_x_RefTrans.log.err.txt
```
Para visualizar os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos foi utilizado o comando abaixo, que seleciona somente as colunas de interesse, no caso a 1 e 2.
`cut -f 1,2 ./Trinity_x_RefTrans.txt`
Uma parte do arquivo gerado está representado na figura abaixo:
Foram encontradas 5 ocorrências de sequencias na base de dados, para a contagem foi utilizado o seguinte comando:
`cut -f 2 ./Trinity_x_RefTrans.txt | sort | uniq -c`
As ocorrências estão representadas na imagem abaixo:
**Discussão:** Observando a tabela resultante, observamos que as ocorrências, os genes mostrados, estão de acordo com os genes escolhidos para ser criado os arquivos de abundância na primeira etapa.
- Trinity Genome Guided:
O mesmo processo foi realizado novamente entretanto para a montagem o qual foram utilizados os alinhamentos obtidos com o programa STAR:
```
blastn -max_hsps 1 \
-max_target_seqs 1 \
-num_threads 8 \
-query ./output/trinity_GG_assembled/Trinity-GG.fasta \
-task blastn \
-db ./RefTrans \
-out ./Trinity-GG_x_RefTrans.txt \
-evalue 1e-5 \
-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
> ./Trinity-GG_x_RefTrans.log.out.txt \
2> ./Trinity-GG_x_RefTrans.log.err.txt
```
- Nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos:
`cut -f 1,2 ./Trinity-GG_x_RefTrans.txt`
A imagem abaixo é uma representação do arquivo gerado:
- Ocorrências de sequências da base de dados (subject) foram encontradas, na imagem abaixo, se encontram os resultados.
`cut -f 2 ./Trinity-GG_x_RefTrans.txt | sort | uniq -c`
**Discussão:** Analisando o arquivo gerado, novamente observamos que os genes encontrados estão de acordo com os genes escolhidos para criação dos arquivos de abundancia na primeira etapa.
Assim, ambas as montagens realizadas estão de acordo com o proposto o inicio da análise.
### Genes diferencialmente expressos na montagem de novo
Para obtenção dos valores de Log Fold-Change e os p-valores ajustados referentes às comparações pareadas entre SAMPLEA e SAMPLEB utilizamos scripts auxiliares do programa Trinity. Além disso foi utilizado o script run-DESeq2.R para executar uma análise comparativa pareada de expressão gênica diferencial.
Essa análise está contida no rnaseq-trinity-abundance.sh. O script rnaseq-trinity-abundance.sh está disponivel no github, para acessá-lo clique [aqui](https://github.com/MMelare/Aula-de-Bioinform-tica-II-An-lise-de-Transcriptomas/blob/master/rnaseq-trinity-abundance.sh)
Para sua execução foi utilizada a linha de comando abaixo:
```
rnaseq-trinity-abundance.sh ./output/renamed/ \
./output/trinity_assembled/ \
./output
```
O resultado pode ser observado utilizando o comando abaixo:
```
cat ./output/abundance/DEG/SAMPLEA-SAMPLEB.txt | column -s$'\t' -t | less -S
```
**Discussão:** Na tabela assim podemos observar a expressão de cada gene em casa condição e réplica. Além disso tem-se o Log Fold Change que mostra a expressão genica diferencial, ou seja a diferença de expressão nas condições A e B. Esse resultado foi utilizado na proxima etapa, a qual envolve a seleção dos genes diferencialmente expressos no programa RStudio.
### Seleção dos genes diferencialmente expressos no R
- Para carregar a biblioteca xlsx para lidar com arquivos do Excel, foi utilizado o comando "library":
```
library("xlsx")
```
- Para carregar a biblioteca da função para gerar o HeatMap - heatmap.2() foi utilizado o seguinte comando:
```
library("gplots")
```
- Similarmente para carregar a biblioteca com a função para criar a paleta de cores:
```
library("RColorBrewer")
```
- Equivalente ao comando "pwd" no linux:
```
getwd()
/home/mmelare
```
- Equivalente ao comando "cd" no linux:
```
setwd("~/")
O "~" foi substituido pelo caminho do diretório desejado:
> setwd("/state/partition1/mmelare/sim/")
```
```
getwd()
[1] "/state/partition1/mmelare/sim"
```
- Para carregar o arquivo texto separado por TAB ("\t") com cabeçalho e sem transformar strings em fatores - Arquivo gerado com o script do programa Trinity para calcular abundancia:
```
deg.df <- read.delim(file="./output/abundance/DEG/SAMPLEA-SAMPLEB.txt",
sep="\t",
header=TRUE,
stringsAsFactors = FALSE)
```
- Para exibir somente as 2 primeiras linhas para checagem foi utilizado o comando abaixo:
```
head(deg.df,2)
X SAMPLEA2 SAMPLEB1 SAMPLEB2 SAMPLEA1 logFC PValue FDR
1 TRINITY_DN0_c0_g1 5179.003 4829.623 5359.392 2533.695 -0.4021695 0.4317459 0.7580055
2 TRINITY_DN1_c0_g1 9721.495 5678.906 6629.650 3761.165 0.1310359 0.8047183 0.8407752
```
- Para saber a dimensão do data.frame, onde o primeiro valor refere-se à quantidade de linhas e o segundo valor à quantidade de colunas foi utilizado o comando abaixo:
```
dim(deg.df)
[1] 5 8
```
- Em seguida foi selecionado um subconjunto de linhas do data.frame "deg.df" e colocado em um outro data.frame denominado "subset.deg.df" contendo somente os genes que possuírem logFC >= 2 ou logFC <= -2, além de possuírem um p-valor corrigido ("FDR") <= 0.05 subset.deg.df <-
```
subset(deg.df, (((logFC >= 2) | (logFC <= -2) ) & (FDR <=0.05)) )
```
```
dim(subset.deg.df)
[1] 1 8
```
- Para visualizar somente as duas primeiras linhas do arquivo gerado:
```
head(subset.deg.df,2)
X SAMPLEA2 SAMPLEB1 SAMPLEB2 SAMPLEA1 logFC PValue FDR
3 TRINITY_DN2_c0_g1 3653.513 419.0021 684.2491 7511.826 3.340099 2.287929e-09 1.143964e-08
```
- Depois foi criado um vetor de nomes de colunas selecionadas, ou seja, colunas que não são "X", "logFC", "PValue" ou "FDR",as quais, sabemos previamente que contém os valores de expressão.
Os nomes das colunas foram ordenados depois de secionados com a função "setdiff", a qual obtém um conjunto de valores (vetor) coma a diferença entre dois conjuntos de
valores: o de todas as colunas ("colnames(subset.deg.df)") menos os nomes das colunas indesejadas.
```
sel_columns <- sort(
setdiff( colnames(subset.deg.df),
c("X", "logFC", "PValue", "FDR")
)
)
```
- Para criar uma nova matriz "expression_data" que contém somente as colunas selecionadas foi utilizado o comando abaixo:
```
expression_data <- as.matrix( subset.deg.df[,sel_columns] )
```
- Foi atribuido para os nomes das linhas da matriz o conteúdo da coluna "X", onde sabemos previamente que contém o identificador dos genes, utilizando o comando abaixo:
```
rownames(expression_data) <- subset.deg.df$X
```
- Para visualizar somente as 2 primeiras linhas:
```
head(expression_data,2)
SAMPLEA1 SAMPLEA2 SAMPLEB1 SAMPLEB2
TRINITY_DN2_c0_g1 7511.826 3653.513 419.0021 684.2491
```
- Foi utilizado a seguinte função para gravar o data.frame em um arquivo Excel, na planilha de nome "DEGS_SAMPLEA-SAMPLEB":
```
write.xlsx(subset.deg.df,
file="./output/abundance/DEG/SAMPLEA-SAMPLEB.xlsx",
sheetName="DEGS_SAMPLEA-SAMPLEB"
)
```
- Para criar uma paleta personalizada de 299 cores do vermelho ao verde, passando pelo amarelo, foi utilizado o seguinte comando:
```
my_palette <- colorRampPalette(c("red", "yellow", "green"))(n = 299)
```
- Para definir as quebras das cores manualmente para transição de cores:
```
col_breaks = c(seq(-1,0,length=100), # for red
seq(0.01,0.8,length=100), # for yellow
seq(0.81,1,length=100)) # for green
```
- Foi gerada uma imagem de tamanho 5 x 5 polegadas utilizando a função png:
```
png("./output/abundance/DEG/heatmap_DEGS_SAMPLEA-SAMPLEB.png", # cria arquivo do tipo PNG for the heat map
width = 5*300, # 5 x 300 pixels de largura
height = 5*300, # 5 x 300 pixels de altura
res = 300, # 300 pixels por polegada
pointsize = 8) # tamanho da fonte
```
- Para criar um heatmap dos genes diferencialmente expressos:
```
heatmap.2(expression_data, # a matriz com os valores de expressão
main = "Correlation", # Título do HeatMap
density.info="none", # desabilita o gráfico de densidade dentro da legenda
trace="none", # desabilita as linhas dentro do HeatMap
margins =c(12,12), # definiação das margens no entorno do gráfico
col=my_palette, # nome do objeto contendo a paleta de cores criada anteriormente
breaks=col_breaks, # pontos de quebra para a transição de cores
dendrogram="both", # desenhar dendrograma para linhas e colunas
distfun = function(x) as.dist(1-cor(t(x))), # distância baseada em correlação
hclustfun = function(x) hclust(x, method="centroid") # método de ligação pelo centróide
)
```
Como foram selecionados somente os genes que apresentam logFC >= 2 ou logFC <= -2, somente um gene foi selecionado e não foi possivel gerar o heatmap.
```
Erro: "Error in heatmap.2(expression_data, main = "Correlation", density.info = "none", :
`x' must have at least 2 rows and 2 columns"
```
Na tentativa de gerar o heatmap a partir da tabela deg.df, um heatmap totalmente verde foi gerado, não mostrando as diferanças de expressão.
Para gerar o heatmap, deve-se conter pelo menos 2 genes.
### Conclusão
Concluindo, os resultados das análises de genes diferencialmente expressos estão todas concordantes com os arquivos de abundancia criados na primeira etapa.
Com os outputs gerados a partir da montagem utilizando o genoma de referência, pode-se observar que as proporções dos genes diferencialmente expressos estão concordantes com as proporções colocadas nos arquivos de abundancia.
Com os outputs gerados a partir da montagem de novo, pode-se observar que os genes que se mostraram diferencialmente expressos também são os mesmos dos genes escolhidos para criação dos arquivos de abundancia.
Assim, ambas as montagens se mostraram adequadas e com resultados concordantes ao esperado.
A discplina de Bioinformática II me proporcionou um maior conhecimento sobre o sistema Linux, Sequenciamento de Nova Geração e Análises de transcriptômica.
Foi visto a diferença da montagem *de novo* e da montagem utilizando um genoma de referência.
Além disso, tivemos um primeiro contato na criação de scripts, o que facilita nas análises pois assim não é necessários fazer um comando para cada amostra e manualmente, principalmente utilizando as estruturas *for* e *find*.
Embora ainda tenho muito que aprender, abriu minha cabeça a novas opções de análises e como e quando utilizar cada programa. Foi discutido no decorrer da disciplina os parâmetros que podem ser utilizados para cada programa, e isso facilitou a entender melhor cada processo. Mesmo o heatmap não ter sido gerado, foi possivel aprender como deve ser feito e as condições necessárias para que seja gerado. A disciplina me forneceu um maior conhecimento em uma área que eu possuia um defict imenso, e irá contribuir para as análises do meu projeto de pós graduação.
### Referências bibliográficas
1.Blumenberg M. Introductory Chapter: Transcriptome Analysis, Transcriptome Analysis. IntechOpen. 2019 DOI: 10.5772/intechopen.85980. Disponível em: < https://www.intechopen.com/books/transcriptome-analysis/introductory-chapter-transcriptome-analysis>. Acesso em 26 de novembro de 2019.
2.Behjati S, Tarpey PS. What is next generation sequencing?. Arch Dis Child Educ Pract Ed. 2013; 98: 236–8.
3.Langmead B, Salzberg S. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012, 9:357-359.
4.Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics.2009. doi:10.1093/bioinformatics/btp120
5.FASTQC Manual. Disponível em: <https://dnacore.missouri.edu/PDF/FastQC_Manual.pdf>. Acesso em 27 de novembro de 2019.
6.Didion JP, Martin M, Collins FS. Atropos: specific, sensitive, and speedy trimming of sequencing reads. PeerJ. 2017; 5:3720
7.PrinSeq Manual. Disponível em: <http://prinseq.sourceforge.net/manual.html> Acesso em 27 de novembro de 2019.
8.STAR Manal. Disponível em <http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STAR.posix/doc/STARmanual.pdf>
9.Introduction to RNA-Seq using high-performance computing. Disponível em <https://hbctraining.github.io/Intro-to-rnaseq-hpc-O2/lessons/03_alignment.html>
10.Moreton J, Izquierdo A, Emes RD. Assembly, assessment, and availability of de novo generated eukaryotic transcriptomes. Front Genet. 2015;6:36
11.Cufflinks.The Cufflinks RNA-Seq workflow. Disponível em: <http://cole-trapnell-lab.github.io/cufflinks/manual/> Acesso em 27 de novembro de 2019.
12.Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.
Sign in with Wallet
Connect another wallet