Relatório Final: Bioinformática II - Análises de transcriptomas

Aluna: Jéssica Luana Souza Cardoso

Introdução

Desde a descoberta da estrutura do DNA em 1953, esforços tem sido direcionados para o desenvolvimento de abordagens que visam a obtenção da sequência de DNA. Essas metodologias tem evoluido continuamente nas últimas décadas, graças a por exemplo ao desenvolvimento da reação em cadeia da polimerase, ao aumento da disponibilidade de enzimas de alta qualidade que agem sob os ácidos nucleicos, entre outros (Levy and Myers, 2016).
Na primeira técnica de sequenciamento de DNA, sequenciamento de Sanger, conhecido como sequenciamento por síntese, a sequência de DNA era determinada, um nucleotideo por vez, em função de uma fita molde de DNA. Anos mais tarde, a Applied Biosystems, automatizou esse processo, surgindo o Projeto Genoma Humano logo após. Nesse projeto, foi necessário a contribuição de varios laboratórios ao redor do mundo, e demorou cerca de uma década para finalizá-lo (Kumar et al., 2019).
Com o desenvolvimento da performance da tecnologia do sequenciamento, emergiu o sequenciamento de nova geração (NGS), abordagem que tem dominado o cenário científico atualmente, pois tem reduzido o tempo e custo consideravelmente, e aumentado a quantidade de dados obtidos a partir do processo quando comparado com a metodologia de Sanger. Enquanto a metodologia de Sanger apresenta a capacidade de produzir uma sequência a partir de um molde por reção, o NGS promove de milhões a bilhões de reações individuais simultaneamente. Atualmente, o genoma humano pode ser sequenciado dentro de 3 dias apenas a um custo muito inferior (Kumar et al., 2019;Levy and Myers, 2016).
Posteriormente, surgiu a necessidade de identificar os elementos funcionais do genoma, o transcriptoma, que é definido como a soma de todos os transcritos de RNA. A informação de um organismo está contida no DNA do genoma, e esse é expresso por meio da transcrição. O transcriptoma captura um cenário instantaneo dos transcritos presentes em uma célula sob determinada condição. Embora a tecnologia de microarray seja também utilizada, ela tem perdido espaço para o sequenciamento de RNA - RNASeq, pois enquanto que a primeira abordagem quantifica apenas um conjunto de genes predeterminados, a segunda consegue capturar praticamente todas as sequencias presentes, sendo muito útil para mensurar a expressão gênica de um organismo sob condições, tecidos e tempos diferentes (Anamika et al., 2016;Lowe et al.,2017).
A análise de dados oriundos do processo de sequenciamento de RNA é realizada por meio da Bioinformática e inclui várias etapas: Avaliação da qualidade, processamento de dados, montagem do transcriptoma que pode ser guiado por uma referência ou não (montagem de novo), quantificação, análises estatisticas e anotação funcional (Figura 1)(Anamika et al.,2016).

Fonte: (Anamika et al.,2016)

Para o presente relatório, buscamos por mRNAs e genoma de referencia em banco de dados públicos, para a realização de uma simulação de dados de RNA-Seq e posterior análise de dados. A seguir serão descritas as etapas.

Instalação do pacote E-direct

Primeiramente foi necessário a instalação do pacote E-Direct com as ferramentas E-utilities. Esse pacote fornece acesso ao banco de dados público NCBI -National Center for Biotechnology Information, a partir de comandos em um terminal UNIX. Para a instalação foi utilizado a seguinte linha de comando:

## Entra no diretório do usuário corrente:
cd ~

## Baixa o arquivo para instalação de uma página web:
  /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");'
  
  ## O comando gunzip é utilizado para descompactar os arquivos baixados, no caso o pacote edirect:
  gunzip -c edirect.tar.gz | tar xf -
  
  ## O comando rm foi utilizado para remover o arquivo compactado:
  rm edirect.tar.gz
  
## Essa linha de comando foi realizada para adicionar o caminho no $PATH :
  builtin exit
  export PATH=${PATH}:$HOME/edirect >& /dev/null || setenv PATH "${PATH}:$HOME/edirect"
  
  ## Por último foi realizada a instalação do pacote por meio do comando acima, solicitando que o script de instalação seja iniciado:
  ./edirect/setup.sh
  

OBS: PATH é uma variável do sistema Linux que indica trajetória dos programas executaveis. Quando um programa é adicionado no PATH ele pode ser executado sem precisar indicar o caminho completo.

Construção do transcriptoma

Para a simulação de dados desse relatório foi escolhido o organismo Arabidopsis thaliana, uma pequena planta herbácea, nativa da Europa, Ásia e África (Noroeste), pertencente a família das Brassicaceae. Trata-se de uma das espécies mais utilizadas na pesquisa científica atualmente, sendo um organismo modelo de plantas em pesquisas na área da genética (Delatorre et al., 2008).

Foi consultado no NCBI os genes pertencentes a espécie Arabidopsis thaliana e foram selecionados algums mRNAs (IDs-NM_) para montagem do transcriptoma utilizando as ferramentas esearch/efetch para fazer o download e construir o arquivo fasta contendo os transcritos referência para a simulação das análises transcriptômicas.

Os genes escolhidos foram:

• NM_123951.3 - DOG1: Gene responsável pelo atraso na germinação
• NM_001345288.1 - PHOT2: Fototropina 2 - Proteina transmembrana que atua como fotoreceptor.
• NM_105102.3 - NPR1: Proteina regulatória - atua na rota de resistência adquirida ao ácido salicílico.
• NM_180129.4 - CCA1: Codifica um repressor transcricional que executa funções sobrepostas com LHY em um loop de feedback regulatório que está intimamente associado ao oscilador circadiano de Arabidopsis.
• NM_125091.4 - EIR1: Codifica uma proteína carregadora de auxina.
• NM_001337167.1 - BRM: Codifica uma ATPase de remodelação da cromatina SWI / SNF que regula positivamente a transcrição dos três genes CUC e está envolvida na formação e / ou manutenção de células de contorno durante a embriogênese. Também medeia a repressão da expressão de proteínas de armazenamento de sementes em tecidos vegetativos.

Para isso o seguinte comando foi utilizado, onde a estrtura de repetição for foi utilizada, significando que para todos os ids colocados faça o seguinte comando para montagem do transciptoma.fa.

