# **RELATÓRIO DE BIOINFORMÁTICA II - ANÁLISE DE TRANSCRITOMA**
> Universidade Estadual Paulista
Faculdade de Ciências Agrárias e Veterinárias (FCAV)
Disciplina: Bioinformártica II - Análise de transcritomas
Docente: Daniel Guariz Pinheiro
Discente: Larissa Graciano Braga
Doutoranda em Ciência Animal
Contato: larissa.graciano@unesp.br
Os scripts utilzados ao longo da disciplina estão disponíveis no [perfil do GitHub do professor](https://github.com/dgpinheiro)
Para a análise, foi escolhido o organismo *Gallus gallus* (frango):
### Criar o diretório da atividade
```
mkdir /state/partition1/lgbraga/project_chicken/data
```
Esse será o diretório principal, todos os arquivos da análise realizada estão neste diretório ou em seus subdiretórios
### Baixar o arquivo fasta e GFF do organismo escolhido
```
mkdir ref
cd ref
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/699/485/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/699/485/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.gff.gz
```
- Recuperar o cabeçalho do arquivo fasta e guardar no arquivo `genome.txt`
```
grep '^>' refs/GCF_000499845.1_PhaVulg1_0_genomic.fna | sed 's/^>\([^.]\+\)\.[0-9]\+ /\1\t/' > refs/genome.txt
```
### Retirar parte do cabeçalho que possam atrapalhar a análise (parte após o .) e colocando isso no arquivo genome.fa
```
sed 's/^>\([^.]\+\).*/>\1/' refs/GCF_000499845.1_PhaVulg1_0_genomic.fna > refs/genome.fa
```
> Antes:
NC_052532.1 Gallus gallus isolate bGalGal1 chromosome 1, bGalGal1.mat.broiler.GRCg7b, whole genome shotgun sequence
Depois:
NC_052532
### Adequar o arquivo GFF original para que os programas possam utilizá-lo
A adequação é feita por meio do script `fixNCBIgff.sh`
O arquivo atualizado passa a ser o genome.gff
```
fixNCBIgff.sh refs/GCF_000499845.1_PhaVulg1_0_genomic.gff
mv refs/GCF_000499845.1_PhaVulg1_0_genomic_fixed.gff refs/genome.gff
```
Observação: esse script é dependente de três outros scripts (`fixGFFgenestruct.pl`, `fixGFFdupID.pl`, `sortGFF.pl`)
- Verificar o número de genes anotados no NCBI para o organismo escolhido
```
grep -c -P '\tgene\t' refs/genome.gff
25403
```
- Checar tamanho máximo e mínimo de íntrons anotados no NCBI do organismo escolhido
```
checkMinMaxIntronSize.sh refs/genome.gff introninfo
1,678773
```
Min size intron: 1
Max size intron: 678773
- Verificar os percentis 0,1% e 99,9% de tamanho de íntrons
```
checkMinMaxIntronSize.sh refs/genome.gff introninfo 0.001 0.999
61,143913
```
# Simulação de dados de sequenciamento de RNA
### Escolha de de alvos (XR, XM) para geração de leituras por simulação
A identificação (ID) de RNA: prefixo XR representa os RNA não codificadores, prefixo XM representa os codificadores
Visualizar o arquivo contendo os RNA:
```
column -t -s$'\t' refs/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.gff | grep "ID=rna" | more
```
ATENÇÃO: a informação alvo (meurna e meugene) está na coluna 9, observe os itens *ID=rna-meurna;Parent=gene-meugene*.
Foram escolhidos seis RNA que estão presentes no cromossomo NC_052532.
>XR_005839958.1
XR_003076330.2
XR_003076320.2
XR_005839940.1
XR_005839933.1
XR_005839931.1
Essa informação será utilizada para produção de leituras por simulação. A partir dessa escolha, um arquivo txt (`ACCS.txt`) com o ID dos RNA escolhidos e seus respectivos genes foi criado:
Importante: O arquivo deve ser separado por tab e não deve conter linhas vazias
### Determinação de abundância da expressão dos RNA alvo previamente escolhidos
Esse procedimento é realizado para obter leituras diferencialmente expressas entre dois grupos (controle e teste). São criados dois arquivos de abundância para os RNA escolhidos (`abundance_A.txt` e `abundance_B.txt`), a diferença de expressão entre os dois grupos fica a critério do aluno
Importante: A abundância deve ser expressa em número decimal, respeitando a soma de 1. Sem linhas vazias e separação deve ser por tab
### Arquivos GFF e fasta contendo apenas os RNA escolhidos e seu cromossomo (**arquivos toy**)
```
grep -P '^(##gff|NC_052532\t)' refs/genome.gff > refs/toygenome.gff
echo -e "NC_052532" | pullseq -N -i refs/genome.fa > refs/toygenome.fa
```
### Produção de arquivos fastq simulados
Produção de três amostras por grupo com 30000 leituras cada amostra
```
cp /usr/local/bioinfo/bioaat/sim.sh ./
mkdir raw
cd raw
./sim.sh 3 30000
```
### Geração de um link simbolico para cada amostra na pasta input
Mudou-se o nome de cada a amostra para facilitar a compreensão da análise
```
mkdir input
cd input
# condicao A (contole): CTR biological replicates
ln -s ../raw/SAMPLEA1_R1.FASTQ CTR01_R1.fastq
ln -s ../raw/SAMPLEA1_R1.fastq CTR01_R2.fastq
ln -s ../raw/SAMPLEA2_R1.fastq CTR02_R1.fastq
ln -s ../raw/SAMPLEA2_R2.fastq CTR02_R2.fastq
ln -s ../raw/SAMPLEA3_R1.fastq CTR03_R1.fastq
ln -s ../raw/SAMPLEA3_R2.fastq CTR03_R2.fastq
# condicao B (teste): TST biological replicates
ln -s ../raw/SAMPLEB1_R1.fastq TST01_R1.fastq
ln -s ../raw/SAMPLEB1_R2.fastq TST01_R2.fastq
ln -s ../raw/SAMPLEB2_R1.fastq TST02_R1.fastq
ln -s ../raw/SAMPLEB2_R2.fastq TST02_R2.fastq
ln -s ../raw/SAMPLEB3_R1.fastq TST03_R1.fastq
ln -s ../raw/SAMPLEB3_R2.fastq TST03_R2.fastq
```
# Pré processamento
O pré processamento foi realizado por meio do script `preprocess5.sh`
```
mkdir -p output
./preprocess5.sh ./input ./output 6
```
Argumentos do script: $1=diretório contendo o arquivos crus, $2=diretório output, $3=número de threads a serem utilizadas
### Programas executados:
#### **1. fastqc (pre)**
Avaliação da qualidade dos arquivos crus (raw reads).
Os resultados desta etapa estão na pasta output/processed/fastqc/pre
Para avaliar os resultados, transferir os arquivos .html para o seu computador local pelo WinSCP e abrir no browser.
#### **2. atropos**
Trimagem dos adaptadores: foram utilizados dois modos do programa atropos (opções insert e adapter)
- atropos insert
```
atropos trim --aligner insert --error-rate 0.1 --times 2 --minimum-length 1 --op-order GAWCQ --match-read wildcards --overlap 3 --quality-cutoff 10 --threads ${THREADS} --correct-mismatches conservative --pair-filter any -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq -p ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq -pe1 ${r1} -pe2 ${r2} --untrimmed-output ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq --untrimmed-paired-output ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq > ${outdir}/processed/atropos/${name}.atropos.log.out.txt
```
- atropos adapter
```
atropos trim --aligner adapter -e 0.1 -n 2 -m 1 --match-read-wildcards -O 3 -q 10 --pair-filter any -pe1 ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq -pe2 ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -g AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -G CAAGCAGAAGACGGCATACGAGATNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT -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
```
> --error-rate/-e (padrão 0.1)
> --times/-n 2 :applies to all the adapters given (padrão 1)
> --minimum-length/-m 1: Throw away processed reads shorter than N bases.
> --op-order GAWCQ: adapter trimming, C = cutting (unconditional), G = NextSeq trimming, Q = quality trimming, W = overwrite poor quality reads (padrão CGQAW)
Os resultados desta etapa estão na pasta `output/processed/atropos`
#### **3. prinseq**
Trimagem de regiões de baixa qualidade, remoção de cauda polyA. Os resultados desta etapa estão na pasta `output/processed/prinseq`
```
prinseq-lite.pl -fastq ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq -fastq2 ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq -out_format 3 -trim_qual_window 3 -trim_qual_step 1 -trim_qual_type mean -trim_qual_rule lt -trim_qual_right 30 -out_good ${outdir}/processed/prinseq/${name}.atropos_final.prinseq -out_bad null -lc_method dust -lc_threshold 50 -min_len 20 -trim_tail_right 5 -trim_tail_left 5 -ns_max_p 80 > ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.out.log
```
As leituras que não possuem pares estão no arquivo com sufixo `_singletons.fastq`
#### **4. fastqc (pos)**
Avaliação da qualidade dos arquivos pré-processados.
Os resultados desta etapa estão na pasta `output/processed/fastqc/pos`
Para avaliar os resultados, refazer procedimento anterior
### Resultado final do pré processamento
Links simbólicos das amostras após a etapa do prinseq foram criados na pasta `output/processed/cleaned`
Esses arquivos serão utilizados na próxima etapa
- Verificar quantas reads foram excluídas no pré-processamento
Exemplo para uma amostra:
```
echo "scale=2; $(cat input/CTR01_R1.fastq | wc -l)/4" | bc
30001.00 # número de leitura antes do pré-processamento
echo "scale=2; $(cat output/processed/cleaned/CTR01.atropos_final.prinseq.cleaned_1.fastq | wc -l)/4" | bc
29366.00 # número de leitura após pré-processamento
```
Observação: em dados reais o R2 possui pior qualidade.
- Verificando qualidade
```
cp /data/tmp/bioaat/stats_preprocess.sh ./
./stats_preprocess.sh ./input ./output > STATS.txt
```
# Alinhamento ao genoma de referência
### Preparação
- Taxonomy ID da organismo escolhido
Verificar qual o número de Taxonomy NCBI do organismo escolhido (https://www.ncbi.nlm.nih.gov/Taxonomy)
```
mkdir home/lgbraga/ex-lista1
cd ex-lista1
wget ftp://ftp.ncbi.nih.gov/gene/DATA/gene_info.gz
gunzip gene_info.gz
```
O arquivo `gene_info` contém informação sobre o organismo escolhido
- Converter GFF3 para GTF
Utilizo o arquivo GTF como guia para o alinhamento
```
gffread toygenome.gff -T -o toygenome.gtf
```
### 1. Alinhamento
Linha de comando utilizada (atualizada para a última versão do `rnaseq-ref.sh` utilizada na disciplina)
```
cp /data/tmp/bioaat/rnaseq-ref.sh ./
./rnaseq-ref.sh star ./input ./output 6 ./refs/toygenome.gff ./refs/toygenome.fa stringtie NA /home/lgbraga/ex-lista1/gene_info 9031
```
#### A. Tophat2
- Indexar genoma (com Bowtie)
```
bowtie2-build --threads ${THREADS} genome.fa genome > bowtie2.out.txt
```
- Alinhamento paired-end:
```
tophat2 --min-anchor 8 --min-intron-length 50 --max-intron-length 500000 --max-multihits 20 --transcriptome-max-hits 10 --prefilter-multihits --num-threads ${THREADS} --GTF ${refgtf} --transcriptome-index ${outdir}/tophat_index/transcriptome --mate-inner-dist 50 --mate-std-dev 50 --coverage-search --microexon-search --b2-very-sensitive --library-type fr-unstranded --output-dir ${outdir}/tophat_out_pe/${name} --no-sort-bam =${outdir}/tophat_index/genome ${r1} ${r2} > ${outdir}/tophat_out_pe/${name}.log.out.txt
```
- Alinhamento single-end:
```
tophat2 --min-anchor 8 --min-intron-length 50 --max-intron-length 500000 --max-multihits 20 --transcriptome-max-hits 10 --prefilter-multihits --num-threads ${THREADS} --GTF ${refgtf} --transcriptome-index ${outdir}/tophat_index/transcriptome --mate-inner-dist 50 --mate-std-dev 50 --coverage-search --microexon-search --b2-very-sensitive --library-type fr-unstranded --output-dir ${outdir}/tophat_out_se/${name} --no-sort-bam ${outdir}/tophat_index/genome ${outdir}/tophat_out_se/${name}/singletons.fastq > ${outdir}/tophat_out_se/${name}.log.out.txt
```
#### **B. STAR**
O professor utiliza as opções com parâmetros usado pelo projeto Encode, entretanto, algumas dessas opções foram alteradas com bases nos critérios da aluna.
- Indexar genoma com STAR
```
STAR --runThreadN ${THREADS} --runMode genomeGenerate --genomeFastaFiles genome.fa --genomeDir ./ --sjdbGTFfile ${absrefgtf} --genomeSAindexNbases 12 --sjdbOverhang 149 > STAR.genomeGenerate.log.out.txt
```
- Alinhamento paired-end
```
STAR --runThreadN ${THREADS} --genomeDir ${outdir}/star_index/ --readFilesIn ${r1} ${r2} --outSAMstrandField intronMotif --outFilterIntronMotifs RemoveNoncanonical --sjdbGTFfile ${refgtf} --outFilterMultimapNmax 20 --outFileNamePrefix ${outdir}/star_out_pe/${name}/ --outSAMtype BAM Unsorted --outFilterType BySJout --outSJfilterReads Unique --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.03 --alignIntronMin 50 --alignIntronMax 500000 --alignMatesGapMax 5000 > ${outdir}/star_out_pe/${name}.log.out.txt
```
- Alinhamento single-end
```
STAR --runThreadN ${THREADS} --genomeDir ${outdir}/star_index/ --readFilesIn ${r1} --outSAMstrandField intronMotif --outFilterIntronMotifs RemoveNoncanonical --sjdbGTFfile ${refgtf} --outFilterMultimapNmax 20 --outFileNamePrefix ${outdir}/star_out_se/${name}/ --outSAMtype BAM Unsorted --outFilterType BySJout --outSJfilterReads Unique --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.03 --alignIntronMin 50 --alignIntronMax 500000 --alignMatesGapMax 5000 > ${outdir}/star_out_se/${name}.log.out.txt
```
Observação: para a execução do cufflinks é necessário: --outSAMstrandField intronMotif e --outFilterIntronMotifs RemoveNoncanonical
#### C. Merge
Caso haja resultado do alinhamento single-end para qualquer um dos alinhadores, os arquivos contendo alinhamento paired-end e single-end são unidos. Primeiro o cabeçalho é recuperado pela opção view e após acontece a união pela opção merge, ambas opções do programa samtools
```
samtools view -H ${outdir}/star_out_pe/${name}/Aligned.out.bam > ${outdir}/star_out_final/${name}/Header.txt
samtools merge -n --threads ${THREADS} -h ${outdir}/star_out_final/${name}/Header.txt ${outdir}/star_out_final/${name}/Aligned.out.bam ${outdir}/star_out_pe/${name}/Aligned.out.bam ${outdir}/star_out_se/${name}/Aligned.out.bam > ${outdir}/star_out_final/${name}.log.out.txt 2> ${outdir}/star_out_final/${name}.log.err.txt
```
#### D. Pós alinhamento
#### Avaliando a qualidade das amostras alinhadas ao genoma de referência
```
./getAlnTab.pl --indir output/stats_tophat/ --aligner tophat --s ./STATS.txt > output/tophat_AlnTab.txt
```
```
./getAlnTab.pl --indir output/stats_star/ --aligner star --s ./STATS.txt > output/star_AlnTab.txt
```
> No comparativo dos dois alinhadores, o STAR obteve:
Maior número de leituras paired-end;
Menor número de alinhamento não concordantes;
Maior número de alinhamentos únicos concordantes;
Maior número de alinhamentos concordantes múltiplos Menor número de leituras singletons;
A partir desse resultado, o alinhamento do **STAR** foi considerado como de melhor resultado e foi escolhido para as próximas etapas
#### Visualizar os alinhamentos
Utilizar bam ordenada (`sorted.bam`) e indexada (`sorted.bam.bai`)
```
samtools sort my.bam > my.sorted.bam
samtools index my.sorted.bam
```
Arquivo .fa e .gtf devem estar indexados
Descubrir a região alvo para visualizar o alinhamento (`grep "my_transcrpts_id" refs/toygenome.gtf`)
Anotar as posições inicial e final dos genes alvo
Passar arquivos .gtf e .fa para a sua máquina local pelo WinSCP
- IGV
Seguir os passos abaixo para visualizar o alinhamento:
*Genomes* (menu superior) - *load genome from file* e selecionar arquivo `toygenoma.fa`;
*File* (menu superior)- *load from file* para selecionar arquivo `toygenome.gtf`;
(menu lateral) - clicar no `toygenome.gtf` com botão direito e selecionar opção *Expanded*;
*File* (menu superior) - *load from file* para slecionar `my.sorted.bam` para visualizar no IGV;
Selecionar as regiões alvo que deseja visualizar abaixo do menu superior e clicar em *Go*;
Ajustar o zoom por meio de uma barra móvel no canto superior direito da tela.
Exemplo:
Opções para gerar imagem: Sashimi plot, salvar como arquivo SVG (imagem vetorial), editar no programa inkscape.
### 2. Sort
Gerar arquivos "sort-coordinated"
```
samtools sort --threads ${THREADS} -o ${outdir}/star_out_final/${name}/Aligned.sorted.bam ${outdir}/star_out_final/${name}/Aligned.out.bam > ${outdir}/star_out_final/${name}/Aligned.sorted.out.txt 2> ${outdir}/star_out_final/${name}/Aligned.sorted.err.txt
```
O resultado final é encontrado na pasta `${outdir}/star_out_final/${name}/Aligned.out.bam`
### 3. Montagem do transcriptoma
#### A. Cufflinks
- Montagem do transcritoma
```
cufflinks --output-dir ${outdir}/${aligner}_cufflinks/${name} --num-threads ${THREADS} --GTF-guide ${refgtf} --frag-bias-correct ${refseq} --multi-read-correct --library-type fr-unstranded --frag-len-mean 300 --frag-len-std-dev 30 --total-hits-norm --min-isoform-fraction 0.25 --pre-mrna-fraction 0.15 --min-frags-per-transfrag 10 --junc-alpha 0.001 --small-anchor-fraction 0.08 --overhang-tolerance 100 --min-intron-length 10 --max-intron-length 1050000 --trim-3-avgcov-thresh 0.05 --trim-3-dropoff-frac 0.01 --max-multiread-fraction 0.75 --overlap-radius 25 --3-overhang-tolerance 600 --intron-overhang-tolerance 50 ${outdir}/align_out_final/${name}/Aligned.sorted.bam > ${outdir}/${aligner}_cufflinks/${name}/cufflinks.out.txt 2> ${outdir}/${aligner}_cufflinks/${name}/cufflinks.err.txt
```
- União com o transcritoma referência pelo cuffmerge
```
cuffmerge -o ${outdir}/${aligner}_cuffmerge --ref-gtf ${refgtf} --ref-sequence ${refseq} --min-isoform-fraction 0.25 --num-threads ${THREADS} ${outdir}/assembly_GTF_list.txt > ${outdir}/${aligner}_cuffmerge/cuffmerge.out.txt 2> ${outdir}/${aligner}_cuffmerge/cuffmerge.err.txt
```
#### **B. StringTie**
- Montagem do transcritoma
```
stringtie ${outdir}/align_out_final/${name}/Aligned.sorted.bam -G ${refgtf} -f 0.15 -m 200 -o ${outdir}/${aligner}_stringtie/${name}/transcripts.gtf -a 8 -j 4 -c 4 -v -g 20 -C ${outdir}/${aligner}_stringtie/${name}/coverages.txt -M 1 -p ${THREADS} -A ${outdir}/${aligner}_stringtie/${name}/abundances.txt -B > ${outdir}/${aligner}_stringtie/${name}/stringtie.out.txt 2> ${outdir}/${aligner}_stringtie/${name}/stringtie.err.txt
```
- União com o transcritoma referência: opção merge do programa Stringtie
```
stringtie --merge -G ${refgtf} -o ${outdir}/${aligner}_stringmerge/merged.gtf -m 200 -c 4 -F 4 -T 4 -f 0.15 -g 20 ${outdir}/assembly_GTF_list.txt > ${outdir}/${aligner}_stringmerge/stringmerge.out.txt 2> ${outdir}/${aligner}_stringmerge/stringmerge.err.txt
```
A montagem realizada pelo **Stringtie** foi escolhida para ser utilizada nas etapas posteriores da análise
### 4. Cuffcompare
```
cuffcompare -r ${refgtf} -s ${refseq} -o ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp ${transcriptomeref} > ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.err.txt
```
Resultado final principal em `${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.combined.gtf`
### 5. Quantificação da expressão
- Cuffquant
```
cuffquant --output-dir ${outdir}/${aligner}_${assembler}_cuffquant/${name} --frag-bias-correct ${refseq} --multi-read-correct --num-threads ${THREADS} --library-type fr-unstranded --frag-len-mean 300 --frag-len-std-dev 30 --max-bundle-frags 9999999 --max-frag-multihits 20 ${transcriptomeref} ${bamfile} > ${outdir}/${aligner}_${assembler}_cuffquant/${name}/cuffquant.log.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffquant/${name}/cuffquant.log.err.txt
```
Resultado final principal em `${outdir}/${aligner}_${assembler}_cuffquant/${name}/abundances.cxb`
### 6. Análise de expressão gênica
#### A. Cuffnorm
Gerar matriz de abundância normalizada
```
cuffnorm --output-dir ${outdir}/${aligner}_${assembler}_cuffnorm --labels $(IFS=, ; echo "${biogroup_label[*]}") --num-threads ${THREADS} --library-type fr-unstranded --library-norm-method geometric --output-format simple-table ${transcriptomeref} ${biogroup_files[*]} > ${outdir}/${aligner}_${assembler}_cuffnorm/cuffnorm.log.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffnorm/cuffnorm.log.err.txt
```
Resultado final principal em `${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.count_table`
#### B. De-normalização
Para ser utilizada em outros programas de análise de expressão gênica (ex: DeSeq2, EdgeR)
B1. De-normalização de transcritos
```
de-normalize-cuffnorm.R --in=${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.count_table --st=${outdir}/${aligner}_${assembler}_cuffnorm/samples.table --out=${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.raw_count_table.txt > ${outdir}/${aligner}_${assembler}_cuffnorm/de-normalize-cuffnorm.isoforms.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffnorm/de-normalize-cuffnorm.isoforms.err.txt
```
Resultado final principal em
`${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.raw_count_table.txt`
- Contagem de transcritos sem normalização
Arquivo `/output/star_stringtie_cuffnorm/isoforms.raw_count_table.txt`)
> Também foram encontradas contagens para transcritos que não foram escolhidos
B2. De-normalização de genes
```
de-normalize-cuffnorm.R --in=${outdir}/${aligner}_${assembler}_cuffnorm/genes.count_table --st=${outdir}/${aligner}_${assembler}_cuffnorm/samples.table --out=${outdir}/${aligner}_${assembler}_cuffnorm/genes.raw_count_table.txt > ${outdir}/${aligner}_${assembler}_cuffnorm/de-normalize-cuffnorm.genes.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffnorm/de-normalize-cuffnorm.genes.err.txt
```
Resultado final principal em `${outdir}/${aligner}_${assembler}_cuffnorm/genes.raw_count_table.txt`
### 7. Análise de expressão diferencial
#### A. Cuffdiff
```
cuffdiff --output-dir ${outdir}/${aligner}_${assembler}_cuffdiff --labels $(IFS=, ; echo "${biogroup_label[*]}") \--frag-bias-correct ${refseq} --multi-read-correct --num-threads ${THREADS} --library-type fr-unstranded --frag-len-mean 300 --frag-len-std-dev 30 --max-bundle-frags 9999999 --max-frag-multihits 20 --total-hits-norm --min-reps-for-js-test 2 --library-norm-method geometric --dispersion-method per-condition --min-alignment-count 10 ${transcriptomeref} ${biogroup_files[*]} > ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.out.txt 2> ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.err.txt
```
Resultado final principal em `${outdir}/${aligner}_${assembler}_cuffdiff/isoform_exp.diff`
- Criando arquivo com valores de expressão diferencial
Nos outputs deste programa, tenho a coluna *teste_id* para os RNA, mas com a Id própria do programa. Isso dificulta a compreensão e visualização do resultados dos RNA transcritos
Então, é necessário realizar uma junção (merge/join) entre o arquivo de maior interesse (`isoform_exp.diff`) e outro arquivo que tenha a correspondência do Id do programa e do Id original do RNA. Isso foi feito pelo script mergeR.R
```
mergeR.R --x=output/star_stringtie_cuffdiff/isoform_exp.diff --y=output/star_stringtie_cuffcompare/TCONS.nearest_ref.txt --by.x="" --by.y="" --all.x --print.out.label --out=output/star_stringtie_cuffdiff/isoform_exp.diff.annot.txt
splitteR.R --x="./gene_exp.diff" --col.x="gene"
```
O arquivo `isoform_exp.diff.annot.txt` contém informação do valor de expressão diferencial e id do RNA transcrito
- Converter arquivo gerado na etapa anterior para formato xls (legível no excel)
```
tsv2xls.pl -i isoform_exp.diff.annot.txt -o isoform_exp.diff.annot -h
```
Transferir arquivo `output/star_stringtie_cuffdiff/isoform_exp.diff.annot.xls` para máquina local com WinSCP para observação
#### B. EdgeR
A pipeline realizada pelo professor informa os valores de expressão normalizados, assim, será preciso desnormalizar os valores e realizar a análise pelo edgeR
- Montar arquivo groups.txt
Arquivo `groups.txt` com informação do id do RNA, nome da amostra e grupo a qual a amostra pertence (controle ou teste). As colunas devem ser nomeadas como *id*, *name* e *group*, respectivamente
- Normalização pelo edgeR pelo script do professor
```
mkdir output/edger
run-edgeR.R --in=output/star_stringtie_cuffnorm/isoforms.raw_count_table.txt --groups=./groups.txt --out=./output/edger
```
- Merge/join de arquivos
Devido ao mesmo motivo da etapa com o Cuffdiff, para facilitar a visualização
```
mergeR.R --x=./output/edger/CTR-TST.txt --by.x="tracking_id" --y="./output/star_stringtie_cuffcompare/TCONS.nearest_ref.txt" --by.y="transcript_id" --all.x --print.out.label --out=./output/edger/CTR-TST.annot.txt
```
- Converter para formato xls para facilitat visualização
```
tsv2xls.pl -i ./output/edger/CTR-TST.annot.txt -o ./output/edger/CTR-TST.annot --header
```
Transferir arquivo `output/edger/CTR-TST.annot.xls` para máquina local
#### C. Deseq2
Seguir os mesmos passos da etapa anterior
Comandos utilizados:
```
mkdir output/deseq2
run-DESeq2.R --in=output/star_stringtie_cuffnorm/isoforms.raw_count_table.txt --groups=./groups.txt --out=./output/deseq2/
```
```
mergeR.R --x=./output/deseq2/CTR-TST.txt --by.x="tracking_id" --y="./output/star_stringtie_cuffcompare/TCONS.nearest_ref.txt" --by.y="transcript_id" --all.x --print.out.label --out=./output/deseq2/CTR-TST.annot.txt
```
```
tsv2xls.pl -i ./output/deseq2/CTR-TST.annot.txt -o ./output/deseq2/CTR-TST.annot --header
```
## Verificar diferença entre resultados da expressão diferencial
- Anotar as id dos RNA escolhidos anteriormente
>XR_005839958.1
XR_003076330.2
XR_003076320.2
XR_005839940.1
XR_005839933.1
XR_005839931.1
Há RNA que não foram escolhidos, mas foram detectados na análise de expressão diferencial
- Montar planilha contendo valores esperados e observados para comparação para os programas utilizados (Cuffdiff, EdgeR e Deseq2)
Para esses dados e com esse processamento, a análise de expressão diferencial pelo **Deseq2** e **EdgeR** foram mais condizentes com o esperado
# Análise com montagem de novo
**Sequenciamento de novo**
Montagem de fragmentos de sequências dos transcritos originais por meio de um consenso
A montagem é feira pelo programa trinity por meio de três etapas:
1. Inchworm (lagarta)
2. Chrysalis (crisálida)
3. Butterfly (borboleta)
### Montagem de novo dos arquivos previamente pré-processados
```
cp /data/tmp/bioaat/rnaseq-denovo.sh ./
mkdir output/trinity
```
`./rnaseq-denovo.sh output/processed/cleaned output/trinity/ 6 12G`
#### Comparando quais transcritos foram simulados e quais foram observados
```
grep '^>' refs/transcriptoma.fa
grep '>' output/trinity/trinity_assembled/Trinity.fasta | wc -l
```
Nota: O arquivo `output/trinity/trinity_assembled/Trinity.fasta` contém as montagens de novo realizadas, e foram realizadas montagens em 20 transcritos (o esperado era apenas seis transcritos) - possível causa: similaridade entre as sequência de nucletídeos dos transcritos
- BLAST
Na pipeline do alinhamento de referência, os transcritos foram identificados pela sequência, arquivo de montagem (.gtf) e pela posição. Na montagem de novo eu não tenho com o que comparar, então para saber o que eu consegui eu faço blast das sequências obtidas com as sequências do RNA escolhidos para simulação
A sequência dos RNA escolhidos para a simulação está no arquivo `refs/transcriptoma.fa`
Esse arquivo será o banco de dados para o blast, e deve estar indexado, para indexar:
```
makeblastdb -title "ReferenceTranscriptome" -dbtype nucl -out RefTrans -in trascriptoma.fa
```
Esse comando para indexação de banco de dados de nucleotídeo gerou os arquivos `RefTrans.nhr`, `RefTrans.nin`, `RefTrans.nsq`, que não são legíveis. Esses arquivos contêm a sequência dos transcritos escolhidos para a simulação
Comparação entre os transcritos que foram montados e os escolhidos para simulação (esperado):
```
blastn -max_hsps 1 -max_target_seqs 1 -query ./output/trinity/trinity_assembled/Trinity.fasta -db ./RefTrans -out ./Trinity_x_RefTrans.txt -evalue 1e-5 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore
```
> Opções: 6 = formato tabular
> Colunas: identificação da sequência query, sequencia referência, porcentagem de identidade, tamanho, quantidade de mismatches, abertura de gaps, começo do alinhamento na query, final do alinhamento na query, início do alinhamento na sequência do banco de dados, final do alinhamento na sequência do banco de dados, e-value, bit score
O arquivo de saída é o `Trinity_x_RefTrans.txt`, ele contem informações sobre 20 transcritos encontrados na minha montagem
### Anotação funcional da montagem - regiões codificadoras
- Verificar quais isoformas são codificadoras de proteínas
A opção LongOrfs do programa TransDecoder seleciona apenas as ORF mais longas (que também são as mais prováveis)
```
TransDecoder.LongOrfs -t ./output/trinity/trinity_assembled/Trinity.fasta -m 90 --gene_trans_map output/trinity/trinity_assembled/Trinity.fasta.gene_trans_map
grep -c '^>' Trinity.fasta.transdecoder_dir/longest_orfs.pep
```
O arquivo de saída `Trinity.fasta.transdecoder_dir/longest_orfs.pep` possui as sequências que são codificadoras de proteínas (ORF). Foram encontradas 150 ORF longas
- Seleção da ORF mais provável por transcrito
Essa predição é baseada nas características esperadas de uma ORF. Foram geradas muitas ORF (N=150) por transcritos (N=20) e isso pode indicar alto número de falsos positivos. Por isso e para facilitar análise, apenas as ORF mais bem avaliadas serão utilizadas para as próximas etapas:
```
TransDecoder.Predict -t ./output/trinity/trinity_assembled/Trinity.fasta --single_best_only
grep -c '^>' Trinity.fasta.transdecoder.pep
```
O arquivo de saída contem uma ORF por transcrito identificado na montagem de novo, foram identifcadas 17 ORF
#### Análise de abundância
```
cp /data/tmp/bioaat/rnaseq-trinity-abundance.sh ./ ./rnaseq-trinity-abundance.sh output/trinity/renamed output/trinity/trinity_assembled output/trinity
```
Umas das etapas do script é rodar o programa Salmon (outra opção é o Calisto). Salmon realiza pseudo-alinhamento para estimar abundância a partir de arquivo `fastq`
É diferente de um alinhador: o alinhamento é feito base-a-base, o pseudo-alinhamento realiza o mapeamento de kmers no arquivo referência. A estimativa de abundância é baseada na cobertura de kmers que ele encontra no alvo
```
less -S output/trinity/abundance/DEG/CTR-TST.txt | more
```
### Anotação funcional dos transcritos montados - proteína
- eggNOGmapper
Anotação e identificação de grupos ortólogos
```
pyenv activate eggnog-mapper
emapper.py -m diamond --no_annot --no_file_comments --cpu 6 -i Trinity.fasta.transdecoder.pep -o Trinity.fasta.transdecoder.pep --dmnd_iterate no
more Trinity.fasta.transdecodeer.pep.emapper.seed_orthologs
```
### Anotação funcional dos transcritos montados - termos Gene Ontology
```
pyenv activate eggnog-mapper
emapper.py --annotate_hits_table Trinity.fasta.transdecoder.pep.emapper.seed_orthologs --no_file_comments -o emapper.out --cpu 6 --override
column -t -s$'\t' emapper.out.emapper.annotations | less -S
```
- Retirando _i1.p1 da identificação (coluna 1):
```
more gene_annotation.txt
cp /data/tmp/bioaat/twocoldescmerger.pl ./
cut -f 1,8 emapper.out.emapper.annotations | sed 's/_i[0-9]\+\.p1//' | ./twocoldescmerger.pl > gene_annotation.txt
```
> Antes TRINITY_DN0_c0_g1_i1.p1
Depois TRINITY_DN0_c0_g1
O arquivo contem `gene_annotation.txt` o gene no qual foi encontrado sítio funcional. apesar de terem sido escolhidos transcritos pertencentes a dois genes, apenas um estão anotado nessa etapa
- Merge para facilitar identificação
```
mergeR.R --x="./output/trinity/abundance/DEG/CTR-TST.txt" --by.x="X" --y="./gene_annotation.txt" --by.y="X.query" --out="./output/trinity/abundance/DEG/CTR-TST.annot.txt" --print.out.label
```
O arquivo `output/trinity/abundance/DEG/CTR-TST.annot.txt` contem informações sobre medida normalizada de cada amostra, Log2FoldChange, pvalue, pvalue ajustado e descrição do gene no qual foi encontrado sítio funcional
# Extra
É possível gerar heatmap entre os grupos avaliados na análise de expressão diferencial. Para isso utilizar o R e comando abaixo:
`heatmap.2( as.matrix( golub.norm.sig[, golub.samps]), na.rm=TRUE, scale="row",density.info="none",trace="none",main="Heatmap",symbreaks=TRUE,col=greenred, distfun = function(x) { dist(x, method="canberra") },hclustfun=function(x) hclust(x, method="average") )`
Sign in with Wallet
Connect another wallet