# **Relatório de Bioinformática Aplicada II: Análise de Transcriptomas**
## 1. Introdução
Ao longo da disciplina foi realizada a simulação, processamento, montagem e a análise de um transcriptoma. Tendo como objetivo a demonstração de todas as etapas necessárias, assim como os diferentes softwares utilizados em cada uma das etapas. Dessa forma, esse relatório tem como objetivo a descrição dos trabalhos práticos realizados, bem como a comparação de cada software utilizado.
## 2. Simulação e análise dos dados simulados
Para melhor descrever os processos realizados, esse relatório será dividido nas etapas de construção da simulação, pré-processamento, montagem com genoma referência, montagem *de novo* e análise dos resultados utilizando o software R. Para realização das práticas, foi utilizado um servidor composto por 4 nodes, com sistema operacional CentOS (versão 6.6), sediado no campus de Jaboticabal da Universidade Estadual Paulista (UNESP). O acesso ao servidor foi realizado de forma remota, utilizando o protocolo SSH, em uma máquina virtual Ubuntu (versão 16.04), construída com o software VirtualBox (versão 6.1).
### 2.1. Construção do transcriptoma simulado
Para a construção do transcriptoma, foi selecionada a espécie de inseto *Halyomorpha halys* (Stål) (Hemiptera: Pentatomidae). Essa espécie de percevejo é considerada praga em lavouras de algodão, milho e soja, causando elevados danos econômicos principalmente nas regiões agrícolas dos Estados Unidos, Canada e Europa. Embora essa espécie tenha grande importância agrícola, estudos ômicos são recentes e restritos. Dessa forma, não há grande quantidade de informações (SPARKS et al. 2014, SPARKS et al. 2020). Ambos transcriptoma (SPARKS et al. 2014) e genôma (SPARKS et al. 2020) estão disponíveis no GeneBank do NCBI, nesse sentido, há informações suficientes para realização das atividades práticas.
#### 2.1.1. Obtenção das sequencias
Para a construção do transcriptoma, inicialmente, foi realizada uma busca de transcritos de genes identificados no transcriptoma de *H. halys* disponível no GeneBank. Apenas sequencias não caracterizadas de mRNA (XM_) e uma sequencia de RNA não codante (NC_) foram selecionadas.
XM_014426076.1 General odorant-binding protein 56d isoform X1
XM_014426082.1 General odorant-binding protein 56d isoform X2
XM_014423650.1 Neural Wiskott-Aldrich syndrome protein
XM_024359338.1 Zinc finger protein 711-like
XM_014438980.2 Cytoplasmic tRNA 2-thiolation protein 1
XM_024358251.1 Uncharacterized protein LOC106692498 isoform X2
XM_024358252.1 Uncharacterized protein LOC106692498 isoform X3
XM_024358250.1 Uncharacterized protein LOC106692498 isoform X1
XM_024358253.1 Uncharacterized protein LOC106692498 isoform X4
XM_024358254.1 Uncharacterized protein LOC106692498 isoform X5
NC_022611.1 Halyomorpha halys virus isolate Beltsville
Para a obtenção desses genes foi criado um arquivo ACCS.txt contendo os números de acesso citados acima (sem as descrições). Além disso, foram criados dois arquivos para simular a abundância dos transcritos de cada gene em duas condições biológicas hipotéticas diferentes. Esses arquivos foram denominados de abundance_A.txt e abundance_B.txt.
abundance_A.txt
XM_014426076.1 0.1
XM_014426082.1 0.1
XM_014423650.1 0.2
XM_024359338.1 0.025
XM_014438980.2 0.055
XM_024358251.1 0.14
XM_024358252.1 0.02
XM_024358250.1 0.02
XM_024358253.1 0.02
XM_024358254.1 0.02
NC_022611.1 0.3
abundance_B.txt
XM_014426076.1 0.1
XM_014426082.1 0.1
XM_014423650.1 0.27
XM_024359338.1 0
XM_014438980.2 0.058
XM_024358251.1 0.002
XM_024358252.1 0.06
XM_024358250.1 0.07
XM_024358253.1 0.07
XM_024358254.1 0.07
NC_022611.1 0.2
#### 2.1.2. Construção da simulação e organização com links simbólicos
Após o preparo dos arquivos ACCS.txt, abundance_A.txt e abundance_B.txt, os arquivos foram checados para remoção de espaços entre os nomes mantendo apenas separações por tabulação. O script [*mksim.sh*](https://github.com/ericmar-santos/bioinfoclass/blob/main/mklink.sh) foi executado no mesmo diretório onde estava presente os arquivos txt informados acima. Esse script faz o download dos acessos e os utiliza em conjunto com os arquivos de abundância para gerar fragmentos aleatórios, simulando dados de sequenciamento NGS paired-end. O output do script é um diretório /raw com arquivos fastq com duas repetições por condição biológica (SAMPLEA e SAMPLEB) e duas reads (paired-end) por repetição, totalizando 8 arquivos fastq.
Para organizar os arquivos de output, foi executado o script [*mklink.sh*](https://github.com/ericmar-santos/bioinfoclass/blob/main/mklink.sh). Esse script pega os arquivos fastq localizados no diretório de input e cria um link simbólico alterando os nomes para CONTROL (no lugar de SEMPLEA) e TEST (no lugar de SEMPLEB), em um diretório indicado como output. Dessa forma, facilitando a identificação dos arquivos para melhor organização.
### 2.2. Pré-Processamento
Os dados de sequenciamento, inicialmente, devem passar por uma etapa de pré-processamento para controle de qualidade das reads, remoção da sequencia de adaptadores e contaminantes. Para realização dessa etapa, foi executada uma pipeline denominada [*preprocess4.sh*](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh). Essa pipeline utiliza o FastQC (ANDREWS, 2010) para controle de qualidade das reads, o Atropos (DIDION et al. 2017) para trimmagem de adaptadores, o PrinSeq (SCHMIEDER e EDWARDS, 2011) para trimmagem de qualidade, o Bowtie2 (LANGMEAD e ALZBERG, 2012) para alinhamento das reads com sequencias de contaminantes em um banco de dados local e o SAMtools (LI et al, 2009) e PullSeq (THOMAS, 2015) para remoção das sequencias contaminantes.
#### 2.2.1. Controle de qualidade
A maioria dos dados brutos de sequenciamento de nova geração são obtidos no formato fastq. Esse formato apresenta duas informações essenciais, a sequencia de nucleotideos das reads e o respectivo valor de qualidade em score phred (phred 33 ou phred 64) de cada uma das bases sequenciadas. De forma geral, scores phred acima de 30 são considerados como qualidade alta. O software FastQC demonstra uma série de características das reads, permitindo a realização do controle de qualidade. Dentre essas caracteristicas, o software produz um gráfico demonstrando a distribuição de scores phred em relação a posição das bases nas sequências das reads. Esse gráfico permite uma visualização geral da qualidade das sequências das reads.
A fim de avaliar a estratégia de pré-processamento utilizada, o controle de qualidade com o FastQC foi feito antes e depois da execução do script [preprocess4.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh). Antes do pré-processamento, é possivel observar que a maioria dos arquivos apresentam uma grande amplitude de qualidade ao longo das reads, com muitas reads com qualidade baixa nas posições mais próximas da terminação (Figura 2).
Figura 2. FastQC report demonstrando a distribuição da qualidade das reads antes do pré-processamento
Após o pré-processamento, onde é feita as trimmagens dos adaptadores e trimmagem de qualidade, é possível observar claramente que houve uma redução na amplitude de qualidade, seendo a maior parte das reads com nucleotídeos com alto phred score na maioria das posições, havendo uma redução na qualidade nas terminações, no entanto em apenas alguns casos a qualidade alcançou scores baixos (Figura 3).
Figura 3. FastQC report demonstrando a distribuição da qualidade das reads após o pré-processamento
Além disso, foi feito o controle da qualidade das reads singletons (Figura 4). Essas reads tiveram seu par excluido, devido a baixa qualidade. É possível observar que elas apresentam qualidade intermediária a baixa em sua maioria, e apresentam comprimento reduzido, devido a trimmagem de qualidade.
Figura 4. FastQC report demonstrando a distribuição da qualidada das reads singletons após o pré-processamento
#### 2.2.2. Trimmagem dos adaptadores
A trimmagem das sequencias dos adaptadores foi realizada com a ferramenta Atropos. O Atropos é um software que se utiliza de uma abordagem hibrida para remoção dos adaptadores, onde é feito um alinhamento semi-global da sequencia dos adaptadores com as reads e também é verificada a sobreposição das reads para identificar a posição dos adaptadores (DIDION et al. 2017). Dessa forma, a trimmagem é realizada duas vezes, sendo a primeira de acordo com o pareamento entre as sequencia dos insertos (insert-match) e a segunda é de acordo com o paremento entre os adaptadores (adapter-match). Após as duas trimmagens, as reads obtidas foram contatenadas nos mesmos arquivos correspondentes. Os parâmetros utilizados para trimmagem dos adaptadores das reads simuladas estão descritos no script [preprocess4.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh).
#### 2.2.3. Trimmagem de qualidade
A segunda trimmagem realizada foi beasada na qualidade do *base-calling* indicada pelo score Phred. Essa etapa foi realizada com o software PrinSeq. O PrinSeq é um software de acesso público, que pode ser usado para filtrar, reformatar ou trimmar dados de sequências genômicas e metagenômicas. Essa trimmagem foi realizada com o objetivo de remover as bases com baixo score Phred de cada read. A baixa qualidade pode ser decorrente de artefatos do processo de sequenciamento, o que pode resultar em erros nas etapas seguintes. Os parâmetros utilizados podem ser encontrados no script [preprocess4.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh).
A fim de comparar a quantidade de reads antes e após as trimmagens, foram usados os seguintes comandos:
```bash=
grep -c "@SAMPLE" CONTROL*
CONTROL1_R1.fastq:25000
CONTROL1_R2.fastq:25000
CONTROL2_R1.fastq:25000
CONTROL2_R2.fastq:25000
grep -c "@SAMPLE" TEST*
TEST1_R1.fastq:25001
TEST1_R2.fastq:25001
TEST2_R1.fastq:25001
TEST2_R2.fastq:25001
grep -c "@SAMPLE" CONTROL*.fastq
CONTROL1.atropos_final.prinseq_1.fastq:24470
CONTROL1.atropos_final.prinseq_1_singletons.fastq:321
CONTROL1.atropos_final.prinseq_2.fastq:24470
CONTROL1.atropos_final.prinseq_2_singletons.fastq:0
CONTROL1.atropos_final.prinseq_singletons.fastq:321
CONTROL2.atropos_final.prinseq_1.fastq:24423
CONTROL2.atropos_final.prinseq_1_singletons.fastq:308
CONTROL2.atropos_final.prinseq_2.fastq:24423
CONTROL2.atropos_final.prinseq_2_singletons.fastq:0
CONTROL2.atropos_final.prinseq_singletons.fastq:308
grep -c "@SAMPLE" TEST*.fastq
TEST1.atropos_final.prinseq_1.fastq:24383
TEST1.atropos_final.prinseq_1_singletons.fastq:353
TEST1.atropos_final.prinseq_2.fastq:24383
TEST1.atropos_final.prinseq_2_singletons.fastq:0
TEST1.atropos_final.prinseq_singletons.fastq:353
TEST2.atropos_final.prinseq_1.fastq:24404
TEST2.atropos_final.prinseq_1_singletons.fastq:359
TEST2.atropos_final.prinseq_2.fastq:24404
TEST2.atropos_final.prinseq_2_singletons.fastq:0
TEST2.atropos_final.prinseq_singletons.fastq:359
```
É possível observar nos resultados acima que houve uma redução de 2,12% na quantidade de reads em CONTROL1_R1 e CONTROL1_R2 e 1,28% de todas as reads nesses dois arquivos resultaram em reads singletons (sem um par correspondente). Em TEST1_R1 e TEST1_RS houve uma redução de 2,31% na quantidade de reads e 1,23% das reads dando origem a reads singletons. No entanto, houve um aumento drástico na qualidade das reads em cada um dos arquivos (Figura 3) e uma qualidade dentro do aceitável entre as reads singletons (Figura 4).
#### 2.2.4. Remoção das sequencias contaminantes
Para remoção de sequencias contaminantes foi criado um diretório /refs/contaminants e foi colocado um arquivo fasta contendo a sequencia NC_022611.1 (Halyomorpha halys virus). Inicialmente, foi feito alinhamento da sequencia contaminante com os arquivos de reads obtidos na etapa anterior, com o software Bowtie2. Na sequência, foi utilizado o comando cat para pegar todo o conteudo de output do alinhamento, o samtools view pega a identificação de todas as reads que alinharam com o banco de dados de contaminantes e cria um arquivo "_excluded.txt" para cada arquivo de reads. Após a obtenção da identificação das reads que devem ser excluidas, foi utilizado o sofware PullSeq para remoção das reads contaminantes em cada um dos arquivos. Como resultado, foram obtidos arquivos livres do contaminante contido no banco de dados. Esse processo está descrito no script [preprocess4.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh).
### 2.3. Montagem com genôma referência
O script [preprocess4.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/preprocess4.sh) foi incorporado na pipeline [rnaseq-ref.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/rnaseq-ref.sh). Essa pipeline foi utilizada para montagem do transcriptoma simulado utilizando um transcriptoma de referência. Nessa pipeline, também foram incorporados os softwares Tophat2 (KIM et al. 2013) e STAR (DOBIM et al. 2013) para alinhamento das reads paired-end e singletons contra o genoma de referência, o software SamTools para fundir os resultados dos alinhamentos paired-end e single-end, os softwares Cufflinks (TRAPNELL et al. 2013) e StringTie (PERTEA et al. 2015) para montagem do transcriptoma, e por fim, comandos construidos no Cufflinks serão utilizados para quantificação da expressão gênica. Portanto, serão avaliadas combinações de dois softwares alinhadores e dois montadores.
#### 2.3.1. Obtenção e preparo do genôma referência
Para montagem do transcriptoma com um genôma de referência, inicialmente é necessária a obtenção, limpeza, padronização e indexação dos arquivos fasta (sequêcia do genôma) e gff (anotação do genôma)do genôma referência. Para obtenção dos arquivos, foi utilizado os seguintes comandos:
```bash=
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/696/795/GCF_000696795.2_Hhal_2.0/GCF_000696795.2_Hhal_2.0_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/696/795/GCF_000696795.2_Hhal_2.0/GCF_000696795.2_Hhal_2.0_genomic.gff.gz
```
Os arquivos obtidos foram descompactados com o *gunzip* da seguinte maneira:
```bash=
gunzip *.gz
```
Após a descompactação, foi feita a limpeza dos cabeçalhos das sequencias contidas no arquivo fasta com o script [cleanfasta.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/cleanfasta.sh). Após a limpeza, apenas os codigos de acesso das sequencias são mantidos (*e.g.* NW_020110170.1). O arquivo genome.fa obtido após a limpeza foi utilizado para realizar a indexação com o seguinte comando:
```bash=
bowtie2-build genome.fa genome
```
Para certificação que a indexação foi bem suscedida, foi utilizado o seguinte comando:
```bash=
bowtie2-inspect genome | grep '^>'
```
O arquivo gff foi submetido a alterações para padronização e remoção de informações não essenciais através da utilização do script [FixNCBIgff.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/fixNCBIgff.sh). O arquivo gff obtido foi utilizado para gerar um arquivo gtf, utilizando o seguinte comando:
```bash=
gffread genome.gff -g genome.fa -T -o genome.gtf
```
O arquivo genome.gtf também contêm coordenadas de genes no genôma, assim como o genome.gff. No entanto, alguns programas reconhecem apenas arquivos gtf. Além disso, o arquivo gtf(Gene Transfer Format) possui maior refinamento das especificações gff(General Feature Format), sendo específico para descrever estruturas gênicas (gene,RNA,CDS).
#### 2.3.2. Montagem
A montagem do transcriptoma utilizando um genôma de referência é obtida a partir do alinhamento das reads com o genôma de referência (feito por alinhadores) e posterior montagem dos transcritos consenso (feito por montadores). Com a montagem do transcriptoma é possível identificar novas isoformas ou novos gênes e/ou utilizar a informação do número de transcritos de cada gene identificado para quantificar a expressão gênica. A expressão gênica pode ser utilizada para descrever o perfil de expressão de um dado organismo, sob diferentes condições biológicas. Possibilitando a identificação dos genes que tem sua expressão afetada por essa condição biológica.
Durante o processo de montagem, existem dois erros comuns que devem ser mitigados. O primeiro é a fragmentação, resultado de uma baixa cobertura de reads em um transcrito. O segundo é a formação de quimeras, resultado da fusão errônea de partes de transcritos diferentes devido a ocorrência de sobreposição ambigua entre esses transcritos. Esses erros podem resultar na identificação de genes/isoformas erradas ou na redução da precisão da quantificação da expressão gênica.
**Alinhamento:**
Como citado anteriormente, foram utilizados os softwares alinhadores Tophat2 e STAR. Os parâmetros utilizados podem ser verificados no script [rnaseq-ref.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/rnaseq-ref.sh). Os resultados do alinhamento estão descritos na Tabela 1 a seguir.
**Tabela 1**. Número de reads alinhadas com TopHat2 e STAR
|Amostras |Alinhador| Input Reads |Uniquely mapped Reads | Multi-mapped Reads|
|--------|---------|-----------|-----------|---|
|CONTROL1|TopHat2 |17707 | 15915 |32 |
|CONTROL1|STAR |17707 | 12200 |38 |
|CONTROL2|TopHat2 | 17718 | 15928 |28 |
|CONTROL2|STAR | 17718 | 12246 |45 |
|TEST1 |TopHat2 | 20032 | 18147 |27 |
|TEST1 |STAR | 20032 | 13964 |45 |
|TEST2 |TopHat2 | 20036 | 18300 |24 |
|TEST2 |STAR | 20036 | 18912 |54 |
Os resultados demonstram que o alinhador Tophat2 mapeou uma maior quantidade de reads em uma unica posição do genôma nas amostras CONTROL1, CONTROL2 e TEST1. Apenas a amostra TEST2 resultou em maior quantidade de reads mapeadas em uma unica posição com o alinhador STAR. De forma geral, o alinhador STAR mapeou uma maior quantidade de reads em mais de uma posição, em comparação com o TopHat2. Dessa forma, com os parâmetros utilizados no script [rnaseq-ref.sh](https://github.com/ericmar-santos/bioinfoclass/blob/main/rnaseq-ref.sh) e com esse grupo de dados simulados, o alinhador TopHat2 se mostrou superior.
**Montagem:**
Após o alinhamento, são obtidas informações das reads que foram mapeadas no genoma de referência. Para junção dessas reads em transcritos, são utilizados softwares montadores. Foram utilizados os softwares Cufflinks e StringTie, a partir dos resultados de alinhamento obtidos com ambos os alinhadores. Tendo sido obtido os resultados descritos na Tabela 2.
|Amostras |Softwares|Nº Transcritos Montados | Nº de Isoformas Mapeadas |
|--------|---------|-----------|--------|
|CONTROL1|TopHat + CuffLinks | 286.731| 28.926|
|CONTROL1|TopHat + StringTie | |
|CONTROL1|STAR + CuffLinks | |
|CONTROL1|STAR + StringTie | |
|CONTROL2|TopHat + CuffLinks |286.731 |28.926|
|CONTROL2|TopHat + StringTie | |
|CONTROL2|STAR + CuffLinks | |
|CONTROL2|STAR + StringTie | |
|TEST1 |TopHat + CuffLinks |286.731 |28.926|
|TEST1 |TopHat + StringTie | |
|TEST1 |STAR + CuffLinks | |
|TEST1 |STAR + StringTie | |
|TEST2 |TopHat + CuffLinks |286.731 |28.926|
|TEST2 |TopHat + StringTie | |
|TEST2 |STAR + CuffLinks | |
|TEST2 |STAR + StringTie | |
O script Cuffmerge do Cufflink, ou a função merge do StringTie foram usadas para utilizar os resultados de ambos os alinhadores para construção de transcriptomas concenso.
#### 2.3.3. Análise da expressão diferencial
### 2.4. Montagem *de novo*
#### 2.4.1. Montagem do transcriptoma
#### 2.4.2. Análise da expressão diferencial
### 2.5. Análise com o software estatístico R
## 3. Conclusões
## Referências bibliográficas
ANDREWS, S. FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. 2010. Disponível em: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
DIDION, J. P.; MARTIN, M.; COLLINS, F. S. Atropos: Specific, sensitive, and speedy trimming of sequencing reads. PeerJ, v. 2017, n. 8, p. 1–19, 2017. Disponível em: https://doi.org/10.7717/peerj.3720
LANGMEAD, B.; SALZBERG, S. L. Fast gapped-read alignment with Bowtie 2. Nature Methods, v. 9, n. 4, p. 357–359, 2012. Disponível em: https://doi.org/10.1038/nmeth.1923
LI, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics, v. 25, n. 16, p. 2078–2079, 2009. Disponível em: https://doi.org/10.1093/bioinformatics/btp352
SCHMIEDER, R.; EDWARDS, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics, v. 27, n. 6, p. 863–864, 2011. Disponível em: https://doi.org/10.1093/bioinformatics/btr026
SPARKS, M. E. et al. Transcriptome of the invasive brown marmorated stink bug, Halyomorpha halys (stål) (Heteroptera: Pentatomidae). PLoS ONE, v. 9, n. 11, 2014. Disponível em: https://doi.org/10.1371/journal.pone.0111646
SPARKS, M. E. et al. Brown marmorated stink bug, Halyomorpha halys (Stål), genome: Putative underpinnings of polyphagy, insecticide resistance potential and biology of a top worldwide pest. BMC Genomics, v. 21, n. 1, p. 1–26, 2020. Disponível em: https://doi.org/10.1186/s12864-020-6510-7
THOMAS, B. C. PullSeq. 2015. Disponível em: https://github.com/bcthomas/pullseq.git
KIM, D. et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology, v. 14, n. R36, 2013. Disponível em: https://doi.org/https://doi.org/10.1186/gb-2013-14-4-r36
DOBIN, A. et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics, v. 29, n. 1, p. 15–21, 2013. Disponível em: https://doi.org/10.1093/bioinformatics/bts635
TRAPNELL, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature protocol, v. 7, n. 3, p. 562–578, 2013. Disponível em: https://doi.org/10.1038/nprot.2012.016.Differential
PERTEA, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature biotechnology, v. 33, n. 3, p. 290–295, 2015. Disponível em: https://doi.org/10.1038/nbt.3122
Sign in with Wallet
Connect another wallet