for acc in NM_123951.3 NM_001345288.1 NM_105102.3 NM_180129.4 NM_125091.4 NM_001337167.1 ; do
    
        esearch -db nucleotide -query ${acc} | efetch \
        -format fasta >> transcriptoma.fa

done

Gerando os arquivos de abundância

Os arquivos abundance_A.txt e abundance_B.txt foram criados manualmente utilizando o comando nano, respectivamente para a condição biológica A e condição biológica B. Os mesmos ids usados na construção do arquivo transcriptoma.fa foram utilizados na coluna 1. Na coluna 2 foi colocada a proporção de cada transcrito levando em consideração que é possível haver diferenças de abundância entre essas duas condições biológicas.

OBS: As duas colunas foram separadas por TAB, sem linhas adicionais e sem espaços.

**abundance_A.txt:**

NM_123951.3     0.15
NM_001345288.1  0.25
NM_105102.3     0.05
NM_180129.4     0.05
NM_125091.4     0.50
NM_001337167.1  0

**abundance_B.txt:**

NM_123951.3     0.15
NM_001345288.1  0.05
NM_105102.3     0.25
NM_180129.4     0.05
NM_125091.4     0
NM_001337167.1  0.50

Gerando os arquivos fastq

Para cada réplica da simulação da corrida para a amostra A, levando em consideração a abundância contida no arquivo abundance_A.txt, foi utilizado o seguinte comando com o objetivo de gerar fragmentos aleatórios utilizando o transcriptoma montado anteriormente “transcriptoma.fa”. No comando o parâmetro -o significa o arquivo de saída, o -t quantos fragmentos aleatórios serão gerados, o -i o tamanho médio dos fragmentos gerados e -s o desvio padrão do tamanho dos fragmentos:

generate_fragments.py -r transcriptoma.fa \
   -a ./abundance_A.txt \
   -o ./tmp.frags.r1 \
   -t 25000 \
   -i 300 \
   -s 30

Em seguida 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). A linha de SED irá substituir a descrição adicional. Para isso foi utilizado o seguinte comando:

cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \
   -if FASTA \
   -of FASTA \
   -p SAMPLEA1 \
   -w 1000 | \
   sed 's/^>\(\S\+\).*/>\1/' \
   > ./frags_A1.fa

Posteriormente, foi realizada uma simulação considerando reads paired end de 151 bp , utilizando o comando 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

Em seguida o comando abaixo foi utilizado a fim de desintercalar as reads em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq:

deinterleave_pairs SAMPLEA1.fastq \
   -o ./raw/SAMPLEA1_R1.fastq \
      ./raw/SAMPLEA1_R2.fastq

Os arquivos desnecessários foram então removidos conforme comando abaixo, pois não os utilizaremos mais

rm -f ./tmp.frags.r1.1.fasta ./frags_A1.fa ./SAMPLEA1.fastq ./SAMPLEA1.err.txt

O mesmo processo foi repetido para réplica 2 da condição biológica A e para as réplicas 1 e 2 da condição biológica B, onde então foi utilizado o arquivo abundanceB.txt.

A fim de verificar o número de sequências geradas o comando
grep -c "@SAMPLE" foi utilizado para todas cada conjunto de reads geradas:

grep -c "@SAMPLE"       SAMPLEA1_R1.fastq = 25001
grep -c "@SAMPLE"       SAMPLEA1_R2.fastq = 25001       
grep -c "@SAMPLE"       SAMPLEA2_R1.fastq = 25001       
grep -c "@SAMPLE"       SAMPLEA2_R2.fastq = 25001       
grep -c "@SAMPLE"       SAMPLEB1_R1.fastq = 25001
grep -c "@SAMPLE"       SAMPLEB1_R2.fastq = 25001
grep -c "@SAMPLE"       SAMPLEB2_R1.fastq = 25001
grep -c "@SAMPLE"       SAMPLEB2_R2.fastq = 25001

Resultado: O total de reads que foram enviadas para processamento foram: 200008

Obtenção dos arquivos de referência para serem utilizados nas análises

Foram baixadas as sequências no formato FASTA e anotações gênicas no formato GFF da espécies Arabidopsis thaliana. Para isso na página do NCBI em genoma da espécie se encontram os arquivos. O link de ambos os arquivos foram copiados e os seguintes comandos foram utilizados:

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.fna.gz

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.gff.gz

Os arquivos baixados foram descompactados utilizando o comando gunzip, onde o asterisco significa todos os arquivos que terminam com .gz

gunzip *.gz

Limpeza do arquivo fasta:

O arquivo fasta contem em seu cabeçalho muitas informações que podem atrapalhar durante o processo de análise de dados, por isso uma limpeza do mesmo foi realizada, mantendo somente as informações necessárias: a identificação de cada sequencia:

O script cleanfasta.sh foi criado utilizando novamente a ferramenta nano.

#!/bin/bash

# Ao chamar o script é necessário passar o arquivo fasta junto:
infile=$1

# Se nenhum arquivo for passado, print na tela que está faltando arquivo fasta:
if [ ! ${infile} ]; then
        echo "Missing input fasta file"
        exit
fi

# Se o arquivo fasta não existe, print na tela que nenhum arquivo foi encontrado: 
if [ ! -e ${infile} ]; then
        echo "Not found input fasta file (${infile})"
        exit
fi

# O comando sed é utilizado para realizar substituições no arquivo:
sed 's/^>\(\S\+\).*/>\1/' ${infile}

OBS: Sempre após finalizar um script, devemos dar permissão para que esse possa ser executado. Isso é feito através do comando chmod + script:

chmod a+x cleanfasta.sh

Para o script ser executado basta adicionar um ./ na frente e dar um nome para o arquivo de saída:

./cleanfasta.sh GCF_000001735.4_TAIR10.1_genomic.fna > genome.fa

O mesmo foi realizado para o arquivo gff, entretanto utilizando o script fixNCBIgff.sh do repositório Github:

./fixNCBIgff.sh GCF_000001735.4_TAIR10.1_genomic.gff > genome.gff

Análise de dados

Para a análise de dados foi utilizado o script rnaseq.sh. Segue o link para acesso rnaseq
As alterações foram realizadas com base na espécie selecionada.
O script rnaseq.sh realiza a chamada do script preprocess.sh, o qual realiza o pré-processamento. Para acessar esse script: preprocess

