Relatório Final: Análise Transcriptômica de Chlamydomonas reinhardtii

# Relatório Final: Análise Transcriptômica de *Chlamydomonas reinhardtii* > ### Relatório apresentado à disciplina de Bioinformática II (PPG-Micro) como requisito para obtenção de nota > URL curta: http://bit.ly/c-reinhardtii-transcriptome ## Sumário > [TOC] ## 1. Seleção do Organismo e Preparação do Ambiente ### 1.1. *Chlamydomonas reinhardtii* Trata-se de uma microalga eucariótica unicelular e flagelada, pertencente ao reino Plantae. Ela é um organismo modelo e se destaca por sua simplicidade estrutural, com uma única célula contendo um núcleo, cloroplastos, mitocôndrias e flagelos, entre outra organelas, como visto na figura abaixo, o que facilita a compreensão de processos biológicos fundamentais. Além disso, essa microalga possui um genoma sequenciado e anotado, o que facilita a análise dos dados de expressão gênica. <center> <img src="https://hackmd.io/_uploads/SJjLx3mP6.jpg" width="280"/> </center> >Fonte: https://doi.org/10.7554/eLife.39233.002 A escolha dessa espécie também está relacionada ao tipo de organismo que está sendo estudado no meu projeto de mestrado, que visa compreender as interações na rizosfera entre plantas de alface, microrganismos associados e microalgas verdes. ### 1.2. Obtenção das Informações Genômicas do Banco de Dados NCBI Ao acessar a página do Taxonomy Browser do NCBI, foi obtido o Taxonomy ID para essa espécie: 3055. Assim, as posteriores consultas no banco de dados foram feitas usando txid3055[Organism:exp]. <center> <img src="https://hackmd.io/_uploads/SyiaA37va.png"/> </center> No caso de transcritos no banco de dados de referência RefSeq os identificadores (IDs) começam com os prefixos "NM_" ou "NR_" para indicar, respectivamente, transcritos codificadores de proteínas (mRNAs) e transcritos não codificadores de proteínas (ncRNAs), ou ainda, podem ser "XM_" ou "XR_" quando as entradas possuem apenas uma previsão computacional, sem terem passado pela fase de curadoria, como os que foram usados no presente projeto. Por orientação, foram selecionados 6 genes do mesmo cromossomo Chr 1 NC_057004.1: XM_001689813.2 CHLRE_01g010400v5 XM_043058287.1 CHLRE_01g010450v5 XM_043058288.1 CHLRE_01g010500v5 XM_043058294.1 CHLRE_01g010800v5 XM_043058296.1 CHLRE_01g010832v5 XM_043058297.1 CHLRE_01g010848v5 ### 1.3. Estruturação dos Diretórios A estruturação de diretórios com os arquivos de referência, dados brutos, entrada de dados de sequênciamento e saída dos dados após o processamento, é necessária antes de começar as análises. Para criar os diretórios usados nesse projeto, usou-se o comando `makedir` ```bash= mkdir /data/bioaat/cdnazareno/C_reinhardtii # Criação do diretório raiz mkdir /data/bioaat/cdnazareno/C_reinhardtii/ref # Criação do diretório para os arquivos de referência mkdir /data/bioaat/cdnazareno/C_reinhardtii/raw # Criação do diretório para os dados brutos de sequenciamento mkdir /data/bioaat/cdnazareno/C_reinhardtii/input # Criação do diretório para os dados de sequenciamento organizados mkdir /data/bioaat/cdnazareno/C_reinhardtii/output # Criação do diretório para os dados pós-processados ``` * `/C_reinhardtii` refere-se ao diretório raiz de onde os demais partirão; * `/C_reinhardtii/ref` refere-se ao diretório que receberá os dados de referência, e que depois passarão pelo processo de minimização das informações genômicas, ou seja, a criação do Toy Genome; * `/C_reinhardtii/raw` refere-se ao diretório que receberá os dados brutos de sequenciamento; * `/C_reinhardtii/input` refere-se ao diretório que dará entrada dos dados de sequenciamento, já devidamente estruturados, para o processamento; * `/C_reinhardtii/output` refere-se ao diretório que recebeu os dados após todo o processamento. Para visualizar a árvore de diretórios criada, usou-se o comando `tree` <center> <img src="https://hackmd.io/_uploads/SyscIpAua.png" width="400"/> </center> ## 2. Entrada dos Dados ### 2.1. Diretório `/ref` Nesse diretório foram armazenados os dados genômicos de referência, para isso foi preciso obter os links para os arquivos do banco de dados NCBI (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/595/GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5/) e utilizar o protocolo FTP. Para fazer o download dos arquivos de referência, usou-se o comando `wget` ```bash= cd ./ref # Mudança de diretório, a partir do diretório raiz, para executar a função wget wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/595/GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5/GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5_genomic.fna.gz wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/595/GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5/GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5_genomic.gff.gz ``` Para descompactar os arquivos baixados, usou-se o comando `gunzip` e depois `tree` ```bash= gunzip *.gz ``` <center> <img src="https://hackmd.io/_uploads/rklhHqTAua.png" width=""/> </center> ### 2.2. Criação do Toy Genome A criação do Toy Genome foi dividida em três etapas. Na primeira, foi feita a limpeza do arquivo *.FNA (FASTA DNA) e, a partir dele, criado o arquivo *.FA (FASTA). Na segunda, foi utilizado um script para ajustar o arquivo *.GFF (Gene Finding Format). Por fim, foi criado o arquivo *.GTF (Gene Transfer Format) a partir da transformação do arquivo *.GFF. Para a etapa de limpeza do arquivo *.FNA, criou-se o script `cleanfasta.sh` com o comando `nano` ```bash= #!/bin/bash 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 sed 's/^>$[^\.]\+$.*/>\1/' ${infile} ``` Com o comando `chmod a+x` o script recebeu permissão para ser executado e o operador `>` redirecionou a saída para o arquivo *.FA ```bash= chmod a+x cleanfasta.sh ./cleanfasta.sh GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5_genomic.fna > genoma.fa ``` Para ajustar o aquivo *.GFF, usou-se o script `fixNCBIgff.sh` disponível no GitHub (https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixNCBIgff.sh) e em seguida isolou-se apenas a parte referente ao cromossomo 1 com o comando `grep` ```bash= fixNCBIgff.sh GCF_000002595.2_Chlamydomonas_reinhardtii_v5.5_genomic.gff > genoma.gff grep -P '^(##gff|NC_057004.1\t)' genoma.gff > genome.gff ``` Para criar um transcriptoma e um proteoma de referência com base no arquivo *.GFF, usou-se `gffread` com os parâmetros a seguir ```bassh= gffread genome.gff -g genome.fa -w transcriptome.fa gffread genome.gff -g genome.fa -y proteome.fa ``` Por fim, para transformar um arquivo *.GFF em *.GTF, usou-se o mesmo comando, mas com outro parâmetro ```bassh= gffread genome.gff -g genome.fa -T -o genome.gtf ``` Os arquivos desnecessários foram excluídos com o comando `rm` <center> <img src="https://hackmd.io/_uploads/rktsXXyYT.png" width=""/> </center> ### 2.3. Diretório `/raw` Nesse diretório foram armazenados os dados brutos de sequenciamento, obtidos por meio de um script, mas antes disso foi preciso criar um arquivo *.TXT contendo os IDs dos 6 transcritos selecionados na primeira coluna e seus respectivos nomes na segunda. Para mudar de diretório a partir do `/ref`, usou-se `cd ../raw`, e depois, para criar o arquivo de texto, usou-se `nano accs.txt` ```bash= XM_001689813.2 CHLRE_01g010400v5 XM_043058287.1 CHLRE_01g010450v5 XM_043058288.1 CHLRE_01g010500v5 XM_043058294.1 CHLRE_01g010800v5 XM_043058296.1 CHLRE_01g010832v5 XM_043058297.1 CHLRE_01g010848v5 ``` Para criar o script que coletou as sequências fasta dos genes, usou-se o comando `nano collectfasta.sh` ```bash= !#/bin/bash for acc in $(cut -f 1 ACCS.txt); do echo "Pegando FASTA para ${acc} ..." esearch -db nucleotide -query ${acc} | efetch \ -format fasta >> transcriptome.fa done ``` Com o comando `chmod a+x` o script recebeu permissão para ser executado ```bash= chmod a+x collectfasta.sh ./collectfasta.sh ``` <center> <img src="https://hackmd.io/_uploads/B1znimkKp.png" width=""/> </center> Na sequência foram criados dois arquivos de abundância para simular duas condições biológicas A e B. Para a condição A, usou-se `nano abundance_A.txt` ```bash= XM_001689813.2 0.50 XM_043058287.1 0.05 XM_043058288.1 0.05 XM_043058294.1 0.25 XM_043058296.1 0.15 XM_043058297.1 0 ``` Para a condição B, usou-se `nano abundance_B.txt` ```bash= XM_001689813.2 0.05 XM_043058287.1 0.05 XM_043058288.1 0.05 XM_043058294.1 0.25 XM_043058296.1 0.10 XM_043058297.1 0.50 ``` Por fim, foram simuladas 1000 sequências para cada tratamento em três repetições e leituras paired-end de tamanho 300pb. Para realizar a simulação, usou-se o script `simLib.pl` ```bash= conda activate python2 # Antes de executar um script em python é preciso ativar o ambiente virtual do python simLib.pl -a abundance_A.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA1 simLib.pl -a abundance_A.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA2 simLib.pl -a abundance_A.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA3 simLib.pl -a abundance_B.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB1 simLib.pl -a abundance_B.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB2 simLib.pl -a abundance_B.txt -i transcriptome.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB3 ``` <center> <img src="https://hackmd.io/_uploads/SynXXN1Y6.png" width=""/> </center> ### 2.4. Diretório `/input` Nesse diretório foram armazenados links simbólicos dos arquivos gerados no diretório anterior, e que posteriormente serão processados, mas antes disso foi preciso compactar e designar novos nomes para esses arquivos. Para simular os dados da mesma forma que receberíamos de uma facility, usou-se o comando para compactação `gzip` ```bash= gzip *.fastq ``` Para renomear os arquivos para TREAT (tratamento) e CNTL (controle) de forma automatizada, usou-se o comando `nano data.txt` para criar um arquivo *.TXT com os novos nomes correspondentes ```bash= SAMPLEA1_R1.fastq.gz CNTRL_B1_R1.fastq.gz SAMPLEA1_R2.fastq.gz CNTRL_B1_R2.fastq.gz SAMPLEA2_R1.fastq.gz CNTRL_B2_R1.fastq.gz SAMPLEA2_R2.fastq.gz CNTRL_B2_R2.fastq.gz SAMPLEA3_R1.fastq.gz CNTRL_B3_R1.fastq.gz SAMPLEA3_R2.fastq.gz CNTRL_B3_R2.fastq.gz SAMPLEB1_R1.fastq.gz TREAT_B1_R1.fastq.gz SAMPLEB1_R2.fastq.gz TREAT_B1_R2.fastq.gz SAMPLEB2_R1.fastq.gz TREAT_B2_R1.fastq.gz SAMPLEB2_R2.fastq.gz TREAT_B2_R2.fastq.gz SAMPLEB3_R1.fastq.gz TREAT_B3_R1.fastq.gz SAMPLEB3_R2.fastq.gz TREAT_B3_R2.fastq.gz ``` Em seguida, para gerar seus links simbólicos, usou-se o script `createSymLink.sh`, executado a partir do diretório raiz ```bash= createSymLink.sh ./input/data.txt ./raw/ ./input/ ``` <center> <img src="https://hackmd.io/_uploads/BJNIySJK6.png" width=""/> </center> ## 3. Pré-processamento dos Dados Todas as etapas do pré-processamento foram automatizadas pelo script fornecido `preprocess.sh`. Os programas executados com o script foram o FastQC, Atropos e Prinseq. ### 3.1. FastQC O FastQC é uma ferramenta usada para avaliar a qualidade de dados de sequenciamento de DNA. Ele fornece métricas detalhadas, como distribuição de qualidade por base e identificação de possíveis problemas, auxiliando na avaliação inicial de dados de sequenciamento. Antes de executar o script foi analisada a qualidade dos dados por meio da interface gráfica de usuário (GUI) do FastQC com o comando `fastqc` <center> <img src="https://hackmd.io/_uploads/Skr7xOyKp.png" width=""/> </center> > Gráfico da qualidade das leituras R1 (foward) para SAMPLEA1 <center> <img src="https://hackmd.io/_uploads/rySPxOkYT.png" width=""/> </center> > Gráfico da qualidade das leituras R2 (reverse) para SAMPLEA1 Abaixo dois trechos do script que avaliam a qualidade das reads R1 e R2 por meio do FastQC ```bash= fastqc -t ${THREADS} \ # Comando para executar a análise de qualidade usando FastQC, com -t especificando o número de threads (ou processadores) a serem utilizados ${r1} \ # Caminho para o arquivo de leitura R1 -o ${outdir}/processed/fastqc/pre/ \ # Diretório de saída para os resultados do FastQC 1> ${outdir}/processed/fastqc/pre/${name}_R1.log.out.txt \ # Redireciona a saída padrão (stdout) para esse arquivo 2> ${outdir}/processed/fastqc/pre/${name}_R1.log.err.txt # Redireciona a saída de erro (stderr) para esse arquivo ``` ```bash= fastqc -t ${THREADS} \ ${r2} \ -o ${outdir}/processed/fastqc/pre/ \ 1> ${outdir}/processed/fastqc/pre/${name}_R2.log.out.txt \ 2> ${outdir}/processed/fastqc/pre/${name}_R2.log.err.txt ``` ### 3.2. Atropos O Atropos é uma ferramenta de pré-processamento de dados de sequenciamento de alta qualidade. Ele é usado para realizar tarefas como remoção de adaptadores, qualidade de filtragem e trimming, contribuindo para preparar dados de sequenciamento para análises subsequentes. Logo após a execução do FastQC, o Atropos é usado para a trimagem (poda) de sequências de baixas qualidade, ou com o Phred Quality Score abaixo de 20 (faixa amarela no FastQC) ```bash= atropos trim --aligner insert \ # Inicia o comando Atropos para a remoção de adaptadores usando o algoritmo de alinhamento 'insert' --error-rate 0.1 \ # Define a taxa máxima de erro permitida durante o alinhamento de adaptadores --times 2 \ # Especifica o número máximo de vezes que a remoção de adaptadores será tentada --minimum-length 1 \ # Define o comprimento mínimo permitido para uma read após o pré-processamento --op-order GAWCQ \ # Especifica a ordem de operações para as etapas de qualidade e adaptação --match-read-wildcards \ # Habilita o uso de curingas (wildcards) nas sequências dos adaptadores --overlap 3 \ # Define o número de bases que as reads devem se sobrepor com os adaptadores para serem consideradas uma correspondência --quality-cutoff 10 \ # Define o limite de qualidade para cortar bases no final de uma read --threads ${THREADS} \ # Especifica o número de threads (ou processadores) a serem usados durante o processo --correct-mismatches conservative \ # Realiza a correção conservadora de mismatches durante o alinhamento de adaptadores --pair-filter any \ # Define o filtro para pares de reads, permitindo qualquer configuração de orientação e estado de trim -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTGAAAAA \ # Especifica as sequências dos adaptadores para as leituras R1 -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTNNNNNNNNNNNNNNNNNNNNGTGGTCGCCGTATCATTAAAAAA \ # Especifica as sequências dos adaptadores para as leituras R2 -o ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \ # Define os caminhos de saída para as reads R1 após o pré-processamento -p ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \ # Define os caminhos de saída para as reads R2 após o pré-processamento -pe1 ${r1} \ # Especifica os arquivos de entrada para as reads R1 -pe2 ${r2} \ # Especifica os arquivos de entrada para as reads R2 --untrimmed-output ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq \ # Define os caminhos de saída para reads que não têm adaptadores removidos --untrimmed-paired-output ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq \ > ${outdir}/processed/atropos/${name}.atropos.log.out.txt \ # Redireciona a saída padrão (stdout) para um arquivo de log de saída 2> ${outdir}/processed/atropos/${name}.atropos.log.err.txt # Redireciona a saída de erro (stderr) para um arquivo de log de erro ``` ```bash= atropos trim --aligner adapter \ # Inicia o comando Atropos para a remoção de adaptadores usando o algoritmo de alinhamento 'adapter' -e 0.1 \ # Define a taxa máxima de erro permitida durante o alinhamento de adaptadores -n 2 \ # Especifica o número máximo de vezes que a remoção de adaptadores será tentada -m 1 \ # Define o número mínimo de bases sobrepostas (match) necessárias para considerar a remoção de adaptadores --match-read-wildcards \ # Habilita o uso de curingas (wildcards) nas sequências dos adaptadores -O 3 \ # Especifica a sobreposição mínima necessária para a remoção de adaptadores -q 20 \ # Define o limite de qualidade para cortar bases no final de uma read --pair-filter any \ # Define o filtro para pares de reads, permitindo qualquer configuração de orientação e estado de trim -pe1 ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq \ # Especifica os arquivos de entrada para as reads R1 que não foram previamente processadas (untrimmed) -pe2 ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq \ # Especifica os arquivos de entrada para as reads R2 que não foram previamente processadas (untrimmed) -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTGAAAAA \ # Especifica as sequências dos adaptadores para as leituras R1, permitindo curingas nas sequências -g TTTTTTAATGATACGGCGACCACNNNNNNNNNNNNNNNNNNNNACACTCTTTCCCTACACGACGCTCTTCCGATCT \ # Especifica as sequências dos adaptadores para as leituras R2, permitindo curingas nas sequências -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTNNNNNNNNNNNNNNNNNNNNGTGGTCGCCGTATCATTAAAAAA \ -G TTTTTCAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT \ -T ${THREADS} \ -o ${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \ -p ${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \ > ${outdir}/processed/atropos/${name}.atropos_adapter.log.out.txt \ 2> ${outdir}/processed/atropos/${name}.atropos_adapter.log.err.txt ``` Durante a execução do Atropos três arquivos foram gerados para cada leitura (com os sufixos insert, adapter e untrimmed). Por fim, o insert e adapter são concatenados em um só arquivo denominado final. ### 3.3. Prinseq O Prinseq é uma ferramenta versátil para a filtragem e análise de dados de sequenciamento. Ele é utilizado para remover sequências de baixa qualidade, duplicatas e contaminantes, além de fornecer estatísticas úteis sobre a composição e qualidade do conjunto de dados. Devido ao enriquecimento feito para a captura de mRNA na fase de montagem da biblioteca para o sequenciamento, as reads acabam apresentando a cauda poli-A. Para lidar com essa situação o programa Prinseq foi inserido no script ```bash= prinseq-lite.pl -fastq ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq \ # Executa o Prinseq nos arquivos de leitura R1 provenientes do resultado do processo Atropos -fastq2 ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq \ # Executa o Prinseq nos arquivos de leitura R2 provenientes do resultado do processo Atropos -out_format 3 \ # Define o formato de saída como formato 3, que geralmente é um formato tabular -trim_qual_window 3 \ # Janela para análise da qualidade -trim_qual_step 1 \ # Passo usado na análise da qualidade -trim_qual_type mean \ # Utiliza a média para determinar a qualidade -trim_qual_rule lt \ # Regra para determinar se a qualidade é menor que o valor determinado -trim_qual_right 30 \ # Valor de corte para a qualidade à direita -out_good ${outdir}/processed/prinseq/${name}.atropos_final.prinseq \ # Define o arquivo de saída para as sequências filtradas e de boa qualidade -out_bad null \ # Descarta as sequências de baixa qualidade -lc_method dust \ # Método para remoção de complexidade baixa -lc_threshold 50 \ # Limiar para a remoção de complexidade baixa -min_len 20 \ # Comprimento mínimo permitido para as sequências após a filtragem -trim_tail_right 5 \ # Remove 5 bases do final de cada sequência -trim_tail_left 5 \ # Remove 5 bases do início de cada sequência -ns_max_p 80 \ # Porcentagem máxima permitida de bases ambiguas (N) em uma sequência -min_qual_mean 28 \ # Define a média mínima de qualidade permitida para uma sequência > ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.out.log \ # Redireciona a saída padrão (stdout) para um arquivo de log de saída 2> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.err.log # Redireciona a saída de erro (stderr) para um arquivo de log de erro ``` Após a execução do Prinseq, retornamos à GUI do FastQC para analisar a qualidade das nossas reads <center> <img src="https://hackmd.io/_uploads/rkWENd1Y6.png" width=""/> </center> > Gráfico da qualidade das leituras R1 (foward) para SAMPLEA1 <center> <img src="https://hackmd.io/_uploads/B1_4Nd1Kp.png" width=""/> </center> > Gráfico da qualidade das leituras R2 (reverse) para SAMPLEA1 Ao final da execução do script o comando `tree -d` foi usado para verificar os diretórios criados dentro de `/output` <center> <img src="https://hackmd.io/_uploads/SJuy7uyF6.png" width=""/> </center> ## 4. Alinhamento das Leituras O alinhamento de leituras é necessário em projetos de sequenciamento para mapear as sequências ao genoma de referência, identificar variantes genéticas e quantificar a expressão gênica. Ele fornece uma base crítica para análises subsequentes, contribuindo para a compreensão abrangente de informações biológicas nos dados de sequenciamento. No presente projeto, toda a etapa de alinhamento, montagem, comparação com o transcriptoma de referência e análise de expressão diferencial foi automatizada pelo script fornecido. Os alinhadores presentes nesse script são o STAR e o TopHat. Em primeiro momento, usou-se o alinhador STAR. Antes de executar o alinhador, é preciso otimizar o código com os resultados obtidos no script `checkMinMaxIntronSize.sh` ```bash= checkMinMaxIntronSize.sh ./ref/genome.gff ./output/ 0.9 ``` Depois de obtido um valor mínimo e máximo de íntrons conforme o arquivo *.GFF, otimizou-se o seguinte trecho do script `rnaseq-ref.sh` ```bash= STAR --runThreadN ${THREADS} \ --genomeDir ${outdir}/star_index/ \ # Indica o diretório onde o índice do genoma de referência para o STAR está localizado --readFilesIn ${r1} ${r2} \ # Especifica os arquivos de leitura (reads) R1 e R2 para o alinhamento --outSAMstrandField intronMotif \ # Inclui informações sobre os motivos dos íntrons no formato SAM --outFilterIntronMotifs RemoveNoncanonical \ # Remove íntrons não canônicos durante o alinhamento --sjdbGTFfile ${refgtf} \ # Fornece um arquivo GTF contendo informações de anotação genômica --outFilterMultimapNmax 20 \ # Define o número máximo de locais de mapeamento para reads multimapeadas --outFileNamePrefix ${outdir}/star_out_pe/${name}/ \ # Define o prefixo para os arquivos de saída gerados pelo STAR --outSAMtype BAM Unsorted \ # Gera saída no formato BAM e não ordena o arquivo de saída --outFilterType BySJout \ # Filtra reads com junções não canônicas --outSJfilterReads Unique \ # Filtra reads com junções exclusivas --alignSJoverhangMin 8 \ # Especificam os requisitos mínimos para a sobreposição de reads nas junções --alignSJDBoverhangMin 1 \ --outFilterMismatchNmax 999 \ # Define o número máximo permitido de mismatchs por read --outFilterMismatchNoverReadLmax 0.04 \ # Define a proporção máxima permitida de mismatchs em relação ao comprimento da read --alignIntronMin 10 \ # Especifica o tamanho mínimo de íntrons permitidos conforme execução do script anterior --alignIntronMax 400000 \ # Especifica o tamanho máximo de íntrons permitidos conforme execução do script anterior --alignMatesGapMax 200 \ # Define o tamanho máximo do gap entre reads emparelhadas > ${outdir}/star_out_pe/${name}.log.out.txt \ 2> ${outdir}/star_out_pe/${name}.log.err.txt ``` E, em seguida, executou-se o script da seguinte forma ```bash= rnaseq-ref.sh star ./input ./output 3 ./ref/genome.gtf ./ref/genome.fa stringtie ``` Após a execução do alinhador, os arquivos finais foram armazenados em `/output/star_out_final/` <center> <img src="https://hackmd.io/_uploads/SJP9j0ltT.png" width=""/> </center> ## 5. Montagem do Transcriptoma A fusão de anotações de transcritos é crucial para integrar dados de montagens distintas, proporcionando uma visão abrangente da expressão gênica. Essa abordagem reduz redundâncias, aprimora a precisão da anotação e facilita análises combinadas, contribuindo para uma compreensão mais completa e precisa do transcriptoma. Para realizar a fusão das anotações dos transcritos, usou-se o `stringtie`, presente no mesmo script da etapa de alinhamento ```bash= stringtie --merge \ # Inicia o comando stringtie para a fusão de anotações de transcritos -G ${refgtf} \ # Especifica o arquivo GTF de referência que contém informações de anotação genômica para guiar a fusão -o ${outdir}/${aligner}_stringmerge/merged.gtf \ # Define o caminho de saída para o arquivo GTF resultante da fusão -m 200 \ # Define o tamanho mínimo de um gene para ser considerado durante a fusão -c 1 \ # Define o nível mínimo de cobertura para aceitar um gene na fusão -F 4 \ # Define o nível mínimo de pontuação FPKM (Fragments Per Kilobase of transcript per Million mapped reads) para aceitar um gene na fusão -T 4 \ # Define o nível mínimo de TPM (Transcripts Per Million) para aceitar um gene na fusão. -f 0.25 \ # Define o fator mínimo de fusão para fusão de transcritos sobrepostos -g 25 \ # Define o tamanho máximo do gap entre genes para permitir a fusão ${outdir}/assembly_GTF_list.txt \ # Especifica um arquivo contendo uma lista de arquivos GTF de montagens individuais a serem fundidas > ${outdir}/${aligner}_stringmerge/stringmerge.out.txt \ 2> ${outdir}/${aligner}_stringmerge/stringmerge.err.txt ``` Após a montagem das reads, os arquivos finais foram armazenados em `/output/star_stringtie/` e `/output/star_stringmerge/` <center> <img src="https://hackmd.io/_uploads/rJLpaRet6.png" width=""/> </center> > Dentro de cada subdiretório existem 10 arquivos, dentre eles arquivos de abundância e cobertura. <center> <img src="https://hackmd.io/_uploads/BJCo0AgK6.png" width=""/> </center> ## 6. Comparação com o Transcriptoma de Refência A comparação de anotações de transcritos é fundamental para validar a qualidade e precisão das montagens, identificando concordâncias e discrepâncias. Isso assegura a confiabilidade das anotações, contribuindo para uma interpretação precisa da expressão gênica e a descoberta de novos insights biológicos. Para realizar a comparação das anotações dos transcritos, usou-se o `cuffcompare`, presente no mesmo script da etapa anterior ```bash= cuffcompare -r ${refgtf} \ # Inicia o comando cuffcompare para a comparação de anotações de transcritos e especifica o arquivo GTF de referência que contém informações de anotação genômica -s ${refseq} \ # Especifica o arquivo de sequência de referência, necessário para a geração de estatísticas -o ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp \ # Define o diretório de saída e o prefixo do arquivo para os resultados da comparação ${transcriptomeref} \ # Especifica o arquivo GTF contendo as anotações de transcritos geradas pelo assembler que serão comparadas > ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.out.txt \ 2> ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.err.txt ``` Na sequência, para realizar a quantificação de expressão gênica, usou-se o `cuffquant` ```bash= cuffquant --output-dir ${outdir}/${aligner}_${assembler}_cuffquant/${name} \ # Inicia o comando cuffquant para a quantificação de expressão gênica e define o diretório de saída onde os resultados da quantificação serão armazenados --frag-bias-correct ${refseq} \ # Utiliza correção de viés de fragmento com base na sequência de referência especificada --multi-read-correct \ # Aplica correção para leituras multimapeadas, aprimorando a precisão da quantificação --num-threads ${THREADS} \ --library-type fr-firststrand \ # Especifica o tipo de biblioteca, indicando a orientação da sequência --frag-len-mean 400 \ # --frag-len-std-dev 50 \ # --max-bundle-frags 9999999 \ # Define o número máximo de fragmentos por agrupamento (bundle) --max-frag-multihits 20 \ # Define o número máximo de multimapeamentos permitidos por fragmento ${transcriptomeref} \ # Especifica o arquivo GTF contendo as anotações de transcritos para os quais a expressão será quantificada ${bamfile} \ # Especifica o arquivo BAM contendo os dados de alinhamento das leituras ao genoma de referência > ${outdir}/${aligner}_${assembler}_cuffquant/${name}/cuffquant.log.out.txt \ 2> ${outdir}/${aligner}_${assembler}_cuffquant/${name}/cuffquant.log.err.txt ``` Após a comparação e quantificação, os arquivos finais foram armazenados em `/output/star_stringtie_cuffcompare/` e `/output/star_stringtie_cuffquant/` <center> <img src="https://hackmd.io/_uploads/rJGsxkbKp.png" width=""/> </center> <center> <img src="https://hackmd.io/_uploads/S1atWJbta.png" width=""/> </center> ## 7. Análise de Expressão Diferencial A análise de diferenciação de expressão gênica é essencial para identificar genes cuja expressão varia entre diferentes condições biológicas, revelando insights sobre processos celulares e respostas a estímulos. Antes da análise de expressão diferencial, foi feita a normalização de expressão gênica com o `cuffnorm`, ainda no mesmo script ```bash= cuffnorm --output-dir ${outdir}/${aligner}_${assembler}_cuffnorm \ # Inicia o comando cuffnorm para a normalização de expressão gênica entre diferentes amostras e define o diretório de saída onde os resultados da normalização serão armazenados --labels $(IFS=, ; echo "${biogroup_label[*]}") \ # Especifica os rótulos (labels) dos grupos biológicos, obtidos a partir da variável biogroup_label. Isso é feito usando a expansão de array e a expressão de substituição de comando echo --num-threads ${THREADS} \ --library-type fr-firststrand \ # Especifica o tipo de biblioteca, indicando a orientação da sequência --library-norm-method geometric \ # Define o método de normalização usando a média geométrica --output-format simple-table \ # Especifica o formato de saída como uma tabela simples ${transcriptomeref} \ # Especifica o arquivo GTF contendo as anotações de transcritos para os quais a expressão será normalizada ${biogroup_files[*]} \ # Especifica os arquivos BAM contendo os dados de alinhamento RNA-Seq para cada grupo biológico, obtidos a partir da variável biogroup_files > ${outdir}/${aligner}_${assembler}_cuffnorm/cuffnorm.log.out.txt \ 2> ${outdir}/${aligner}_${assembler}_cuffnorm/cuffnorm.log.err.txt ``` Por último, para realizar a análise de diferenciação de expressão gênica, usou-se o `cuffdiff` ```bash= cuffdiff --output-dir ${outdir}/${aligner}_${assembler}_cuffdiff \ # Inicia o comando cuffdiff para a análise de diferenciação de expressão gênica entre diferentes condições biológicas e define o diretório de saída onde os resultados da análise serão armazenados --labels $(IFS=, ; echo "${biogroup_label[*]}") \ # Especifica os rótulos (labels) dos grupos biológicos, obtidos a partir da variável biogroup_label. Isso é feito usando a expansão de array e a expressão de substituição de comando echo --frag-bias-correct ${refseq} \ # Utiliza correção de viés de fragmento com base na sequência de referência especificada --multi-read-correct \ # Aplica correção para leituras multimapeadas, aprimorando a precisão da análise --num-threads ${THREADS} \ # --library-type fr-firststrand \ # Especifica o tipo de biblioteca, indicando a orientação da sequência --frag-len-mean 400 \ # Define o tamanho médio do fragmento de leitura --frag-len-std-dev 50 \ # Define o desvio padrão do fragmento de leitura --max-bundle-frags 9999999 \ # Define o número máximo de fragmentos por agrupamento (bundle) --max-frag-multihits 20 \ # Define o número máximo de multimapeamentos permitidos por fragmento --total-hits-norm \ # Normaliza o número total de leituras para cada biblioteca --min-reps-for-js-test 2 \ # Especifica o número mínimo de réplicas necessárias para o teste de junção suave (js-test) --library-norm-method geometric \ # Define o método de normalização usando a média geométrica --dispersion-method per-condition \ # Define o método de dispersão como "per-condition" --min-alignment-count 10 \ # Define o número mínimo de leituras alinhadas necessárias para incluir um transcrição no cálculo da expressão ${transcriptomeref} \ # Especifica o arquivo GTF contendo as anotações de transcritos para os quais a expressão será analisada ${biogroup_files[*]} \ # Especifica os arquivos BAM contendo os dados de alinhamento RNA-Seq para cada grupo biológico, obtidos a partir da variável biogroup_files > ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.out.txt \ 2> ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.err.txt ``` Após a normalização e análise de expressão diferencial, os arquivos finais foram armazenados em `/output/star_stringtie_cuffnorm/` e `/output/star_stringtie_cuffdiff/` <center> <img src="https://hackmd.io/_uploads/HJLF_1bYa.png" width=""/> </center> <center> <img src="https://hackmd.io/_uploads/BJ5Fuy-Ka.png" width=""/> </center> ## 8. Heatmap A geração de um heatmap (mapa de calor) serve para visualizar de maneira eficaz os padrões de expressão gênica diferencial entre tratamentos e controles. A representação gráfica destaca clusters de genes com perfis semelhantes, proporcionando uma compreensão visual e rápida das diferenças nas respostas genéticas às condições experimentais. Essa representação facilita a identificação de grupos de genes co-regulados, permitindo insights valiosos sobre os processos biológicos subjacentes. No presente projeto, o heatmap foi gerado no RStudio com o seguinte script ```bash= library('xlsx') library('pheatmap') # Função para leitura de dados read_data <- function(file_path) { read.delim(file_path, header = TRUE, sep = "\t", dec = ".", stringsAsFactors = FALSE) } # Função para seleção de genes diferencialmente expressos select_significant_genes <- function(genes_diff_df, threshold = 0.05) { subset(genes_diff_df, q_value <= threshold) } # Função para união de dataframes merge_dataframes <- function(genes_diff_df, c_reinhardtii_df) { merge(genes_diff_df, c_reinhardtii_df, by.x = 'test_id', by.y = 'tracking_id') } # Leitura dos dados c_reinhardtii_df <- read_data("./output/star_stringtie_cuffnorm/genes.fpkm_table") genes_diff_df <- read_data('./output/star_stringtie_cuffdiff/gene_exp.diff') # Seleção de genes diferencialmente expressos signif_genes_diff_df <- select_significant_genes(genes_diff_df) # União dos dataframes signif_c_reinhardtii_df <- merge_dataframes(signif_genes_diff_df, c_reinhardtii_df) # Escrita do arquivo Excel write.xlsx(signif_c_reinhardtii_df[, c('locus', 'sample_1', 'value_1', samples.CNTRL, 'sample_2', 'value_2', samples.TREAT, 'log2.fold_change.', 'q_value')], file = "./signif_genes.xlsx", row.names = FALSE) # Criação do dataframe para o heatmap dfh <- data.frame(group = as.character(c_reinhardtii.samps)) dfh$group <- ifelse(startsWith(as.character(c_reinhardtii.samps), "TREAT"), "Tratamento", "Controle") dfh$group <- as.factor(dfh$group) rownames(dfh) <- as.character(c_reinhardtii.samps) # Definição de cores para o heatmap ann_colors <- list(group = c('Tratamento' = 'black', 'Controle' = 'gray50')) # Configuração do heatmap e salvamento da imagem png(filename ="./c_reinhardtii_heatmap.png", res=300, width=3000, height=3000) pheatmap(as.matrix(signif_c_reinhardtii_df[, c_reinhardtii.samps]), scale = "row", annotation_col = dfh, color = colorRampPalette(c("cyan", "gold", "magenta"))(50), annotation_colors = ann_colors, annotation_names_row = TRUE, show_rownames = TRUE, cluster_rows = FALSE) dev.off() ``` O seguinte dataframe foi gerado a partir da tabela FPKM (Fragments Per Kilobase Million) armazenada no diretório `/output/star_stringtie_cuffnorm/` <center> <img src="https://hackmd.io/_uploads/Byyb0pWYp.png" width=""/> </center> Na sequência, o arquivo `/output/star_stringtie_cuffdiff/gene_exp.diff` passou pela seleção de genes diferencialmente expressos `select_significant_genes` e gerando outro datadrame <center> <img src="https://hackmd.io/_uploads/Bke2Ca-F6.png" width=""/> </center> Por fim, ambos os dataframes foram unidos. Como pode ser visto abaixo, somente o gene CHLRE_01g010400v5 (XM_001689813.2) foi significativo <center> <img src="https://hackmd.io/_uploads/H107eCZKT.png" width=""/> </center> O heatmap foi gerado em formato *.PNG <center> <img src="https://hackmd.io/_uploads/B1b-p9Wta.png" width=""/> </center> ## 9. Estrutura Final de Diretórios Ao final do projeto, 87 diretórios e 763 arquivos foram criados <center> <img src="https://hackmd.io/_uploads/B1WkHyftp.jpg" width="400"/> </center>