O script rnaseq.sh é referente a etapa de montagem do transcriptoma utilizando um genoma de referencia. Assim, contem as etapas de alinhamento e montagem; Além disso, o script finaliza realizando uma análise de expressão gênica diferecial.

Controle de qualidade

Os sequenciadores de próxima geração utilizam um score de qualidade denonominao Phred, o qual é a probabilidade de uma chamada de base não ter sido realizada com acurácia. Baixos scores de qualidade, menores que 30, indicam bases de baixa qualidade. Isso, pode afetar a interpretação na analise de dados.
Para checar a qualidade dos dados foi utilizado o software FASTQC.
Uma vez que as reads foram checadas, elas foram processadas para remover bases de baixa qualidade, sequencias de adaptadores, e outras sequencias contaminantes.
O arquivo de entrada para o FASTQC pode ser nos formatos Fastq, SAM ou BAM, compactados ou não.
O FASTQC reporta em uma interface gráfica, estatisticas básicas das reads, qualidade das bases, sequencias super representadas, conteúdo de k-mers, conteúdo de adaptadores, duplicação das reads, entre outros. O controle de qualidade foi realizado antes e após o processamento das reads:

Exemplo: Para a amostra A1, read 1:

• Estatisticas básicas:
  – Pre processamento:
  
  
  Figura 2. Estatisticas básicas do programa fastqc para a amostra A1, read 1, antes do processamento.
  – Pós processamento:
  
  
  Figura 3. Estatisticas básicas do programa fastqc para a amostra A1, read 1, após o processamento.

Pode-se observar que após o processamento de dados o número e o comprimento das sequências variou, devido a trimagem realizada.

• Qualidade por base na sequencia:
  – Pré - processamento:
  
  
  Figura 4. Qualidade por base da sequencia do programa fastqc para a amostra A1, read 1, antes do processamento.
  – Pós - processamento:
  
  
  Figura 5. Qualidade por base da sequencia do programa fastqc para a amostra A1, read 1, após o processamento.

Pode-se observar que as bases de baixa qualidade foram trimadas durante a etapa de processamento, melhorando o score de qualidade das reads.

Processamento de dados

Na etapa de processamento de dados é realizada a trimagem de sequências de adaptadores, sequências de baixa qualidade e contaminantes.
Há três tipos de softwares de trimagem (Didion et al., 2017):

• Aqueles que dependem unicamente da correspondência da sequência do adaptador, usando alinhamento semi-global (que é a única opção disponível para leituras de extremidade única; single-end). Ex:cutadapt
• Aqueles que utilizam a sobreposição entre as reads pair-end para identificar posições iniciais dos adaptadores (match-insert). Ex: SeqPurge.
• E aqueles hibridos que utilizam ambas as abordagens.Ex:AdapterRemoval e Atropos.

A primeira etapa do processamento de dados foi realizada pelo softaware Atropos.

Atropos

O Atropos é um software de trimagem de adaptadores e de sequencias de baixa qualidade de reads sequenciadas. O processo de trimagem é importante para atenuar os efeitos de contaminação por adaptadores e erros de sequenciamento nos dados obtidos do sequenciamento de nova geração (NGS). Contaminação por adaptadores e erros de sequenciamento podem levar ao aumento de taxas de leituras não alinhadas, ou alinhamento erroneo, resultando em erros de analise downstream (Didion et al., 2017).
O Atropos foi desenvolvido a partir de outro software de trimagem, o cutadapt, que possui licença MIT e portanto tem código fonte aberto que pode ser melhorado. Ele foi melhorado em 3 pontos comparado com o cutadapt (Didion et al., 2017):

• Implementação do algoritmo insert-match, melhorando assim sua acurácia.
• Adição de suporte de multiprocessamento, comparado com o cutadapt um único processado é disponibilizado, resultando em um aumento em sua performance.
• Adição de funções como por exemplo trimagem automatizada de sequências metiladas, detecção automatizada de sequências de adaptadores em reads onde os protocolos experimentais não sao conhecidos previamente pelo analista de dados, estimação da taxa de erro de sequenciamento e geração de medidas de controle de qualidade.

São necessários 2 comandos de trimagem, relacionados aos dois tipos de trimagem que o algoritmo realiza: O primeiro de acordo com a correspondencia das reads (insert-match) e o segundo, as reads que não foram trimadas na primeira etapa são submetidas a essa segunda etapa e são trimadas de acordo a correspondencia entre os adaptadores. A figura 6 representa um esquema de como o Atropos funciona.
Ao final, as reads trimadas oriundas de ambos os algoritmos são concatenadas em um arquivo final.

Os parâmetros utilizados estão comentados no script preprocess.sh.

Figura 6. Detecção dos adaptadores para trimagem pelo Atropos. (A)Quando um fragmento é mais curto que a read, a sequencia da read irá conter sequencias de adaptadores (Sequencias azuis e roxas). (B,C)Algoritmos utilizados para detectar os adaptadores.(B) Algoritmo do tipo "adapter-match": Identifica o melhor alinhamento entre os adaptadores e o final, corresponde a read. Ou seja, é feito por correspondencia entre os adaptadores. © Algortimo do tipo "insert-match": Identifica o melhor alinhamento entre a read 1 e a fita reversa da read 2, o que sobrar são os adaptadores. (D)O alinhamento com a maior média de qualidade é escolhido. (E) Primeiramente é feito o algoritmo do tipo insert, as reads que não forem trimadas são então submetidas ao algoritmo do tipo adapter. As reads que forem muito pequenas são descartadas.
Fonte:(Didion et al., 2017)

PrinSeq

A segunda etapa do processamento foi realizada pelo software PrinSeq.
Trata-se de uma ferramenta pública, que tem a finalidade de filtrar, reformatar e trimar sequencias, fornecendo um resumo estatistico sobre seus dados, similarmente ao fastqc.
No presente relatório o PrinSeq foi utilizado para filtrar e trimar as sequencias provenientes da etapa anterior, a fim de assegurar que as reads a serem utilizadas para a analise dados não apresenta sequencias de baixa qualidade ou artefatos que podem levar a conclusões errôneas.

Os parâmetros utilizados estão comentados no script preprocess.sh.

Para contabilizar o numero de sequencias que foram processadas:

grep -c "@SAMPLE" SAMPLEA1.atropos_final.prinseq_1.fastq
23811
grep -c "@SAMPLE" SAMPLEA1.atropos_final.prinseq_1_singletons.fastq 
91
grep -c "@SAMPLE" SAMPLEA1.atropos_final.prinseq_2_singletons.fastq 
0
grep -c "@SAMPLE" SAMPLEA1.atropos_final.prinseq_2.fastq 
23811
grep -c "@SAMPLE" SAMPLEA2.atropos_final.prinseq_2.fastq 
23750
grep -c "@SAMPLE" SAMPLEA2.atropos_final.prinseq_1.fastq 
23750
grep -c "@SAMPLE" SAMPLEA2.atropos_final.prinseq_1_singletons.fastq 
108
grep -c "@SAMPLE" SAMPLEA2.atropos_final.prinseq_2_singletons.fastq 
0
grep -c "@SAMPLE" SAMPLEB1.atropos_final.prinseq_2_singletons.fastq 
0
grep -c "@SAMPLE" SAMPLEB1.atropos_final.prinseq_1_singletons.fastq 
89
grep -c "@SAMPLE" SAMPLEB1.atropos_final.prinseq_1.fastq 
23776
grep -c "@SAMPLE" SAMPLEB1.atropos_final.prinseq_2.fastq 
23776
grep -c "@SAMPLE" SAMPLEB2.atropos_final.prinseq_2.fastq 
23707
grep -c "@SAMPLE" SAMPLEB2.atropos_final.prinseq_1.fastq 
23707
grep -c "@SAMPLE" SAMPLEB2.atropos_final.prinseq_1_singletons.fastq 
97
grep -c "@SAMPLE" SAMPLEB2.atropos_final.prinseq_2_singletons.fastq 
0

Resultado: O número total de reads que foram processadas foi igual a: 190473

Montagem do transcriptoma

A montagem de transcriptomas é uma ferramenta poderosa para detecção de variações na expressão gênica e sequências entre condições, tecidos ou espécies. Entretanto a capacidade de se realizar tais analises com acurácia depende da qualidade da montagem. Por exemplo, para a detecção de isoformas e quantificação de transcritos, genes de interesses que não forem montados pode aumentar a taxa de falsos positivos e negativos (Voshall et al., 2018).
Dois dos maiores erros na montagem são fragmentação, onde um único transcrito é montado com um único ou contigs menores, e quimeras, conde um contig é montando utilizando partes de mais de um transcrito. A fragmentação geralmente ocorre por baixa cobertura da read ao longo do transcrito. Enquanto que quimeras ocorrem por sobreposições ambiguasentre as reads. Ambos os erros pode implicar na identicação de genes de maneira errônea (Voshall et al., 2018).
Além disso, transcritos podem ser colapsados em um único contig ou pode ocorrer alterações de nucleotídeos na sequência contig, devido a falta de reads sequenciadas, ou também a ambiguidade de reads apresentarem alta similaridade, por exemplo em genes heterozigotos e duplicados. Esses erros também podem levar a uma montagem incompleta (Voshall et al., 2018).
A montagem pode se dar de duas maneiras:Pode ser guiada por um genoma de referência, onde se realiza um alinhamento contra um genoma de referência para então posterior construção do transcriptoma, ou o transcriptoma pode ser montado sem ajuda de um genoma de referência, processo denominando montagem de novo.(Figura 7).

Figura 7. Montagem de transcriptomas.

Montagem de novo

Montadores de novo geram transcritos baseando-se somente nos dados de RNASeq, ou seja não utilizam um genoma de referência como guia no processo de montagem.
O montador mais popular atualmente, Trinity, realiza a montagem utilizando o grafo de Bruijn, gerados a partir de Kmers (subsequências) provenientes das reads dos dados de RNASeq. A reconstrução ocorre a partir da sobreposição dessas sequências de kmers. A limitação desse tipo de grafo, é a necessidade de um kmer para iniciar cada posição de uma sequencia original. E pode-se ter kmers muito pequenos que podem ser altamente ambiguos e serem correspondentes a mais de um transcrito, e kmers muito longos, que embora a ambiguidade seja evitada, eles não cobrem a sequencia inteira de alguns transcritos, resultando em fragmentação (Anamika et al., 2016; Voshall et al., 2018).
O montador Trinity usa 3 passos para a montagem referentes aos estágios de desenvolvimento de uma borboleta (lagarta,crisálida e borboleta) (Figura 8)(Anamika et al., 2016; Grabherr et al., 2011).

Figura 8. Montagem de novo utilizando Trinity. (a)Inchworm: constrói os conjuntos iniciais de contigs usando os kmers.
(b) Chrysalis: Agrupa esses contigs e constrói grafos de bruijn a partir desses componentes.
© Butterfly: Simplifica e resolve os grafos, gerando um conjunto final de transcritos contendo variantes de splicing e isoformas.
Fonte: (Grabherr et al., 2011).

Montagem utilizando genoma de referência

Montadores que utilizam um genome de referência pra realizar a montagem evitam a ambiguidade de kmers, mapeando os dados de RNASeq ao genoma de referência. Assim, a complexidade da montagem dos transritos é reduzida por clusterização das reads de acordo com a localização genômica de sobreposição das reads com o genoma. Esse alinhamento é realizado por ferramentas como Bowtie2, TopHat, STAR, entre outros.
Nesse metodologias, também pode ocorrer ambiguidades como por exemplo, as reads podem mapear em mais de um lugar dentro do genoma. Como o algortimo de cada alinhador lida com a escolha de qual a melhor posição pode impactar o resultado final (Voshall et al., 2018).
Alguns exemplos de montadores que utilizam um genoma de referência incluem Cufflinks e StringTie.

Alinhamento

Para se realizar um alinhamento das reads contra o genoma de referência com acurácia, deve-se considerar duas situações (Dobin et al., 2013):

• Um alinhamento de reads que contem mismatches, inserções e deleções, causadas por variações e erros de sequenciamento;
• Mapeamento de sequencias derivadas de regiões genômicas não contiguas que compreendem sequencias que sofreram splicing.
  Esses desafios de alinhamento são adicionados ao fato de que há presença de multiplas cóias de sequencias de idendicas ou relacionadas, o que dificulta o mapeamento de maneira precisa (Dobin et al., 2013).
  Varios alinhadores tem sido desenvolvidos visando superar esses desafios. Entretanto no presente relatório, foi utilizado o alinhador STAR - Spliced Transcripts Alignment to a Reference.

STAR

O alinhador STAR funciona em 2 etapas (Dobin et al., 2013):

• Primeiro ele gera um index do genoma de referência, utilizando tanto o arquivo FASTA, quanto a anotação em gff. Esse index, é como se fosse um sumário do genoma, que facilita o mapeamento das reads.
• Mapeamento das reads em formato FASTA ou FASTQ, no genoma de referência, utilizando o index criado no passo anterior.
  É gerado vários arquivos de saída: alinhamento SAM/BAM, resumo estatistico do mapeamento, junções de splicing, etc; que são controlados pelas opções de parâmetros na linha de comando.
  Os parâmetros utilizados podem ser visualizados no script rnaseq.sh.

Resultado do alinhamento:
Segue os resultados estatisticos do alinhamento, resumidamente:

**SAMPLEA1**
23902 aligned fragments; of these:
  23811 were paired; of these:
    3 aligned concordantly 0 times
    23800 aligned concordantly exactly 1 time
    8 aligned concordantly >1 times
Overall,  99.61% of aligned fragments aligned as proper pairs
  91 were single-end reads; of these:
      91 single-end reads mapped uniquely
      0 single-end reads mapped >1 times

SAMPLEA2

23857 aligned fragments; of these:
  23749 were paired; of these:
    0 aligned concordantly 0 times
    23746 aligned concordantly exactly 1 time
    3 aligned concordantly >1 times
Overall,  99.55% of aligned fragments aligned as proper pairs
  108 were single-end reads; of these:
      108 single-end reads mapped uniquely
      0 single-end reads mapped >1 times

SAMPLEB1

23865 aligned fragments; of these:
  23776 were paired; of these:
    1 aligned concordantly 0 times
    23769 aligned concordantly exactly 1 time
    6 aligned concordantly >1 times
Overall,  99.62% of aligned fragments aligned as proper pairs
  89 were single-end reads; of these:
      89 single-end reads mapped uniquely
      0 single-end reads mapped >1 times

SAMPLEB2

23804 aligned fragments; of these:
  23707 were paired; of these:
    0 aligned concordantly 0 times
    23702 aligned concordantly exactly 1 time
    5 aligned concordantly >1 times
Overall,  99.59% of aligned fragments aligned as proper pairs
  97 were single-end reads; of these:
      97 single-end reads mapped uniquely
      0 single-end reads mapped >1 times

Resultado:

Total de reads alinhadas: 95428

Total de reads alinhadas Paired End: 95043

Total de reads alinha Single-end: 385

Montagem dos transcritos com Cufflinks

Para quantificar a expressão gênica a partir de dados de RNA-Seq, é necessário, primeiramente, a identificação de quais isoformas de um dado gene foram produzidas por cada read. Para isso, é preciso identificar todas as isoformas daquele gene.
O software cufflinks realiza a montagem de transcritos individuais provenientes das reads dos dados de RNA-Seq, que apresentaram alinhamento contra o genoma de referência. Como uma amostra pode conter reads que são oriundas de variantes de splicing de um gene, o programa precisa inferir a estrutura de splicing de cada gene. Entretanto, os genes podem apresentar muitos eventos de splicing alternativo, resultando em diferentes reconstruções possiveis do gene. Sendo assim, o programa cufflinks realiza a montagem do transcriptoma utilizando a parcimonia. O algoritmo do programa relata o comprimento dos fragmentos dos transcritos (transfrags) com o objetivo de explicar os eventos de splicing nos dados (Trapnell, C. et al.2012).
Os outputs do programa são: genes.fpkm_tracking,isoforms.fpkm_tracking,skipped.gtf,transcripts.gtf para cada amostra.

Para testar outros montadores, também foi utilizado o programa StringTie, que utiliza um fluxograma como algortimo como passo de uma montagem opcional de novo, visando montar e quantificar transcritos completos, representando as variantes causadas por splicing de cada locus gênico. Como input também utiliza as reads processadas. Seu output foi utilizado pelo cuffdiff para realizar análise de expressão gênica diferencial.

Após a montagem dos transfrags de cada amostra, deve-se juntar os mesmo utilizando o programa Cuffmerge e Stringmerge, os quais podem utilizar uma anotação do genoma de referencia (formato gtf) na montagem final (Figura 9). Genes de baixa expressão pode apresentar baixa profundidade para sua reconstrução total, entretanto, ao se unir as montagens, pode-se recuperar o gene completo. O output é um unico arquivo de anotação no formato gtf que pode ser utilizado na analise de genes diferencialmente expressos.

Figura 9. Junção das montagens de cada amostra utilizando uma anotação de referência.
Fonte: (Trapnell, C. et al.2012)

A identificação de novos genes e transcritos pode ser dificil distinguir um transcrito novo completo de fragmentos parciais. Gaps do processo de sequenciamento pode causar falhas na reconstrução dos transcritos. Para ajudar, o programa cuffcompare realiza a comparação das montagens do cufflinks e arquivos de anotação do genoma de referência para ajudar a selecionar novos genes dos já anotados (Trapnell, C. et al.2012).

Cada uma das réplicas foram quantificadas utilizando o programa cuffquant. O output são arquivos de abundancia no formato cxb que são utilizados no programa cuffdiff para análise de expressão diferencial, e ao programa cuffnorm, visando a obtenção de tabelas de valores de expressão normalizados para que fiquem em uma mesma escala, nesse caso a normalização é feita pelo tamanho da biblioteca.

Análise de genes diferencialmente expressos - Cuffdiff

Em seguida, o nivel de expressão de cada transfag foi calculado utilizando o Cuffdiff, o qual calcula a expressão das amostras e testa a significancia estatistica de cada alteração observada entre elas.

O input utilizado no presente relatório foi o output do montador StringTie.

Para isso, é assumido que o número de reads produzidas por cada transcrito é proporcional a sua abundancia (Trapnell, C. et al.2012).
Quando se tem multpiplas réplicas, o programa aprende como a contagem de reads varia para cada gene e usa essa variância para calcular a significancia das alterações observadas na no perfil de expressão.
Como output do programa, tem-se vários arquivos contendo o resultado da análise de expressão genica diferencial. Em conjunto com os niveis de expressão dos transcritos e genes, tem-se estatisticas como fold change e p valor. o cuffdiff também reporta analise diferencial de splicing alternativo (genes), ou genes diferencialmente regulados por promotores.
Além disso, o programa agrupa isoformas de um gene que apresentem o mesmo TSS (sitio de inicio de transcrição), ou seja, isoformas derivadas de um mesmo pre-mRNA. Logo, alterações na abundancia de um reflete no splicing diferencial do seu pre-mRNA. É calculado também o nível de expressão do grupo TSS, adicionando os niveis de expressão das isoformas. Genes que apresentam mais de um TSS, é analisado a abundancia relativa entre eles, o que reflete em alterações da preferencia de TSS e promotores entre as condições analisadas. A figura 10 mostra como o Cuffdiff agrupa os TSS e os utiliza para inferir regulação gênica diferencial (Trapnell, C. et al.2012).

Figura 10. Inferência de regulação gênica diferencial a partir de agrupamentos de TSS. (a) Os genes produzem variantes a partir de splicing alternativo em diferentes quantidades por meio de diferentes TSSs. (b)As isoformas são então agrupadas de acordo com seu TSS e então é realizada uma análise para buscar abundancia relativa entre e dentro desses agrupamentos, que pode levar a pistas da diferença de regulação entre eles. © Nesse gene hipotético, diferentes abudancias entre as isoformas A e B dentro de um agrupamento TSS pode ser devido a splicing diferencial do transcrito primários. (d) A adição de seus níveis de expressão gera um
valor da expressão para este transcrição primário. (e) Alterações nesse nivel de expressão, por ex. isoforma C, indica possivel preferencia diferencial de promotores entre as condições. (f,g) Genes com multiplos CDS (f) também podem ser analisados para sequencias codificantes de proeina diferenciais (g).
Fonte: (Trapnell, C. et al.2012)

Como output tem-se também vários arquivos, entre eles tabelas que mostram a analise de expressão diferencial de genes e isoformas. O arquivo que contem a expressão diferencial de genes foi utilizada para verrificar a concordância dos genes diferencialmente expressos, com os genes que escolhemos no inicio para montar a simulação (Arquivos abundance A e B).

Comparando os resultados do cuffdiff com os arquivos abundance A e B:

Para realizar a comparação, teve-se que buscar qual gene_ID do cuffdiff é referente a cada gene dos arquivos abundance, e em seguida avaliar o log fold change ver se são concordantes com os valores inicias. Foi costruída uma tabela com os respectivos IDS e Genes do arquivo de abundancia e os valores dados a cada um no arquivo, e por ultimo o valor de expressão log fold change a titulo de comparação:

Resposta: Comparando os resultados log2 com os valores dos arquivos abundance A e B, pode-se observar que os resultados estão de acordo com o esperado. A relação está concordante com o esperado, pois por exemplo, o gene NM_125091.4, era esperado que ele estivesse mais expresso na condição A, e analisando o log2 fold change, é isso o que acontece, com um valor de -13.2891. O reverso acontece com o gene NM_001337167.1, ele se encontra mais expresso na condição B, com um log2 igual a 13.1845, novamente concordando com os arquivos Abundance A e B. O valor log 2 se dá pela razão entre value_2(referente a amostra 2) w o value_1 (referente a amostra 1). Observa-se também que quando colocamos valores similares em ambos os arquivos de abundancia, o log2 resulta em torno de zero, comprovando a concordancia entre o resultado e a simulação proposta.

Fusão de dados

Em seguida, os dados oriundos do programa cuffdiff, mais especificamente, o arquivo gene_exp.diff foi juntado com descrições dos genes obtidos no banco de dados do NCBI. A fusão foi realizada baseando-se em uma coluna comum, que apresente um identificador comum entre os dois arquivos.
Foi utilizado o comando abaixo para obtenção do arquivo contendo os genes de arabdopsis thaliana no banco de dados do NCBI:

esearch -db gene -query "3702[Taxonomy ID]" \
		| efetch -format tabular > ./Ref/gene_result.txt

A seguir, foram selecionadas somente as colunas de interesse a partir do arquivo gene_result.txt, originando o arquivo gene_info.txt. As colunas selecionadas foram: GeneID, Symbol e description.

cut -f 3,6,8 ./refs/gene_result.txt > ./refs/gene_info.txt

Ocasionalmente pode ocorrer a fusão de genes que possuem ids distintos durante a montagem pelo cufflinks/cuffmerge, durante o seu processamento. Então os mesmos foram separados utilizando o comando abaixo:

splitteR.R  --x="./output/cuffdiff/gene_exp.diff" \
            --col.x="gene" \
            --by.x="," \
            --out="./output/cuffdiff/gene_exp.diff.splitted.txt"	

Por fim, a fusão dos arquivos foi realizada utilizando o comando abaixo:

-all.x = todas as linhas do arquivo gene_exp.diff.spplitted.txt são mantidas no novo arquivo gene_exp.diff.spplitted.merged.txt.
x e y = arquivos que serão juntados.
by.x e by.y = colunas a serem juntadas.

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"

Como resultado obtivemos uma tabela contendo todas as colunas do arquivo gene_exp.diff.spplitted.txt, mais as do arquivo gene_info.txt, sendo que as colunas gene e symbol são agora uma só, pois são idênticas.

Montagem de novo

O mesmo conjunto de reads utilizado para a análise transcriptômica, já processadas foi utilizado para realizar a montagem de novo utilizando o programa Trinity.

Para a realização da montagem de novo foi utilizado o script rnaseq-denovo.sh.
Para acessar o script: rnaseq-denovo

Para execução do script foi utilizada a linha de comando abaixo:
./rnaseq-denovo.sh ./output/processed/prinseq ./output

Nesse script inicialmente foi realizado um processo de renomeação das reads para que contenham o sufixo /1 e /2, para as reads 1 e 2 respectivamente, o que é uma exigência do software trinity.

O montador Trinity também apresenta uma segunda opção, de realizar a montagem guiada por genoma de referência. Para isso, são utilizados os alinhamentos das reads contra o genoma referência, arquivos no formato BAM, outputs do alinhador STAR. Para essa montagem foi utilizado o script
rnaseq-ref-trinity.sh. Para acessá-lo:rnaseq-ref-trinity

O script foi executado com a seguinte linha de comando:
./rnaseq-ref-trinity.sh ./output/star_out_final ./output

Primeiramente é realizado uma etapa de junção de todos os alinhamentos obtidos no programa STAR em um só arquivo: All.sorted.bam.

Os parametros de ambas as montagens estão comentados nos respectivos scripts.

Avaliação das montagens

Primeiramente foi realizada a indexação do transcriptoma utilizando o seguinte comando:
Indexação do transcriptoma utilizado para a simulação:

makeblastdb -title "RefTrans" \
			-dbtype nucl \
			-out RefTrans \
			-in transcriptoma.fa \
			 > makebleastdb.log.out.txt \
			2> makeblastdb.log.err.txt

onde:
-title: Titulo do banco de dados que escolher;
-dbtype: tipo do banco de dados, aminoacidos ou nucleotideos.
-out:output;
-in: Arquivo de entrada, no caso o transcriptoma montado.

Como a indexação foi bem sucedida, apareceram os seguintes arquivos:

• RefTrans.nhr
• RefTrans.nin
• RefTrans.nsq

Em seguida foi realizado a comparação entre os arquivos utilizando blastn:

• Utilizando a montagem com Trinity de novo:

blastn 	-max_hsps 1 \
		-max_target_seqs 1 \
		-num_threads 8 \
		-query ./output/trinity_assembled/Trinity.fasta \
		-task blastn \
		-db ./RefTrans \
		-out ./Trinity_x_RefTrans.txt \
		-evalue 1e-5 \
		-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
		 > ./Trinity_x_RefTrans.log.out.txt \
		2> ./Trinity_x_RefTrans.log.err.txt

-max_hsps: Número maximo de HSPs(Segmentos de alto score de alinhamento) para manter para qualque query single - subject pair.
-max_target_seqs: Número de sequências alinhadas a se manter
-num_threads: número de processadores a serem utilizados paralelamente
-query: Sequencia de entrada
-task: Qual blast será realizado
-db: Nome do banco de dados
-out: arquivo de saída
-evalue: numero esperado de alinhamentos melhores ou iguais ao que é encontrado.
-outfmt: Opções de visualização do alinhamento. O 6 indica formato tabular. qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore > default, equivalente a std.

Para visualizar somente os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos foi utilizado o comando abaixo:

cut -f 1,2 ./Trinity_x_RefTrans.txt

Resposta: total de HITs: 31

Em seguida foram contadas quantas ocorrências de sequências da base de dados (subject) foram encontradas:

cut -f 2 ./Trinity_x_RefTrans.txt | sort | uniq -c

**Resposta: Foram encontradas 6 ocorrências de sequencias na base dados, concordando com os genes escolhidos no inicio para criação do arquivo abundance A e B: **

5 NM_001337167.1
7 NM_001345288.1
6 NM_105102.3
1 NM_123951.3
6 NM_125091.4
6 NM_180129.4

Comparando com o arquivo abundanceA.txt:

cat /home/jcardoso/sim/abundanceA.txt
NM_123951.3     0.15
NM_001345288.1  0.25
NM_105102.3     0.05
NM_180129.4     0.05
NM_125091.4     0.50
NM_001337167.1  0

Resposta: Pode-se observar que os genes montados são concordantes com os genes que se encontram nos arquivos abundance A e B.

• Para a montagem do Trinity que utilizou o Genome Guided:

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

Similarmente aos passos anteriores:

cut -f 1,2 ./Trinity-GG_x_RefTrans.txt
GG1|c0_g1_i1 NM_001337167.1
GG2|c0_g1_i1 NM_180129.4
GG3|c0_g1_i1 NM_001345288.1
GG4|c0_g1_i1 NM_125091.4
GG5|c0_g1_i1 NM_123951.3
GG6|c0_g1_i1 NM_105102.3

cut -f 2 ./Trinity-GG_x_RefTrans.txt | sort | uniq -c
1 NM_001337167.1
1 NM_001345288.1
1 NM_105102.3
1 NM_123951.3
1 NM_125091.4
1 NM_180129.4

Novamente os genes são concordantes com o arquivo criado incialmente abundance A e B.

Análise de expressão gênica diferencial referente a montagem de novo:

Com o objetivo de obter os valores de Log Fold-Change e os p-valores ajustados referentes às comparações pareadas entre SAMPLEA e SAMPLEB, foi utilizado o um script auxiliar do montador Trinity:

Primeiramente o caminho do diretório base do Trinity, ou seja, onde se encontra o programa e todas as suas dependencias foi buscado pr meio do comando abaixo:

echo ${TRINITY_HOME}

O caminho é : /usr/local/bioinfo/trinityrnaseq

Em seguida foi verificado os parâmetros referentes aos scripts auxiliares align_and_estimate_abundance.pl e abundance_estimates_to_matrix.pl:

${TRINITY_HOME}/util/align_and_estimate_abundance.pl

${TRINITY_HOME}/util/abundance_estimates_to_matrix.pl

Para realizar a análise de expressão gênica diferencial entre as condições A e B foi utilizado o script run-DESeq2.R

Para executar esse protocolo de análise de expressão gênica diferencial de uma forma automatizada, foi criado o script rnaseq-trinity-abundance.sh. Esse script está disponível aqui.

Para a análise foi utilizado comando de execução abaixo:

rnaseq-trinity-abundance.sh	./output/renamed/ \
				./output/trinity_assembled/ \
				./output

Para visualizar o resultado da anális com DESeq2 de modo alinhado:

cat ./output/abundance/DEG/SAMPLEA-SAMPLEB.txt | column -s$'\t' -t | less -S

A tabela abaixo foi aberta por meio do programa RStudio:

Os genes diferencialmente expressos novamente estão concordantes com o arquivo de abundancia gerado no inicio. Os genes TRINITY_DN1_C0_g1, TRINITY_DN5_c0_g1 e TRINITY_DN0_c0_g1 se encontram diferencilamente expressos, pois apresentam log Fold Change maior que 2/-2. O restante a diferença de expressão não é estatisticamente significativa.

R

Foi realizado uma etapa de visualição dos resultados no programa R. Primeiramente os programas as bibliotecas a serem utilizadas foram carregadas:

library("xlsx")

library("gplots")

library("RColorBrewer")

• Para visuzalizar o caminho de onde vc está: (Semelhante ao pwd no linux)

getwd()

/home/jcardoso

• Para entrar no diretório em que estão contidas as análises de DEG:

> setwd("/state/partition1/jcardoso/Relatorio")

> getwd()

"/tank/partition1/jcardoso/Relatorio"

• Em seguida o arquivo texto gerado na análise de DEG foi carregado, exigindo separação por TAB, com cabeçalho e sem transformar strings em fatores:

deg.df <- read.delim(file="./output/abundance/DEG/SAMPLEA-SAMPLEB.txt", 
                     sep="\t", 
                     header=TRUE, 
                     stringsAsFactors = FALSE)

Para ver a dimensão do data.frame, onde o primeiro valor refere-se á quantidade de linhas e o segunda a quantidade de colunas:

dim(deg.df)
[1] 6 8

• Em seguida foi selecionado um subconjunto de linhas desse data frama em outro data.frama denominado subset.deg.df, contendo somente os genes que apresentarem valor de logFC >= 2 ou logFC <= -2 e p-valor corrigido ("FDR") <= 0.05:

subset.deg.df <- subset(deg.df, (((logFC >= 2) | (logFC <= -2) ) & (FDR <=0.05)) )

• Avaliando a dimensão desse subconjunto gerado:
  dim(subset.deg.df)
[1] 3 8

• O próximo passo foi a criação de um vetor de nomes de colunas selecionadas, essas colunas serão ordenadas depois de secionados utilizando 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") 
                      )
                    )

• Foi gerada uma nova matriz denominada "expression_data" que contém somente as colunas selecionadas:

expression_data <- as.matrix( subset.deg.df[,sel_columns] )

• Foram atribuidos o conteúdo da coluna "X" para os nomes das linhas da matriz (Identificador do gene):

rownames(expression_data) <- subset.deg.df$X

head(expression_data,2)
                  SAMPLEA1 SAMPLEA2 SAMPLEB1 SAMPLEB2
TRINITY_DN0_c0_g1 869.9345 1588.315 6013.741 7706.994
TRINITY_DN1_c0_g1   0.0000    0.000 4014.274 3140.861

• 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"
           )

• Gerando uma paleta personalizada de 299 cores do vermelho ao verde, passando pelo amarelo:

my_palette <- colorRampPalette(c("red", "purple", "blue"))(n = 299)

• As quebras das cores podem ser definidas manualmente (transição das cores):

col_breaks = c(seq(-1,0,length=100),        # for red
               seq(0.01,0.8,length=100),    # for purple
               seq(0.81,1,length=100))      # for blue

• Visando a geração de uma imagem de tamanho 5x5 polegadas, foi utilizado o comando abaixo:

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 a geração do heatmap:

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
          )

No heatmap gerado pode ser observado que o gene TRINITY_DN0_c0_g1 se encontra expresso igualmente em ambas as condições e réplicas.
Enquanto que o gene TRINITY_DN1_c0_g1 está sendo mais expresso na condição B do que na condição A.
Por último, o gene TRINITY_DN5_c0_g1 teve uma leve diferença entre as réplicas da condição B. Mais ele se encontra mais expresso na condição A do que na condição B.

Conclusão

A tecnologia de RNA-Seq tem provado ser uma ferramenta valiosa para realizar analises transcriptomicas, pois fornece novos *insights * sobre a expressão gênica, análise de expressão diferencial, eventos de splicing alternativo, alterações na expressão na progressão de doenças, edição de RNA, expressão diferencial de TSS, isoformas e regulação diferencial (preferencia por promotores dada uma condição).

No presente relatório, foi realizada uma análise transcriptômica, entretanto de apenas uma pequena simulação de dados, genes extraidos do banco de dados NCBI, que foram fragmentados, simulando dados do RNA-Seq. Entretanto, mesmo sendo apenas uma simulação, pode-se observar uma concordância entre os arquivos de abundancia gerados no inicio, com os resultados finais de expressão gênica diferencial. Sendo que as proporção escolhidas são equivalentes as diferenças da expressão.
Também foi observado concordância na montagem de novo, e na montagem guiada por um genoma do montador Trinity. Em ambas as situações os genes escolhidos foram identificados e montados.

A análise completa contribuiu para agregar conhecimento de cada etapa de uma análise de transcriptoma, de quais softwares utilizar para cada caso, e por que. Além do inicio de um aprendizado em montar scritps, em ferramentas da bioinformática e mesmo do sistema Linux. Contribui também para entender sobre quais parametros utilizar em cada situação e as possibilidades de análises que se pode submeter os dados de RNA-Seq. Concluindo a disciplina, principalmente as aulas práticas foram de extrema importância para minha formação, e agregou conhecimentos básicos úteis para análises de dados de sequenciamento.

Referências Bibliográficas:

Anamika K, Verma S, Jere A, Desai A. Transcriptomic Profiling Using Next Generation Sequencing - Advances, Advantages, and Challenges. Next Generation Sequencing - Advances, Applications and Challenges. 2016;9:7355–7365.

DELATORRE, Carla Andréa e SILVA, Adriano Alves da. Arabidopsis thaliana: uma pequena planta um grande papel. Rev. de Ciências Agrárias. 2008, vol.31, n.2 [citado 2019-11-23], pp.58-67.

Didion et al., Atropos: specific, sensitive, and speedy trimming of sequencing reads.2017. PeerJ 5:e3720.

Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR: STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013, 29: 15-21. 10.1093/bioinformatics/bts635.

Grabherr, M.G. et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat. Biotechnol.2011. 29, 644–652.

Kumar KR, Cowley MJ, Davis RL. Next generation sequencing and emerging technologies. Semin Thromb Hemost 2019; 45 (07) 661-673.

Levy S.E., Myers R.M. Advancements in next-generation sequencing. Annu. Rev. Genomics Hum. Genet. 2016; 17:95–115.

Lowe R, Shirley N, Bleackley M, Dolan S, Shafee T. Transcriptomics technologies. PLOS Comput Biol 2017;13(5).

PrinSeq.PrinSeq Manual. Disponível em: http://prinseq.sourceforge.net/manual.html. Acesso em 24 de nov de 2019.

Trapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 2012.7, 562–578.

Voshall A, Moriyama EN: Next-generation transcriptome assembly: strategies and performance analysis. In Bioinformatics in the Era of Post Genomics and Big Data Edited by Abdurakhmonov I: IntechOpen; 2018.