---
# System prepended metadata

title: 'Relatório Final - Isabella '

---

# Relatório Final

**Disciplina:** Bioinformática Aplicada II: Análise de Transcritomas
**Discente:** Isabella Cardeal Campos
**Docente:** Daniel Guariz Pinheiro

## Introdução
A transcritoma refere-se ao conjunto de transcritos de um organismo, células, tecidos ou órgãos sob determinadas circunstâncias, ou seja, é o reflexo direto da expressão dos genes em resposta a diferentes condições. 

## Requisitos para a simulação
Seguir as instruções para instalação do e-direct.

Caso precise atualizar o conjunto de programas do e-direct, remova primeiro a pasta correspondente
```python=
cd ~
rm -fr ./edirect
```
Copiar e colar o código a seguir
```python=
cd ~
sh -c "$(curl -fsSL ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/install-edirect.sh)"
```
Adicionar o caminho no $PATH caso não tenha solicitado que seja feito automaticamente no final do passo anterior

```python=
echo "export PATH=\$PATH:\$HOME/edirect" >> $HOME/.bash_profile
```
### Estrutura de diretórios utilizada
![Untitled diagram-2023-12-21-164916](https://hackmd.io/_uploads/ry8MTJzvT.svg)


Link para visualização da imagem acima: 
https://www.mermaidchart.com/raw/13bbe15f-0eb2-45c9-ad34-523dc4e8b9e5?theme=light&version=v0.1&format=svg

## Construção do transcritoma para simulação
### Download dos arquivos FASTA e GFF

Nessa simulação trabalhamos com a bactéria ***Salmonella* *enterica* subsp. *enterica* serovar Gallinarum str. 287/91** (NC_011274.1). 

Busca de informações sobre a sequência NC_011274.1 no banco de dados de nucleotídeos do NCBI. Em seguida, o comando ```efetch``` realiza o download da sequência no formato FASTA:

```python=
esearch -db nucleotide -query NC_011274.1 | efetch -format fasta
```

Também foi acessado o banco de dados disponível no NCBI FTP site:
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/525/GCF_000009525.1_ASM952v1/

Foi realizado o download dos arquivos FASTA e GFF compactados, seguidos de descompactação
```python=
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/525/GCF_000009525.1_ASM952v1/GCF_000009525.1_ASM952v1_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/525/GCF_000009525.1_ASM952v1/GCF_000009525.1_ASM952v1_genomic.gff.gz
gunzip *.gz
```

### Gerando os arquivos .fastq
```python=
simLib.pl -a ../ref/abundance_A.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA1
simLib.pl -a ../ref/abundance_A.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA2
simLib.pl -a ../ref/abundance_A.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA3
simLib.pl -a ../ref/abundance_B.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB1
simLib.pl -a ../ref/abundance_B.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB2
simLib.pl -a ../ref/abundance_B.txt -i ../ref/transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB1
```

Para escolher os genes que serão utilizados para a análise clicar na página do genoma no NCBI em grafics>download>CSV e prosseguir para a criação do arquivo **ACCS.txt** que será utilizado para as análises

**Obs.:** O número de início (start) e fim (stop) têm que ser ajustados. Devemos considerar 1 base a menos no início e uma base a menos no fim.

A seguir o arquivo construído com os IDs dos transcritos e o nome dos genes na segunda coluna (as colunas devem ser separadas por <TAB>. Isso auxiliará no futuro a identificar quais transcritos relacionam-se com quais genes.

```
tal     NC_011274.1     7554    8507    plus
satP    NC_011274.1     9265    9831    minus
bcfB    NC_011274.1     25000   25686   plus
nhaA    NC_011274.1     46074   47240   plus
lspA    NC_011274.1     56572   57072   plus
citS    NC_011274.1     66929   68267   minus
```
Em seguida foi utilizado o seguinte script getTranscriptome.sh

```python=
#!/bin/bash
rm -f transcriptoma.fa

IFS=$'\n'
for accline in $(cat ./ACCS.txt); do
        acc=`echo ${accline} | cut -f 1`
        seqref=`echo ${accline} | cut -f 2`
        chr_start=`echo ${accline} | cut -f 3`
        chr_stop=`echo ${accline} | cut -f 4`
        strand=`echo ${accline} | cut -f 5`

        echo "Pegando FASTA para ${acc}  [${seqref}:${chr_start}-${chr_stop}(${strand})] ..."

        efetch -db nucleotide -id ${seqref}  -format fasta \
        -chr_start ${chr_start} \
        -chr_stop ${chr_stop} \
        -strand ${strand} | \
        sed "s/^>.*/>${acc}/" \
        >> transcriptoma.fa
done
```
Executar o script e visualizar o arquivo **transcriptoma.fa** gerado
```python=
chmod a+x analise.sh
./analise.sh
less transcriptoma.fa
grep '^>' transcriptoma.fa
```


### Criar arquivos de abundância

abundance_A.txt
```
tal     0.1
satP    0.2
bcfB    0.35
nhaA    0.05
lspA    0.15
citS    0.15
```

abundance_B.txt
```
tal     0.2
satP    0.1
bcfB    0.15
nhaA    0.15
lspA    0.35
citS    0.05
```

Obs.: sempre separar colunas por TAB.

Execute os comandos abaixo, o resultado da soma dos valores da coluna 2 tem que ser 1 para ambos os arquivos:
```python=
perl -F"\t" -lane 'INIT{$sum=0;} $sum+=$F[1]; END{print "SOMA: $sum"; } ' 
 abundance_A.txt
perl -F"\t" -lane 'INIT{$sum=0;} $sum+=$F[1]; END{print "SOMA: $sum"; } ' 
 abundance_B.txt
```
## Limpeza dos arquivos referência
### Limpeza do cabeçalho FASTA
Permite a padronização dos cabeçalhos, visto que estes podem conter informações como identificadores de sequências, informações sobre a fonte da sequência, descrições, caracteres especiais, espaços em branco desnecessários, entre outros. Além disso, algumas ferramentas de bioinformática e pipelines podem ser sensíveis ao formato dos cabeçalhos. 

Criar arquivo com o script (cleanfasta.sh) abaixo:
```python=
#!/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}
```

Permissão seguido da execução
```python=
chmod a+x cleanfasta.sh
./cleanfasta.sh GCF_000009525.1_ASM952v1_genomic.fna > genome.fa
```

### Ajustes do arquivo GFF
Executar o script fixNCBIgff.sh do repositório Github:
https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixNCBIgff.sh

```python=
fixNCBIgff.sh GCF_000009525.1_ASM952v1_genomic.gff genome.gff
```
Obs.: deixa o arquivo GFF adequado para rodar, só não resolve situações de *trans-splicing*.

* Conferir os cabeçalhos dos arquivos FASTA E GFF
```python=
grep ‘^>’ genome.fa
cut -f 1 genome.gff | sort -u
```
* Conversão de um arquivo GFF em GTF (algumas análises precisam dessa conversão)
```python=
gffread -g genome.fa --gtf -o genome.gtf genome.gff
```
* Para visualizar em colunas
```python=
column -t -s$'\t' genome.gtf | less -S
```
### Testes com GFF (gerar o arquivo do transcritoma e proteoma)
```python=
gffread -g genome.fa -w fulltranscriptome.fa genome.gff
gffread -g genome.fa -y protein.fa genome.gff
```

## Verificar a presença de íntrons 
No caso de procariotos não há a presença de íntrons, portanto, quando rodarmos essa análise o valor obtido deverá ser zero. 

A verificação da presença de íntrons é uma etapa importante durante a análise, pois fornece informações sobre a estrutura gênica e o processamento do RNA, como a identificação de *splicing* alternativo, variação na estrutura dos genes, detecção de fusões de genes e rearranjos cromossômicos, entre outros. 

Entrar no ambiente do python2
```python=
source /usr/local/bioinfo/anaconda3/etc/profile.d/conda.sh
conda activate python2
```
Executar dentro do diretório ```ref/``` 
```python=
checkMinMaxIntronSize.sh genome.gff . 0.25 0.7 5
```
Visualizar dados do arquivo bam
```python=
samtools view tophat_out/accepted_hits.bam | cut -f 1 | sort | uniq -c | grep -P '^ *2 '
```

Visualizar sequência FASTA
```python=
samtools faidx genome.fa
samtools faidx genome.fa NC_011274.1:66929-68267>NC_011274.1:66929-68267
```

## Atropos 
Utilizado para clivagem de adaptadores 

Criação dos diretórios
```python=
mkdir -p  output/processed/fastqc/pre
tree output/
```
Analisar a qualidade das leituras (*reads*) obtidas no sequenciamento
```python=
fastqc --threads 5 -o output/processed/fastqc/pre/ raw/SAMPLEA1_R1.fastq
fastqc --threads 5 -o output/processed/fastqc/pre/ raw/SAMPLEA1_R2.fastq
```

Execução do atropos
```python=
atropos trim --aligner insert \
             -e 0.1 \
             -n 2 \
             -m 1 \
             --op-order GAWCQ \
             --match-read-wildcards \
             -O 20 \
             -q 25 \
             -T 2 \
             --correct-mismatches conservative \
             --pair-filter any \
             -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
             -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT  \
             -o ./output/processed/atropos/SAMPLEA1_R1.atropos_insert.fastq \
             -p ./output/processed/atropos/SAMPLEA1_R2.atropos_insert.fastq \
             -pe1 ./raw/SAMPLEA1_R1.fastq \
             -pe2 ./raw/SAMPLEA1_R2.fastq \
             --untrimmed-output ./output/processed/atropos/SAMPLEA1_R1.atropos_untrimmed.fastq \
             --untrimmed-paired-output ./output/processed/atropos/SAMPLEA1_R2.atropos_untrimmed.fastq \
               > ./output/processed/atropos/SAMPLEA1.atropos.log.out.txt \
              2> ./output/processed/atropos/SAMPLEA1.atropos.log.err.txt
```
Criar o diretório do pós-processamento
```python=
mkdir -p ./output/processed/fastqc/pre/ ./output/processed/fastqc/pos/ ./output/processed/atropos/
```
A saída do atropos no modo insert das *reads* que não foram podadas serão as entradas do atropos no modo adapter:
```python=
atropos trim    --aligner adapter \
                -e 0.1 \
                -n 2 \
                -m 1 \
                --match-read-wildcards \
                -O 3 \
                -q 20 \
                --pair-filter both \
                -pe1 ./output/processed/atropos/SAMPLEA1_R1.atropos_untrimmed.fastq \
                -pe2 ./output/processed/atropos/SAMPLEA1_R2.atropos_untrimmed.fastq \
                -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG \
                -g AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT  \
                -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
                -G CAAGCAGAAGACGGCATACGAGATNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT \
                -T 2 \
                -o ./output/processed/atropos/SAMPLEA1_R1.atropos_adapter.fastq  \
                -p ./output/processed/atropos/SAMPLEA1_R2.atropos_adapter.fastq \
                 > ./output/processed/atropos/SAMPLEA1.atropos_adapter.log.out.txt \
                2>  ./output/processed/atropos/SAMPLEA1.atropos_adapter.log.err.txt            
```

Concatenando os resultados de ambos os modos de execução:
```python=
cat     ./output/processed/atropos/SAMPLEA1_R1.atropos_insert.fastq \
        ./output/processed/atropos/SAMPLEA1_R1.atropos_adapter.fastq \
   > ./output/processed/atropos/SAMPLEA1_R1.atropos_final.fastq

cat     ./output/processed/atropos/SAMPLEA1_R2.atropos_insert.fastq \
        ./output/processed/atropos/SAMPLEA1_R2.atropos_adapter.fastq \
   > ./output/processed/atropos/SAMPLEA1_R2.atropos_final.fastq
```
Resultados obtidos após a execução:
![image](https://hackmd.io/_uploads/B172Xa2La.png)
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1_R1.atropos_final_fastqc.html

![image](https://hackmd.io/_uploads/ByiWNTn86.png)
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1_R2.atropos_final_fastqc.html

## Testes com prinseq-lite.pl
O prinseq-lite.pl é um script Perl que faz parte da ferramenta de análise de qualidade de sequências chamada PRINSEQ (Primer INdependent SEQuence analysis). É utilizado para realizar uma análise abrangente de dados de sequenciamento, especialmente para dados de sequenciamento de nova geração (NGS) como Illumina, 454, e outros. 

Tem como função específica a realizar a filtragem, limpeza e estatísticas básicas sobre dados de sequenciamento.

Trata-se de uma poda de qualidade baseada em janela deslizante (*sliding window trimming*), isto é, um processo de pré-processamento de dados de sequenciamento de DNA ou RNA que envolve a avaliação da qualidade das bases ao longo de uma leitura (sequência) utilizando uma janela deslizante. Esse método é frequentemente usado para eliminar partes de leituras que apresentam quedas significativas na qualidade das bases.

:::success
A qualidade de uma base em uma *read* é representada por um valor de qualidade associado a cada nucleotídeo (*Phred score*), que indica a probabilidade de erro na base correspondente. 

Quedas abruptas na qualidade ao longo de uma leitura podem sugerir a presença de áreas problemáticas, como regiões onde a qualidade do sequenciamento diminui. Isso também pode ocorrer em áreas próximas ao final da leitura, onde a qualidade frequentemente se reduz.
:::

```python=
mkdir -p ./output/processed/prinseq

prinseq-lite.pl -fastq ./raw/SAMPLEA1_R1.fastq \
                -fastq2 ./raw/SAMPLEA1_R2.fastq \
                -out_format 3 \
                -trim_qual_window 3 \
                -trim_qual_step 1 \
                -trim_qual_right 30 \
                -trim_qual_type mean \
                -trim_qual_rule lt \
                -out_good ./output/processed/prinseq/SAMPLEA1.prinseq \
                -out_bad ./output/processed/prinseq/SAMPLEA1.prinseq.bad \
                -lc_method dust \
                -lc_threshold 30 \
                -min_len 20 \
                -trim_tail_right 5 \
                -trim_tail_left 5\
                -ns_max_p 80 \
                -noniupac \
                > ./output/processed/prinseq/SAMPLEA1.prinseq.out.log \
                2> ./output/processed/prinseq/SAMPLEA1.prinseq.err.log
```
![image](https://hackmd.io/_uploads/S1eNwhpUp.png)
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1.prinseq_1_fastqc.html

![image](https://hackmd.io/_uploads/ryWIw2TUp.png)
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1.prinseq_2_fastqc.html

## Testes com fastqc (pós)
```python=
mkdir -p output/processed/fastqc/pos

fastqc -t 2 \
   ./output/processed/prinseq/SAMPLEA1.prinseq_1.fastq \
   -o ./output/processed/fastqc/pos/
fastqc -t 2 \
   ./output/processed/prinseq/SAMPLEA1_2.prinseq.fastq \
   -o ./output/processed/fastqc/pos/

## SE EXISTIR SAMPLEA1.prinseq_1_singletons.fastq
fastqc -t 2 \
   ./output/processed/prinseq/SAMPLEA1.prinseq_1_singletons.fastq \
   -o ./output/processed/fastqc/pos/
## SE EXISTIR SAMPLEA1.prinseq_2_singletons.fastq
fastqc -t 2 \
   ./output/processed/prinseq/SAMPLEA1.prinseq_2_singletons.fastq \
   -o ./output/processed/fastqc/pos/
```
Executando o script getTranscriptome-prok.sh
```python=
chmod a+x getTranscriptome-prok.sh
./getTranscriptome-prok.sh
```
Simulação 
```python=
runSim.sh . transcriptoma.fa . 1000
```
> Generaring simulation for group A ...
        Replicate 1
        Replicate 2
        Replicate 3
Generaring simulation for group B ...
        Replicate 1
        Replicate 2
        Replicate 3

Criar e entrar no diretório ```input/``` 
```python=
mkdir input
cd input/
```
Entrar no diretório **raw** e comprimir os arquivos fastq
```python=
cd ..
cd raw/
gzip *.fastq
```

Voltar para o diretório ```input/``` e copiar o arquivo **data.txt**
```python=
cd ..
cd input/
cp /data/bioaat/info/data.txt /data/bioaat/icampos/salmonella/input/
```
Executar createSymLink.sh no diretório ```salmonella/``` 
```python=
cd ../
createSymLink.sh ./input/data.txt ./raw/ ./input/
```
## Alinhamentos
Baixar o Integrative Genomics Viewer (IGV) para visualização dos alinhamentos

### Testes de alinhamento com Bowtie2
Uso eficiente da memória no alinhamento de sequências curtas (50 a 1000 pb) em genomas relativamente grandes (e.g. genoma humano = consumo de memória ~3,2GB). 

* É mais eficiente (rápido, mais eficiente e consome menos memória) geralmente para sequências maiores que 50 pb;
* Suporta alinhamento com GAPs;
* Alinhamento ponta-a-ponta ou local;
* Não há limite para tamanho das sequências;
* Sobrepõe com caracteres ambíguos (tabela IUPAC) (e.g. N) na referência;
* Score de alinhamento similar ao de Smith-Waterman e Needleman-Wunsch;
* Alinhamentos *paired-end* que não alinham de forma pareada podem ser alinhados de forma não pareada;
* Não alinha sequências em *colorspace*;

Obs.: Bowtie2 não considera que existe *splicing*, por isso é melhor para procarioto, já para eucarioto não é recomendado para RNA-seq

Criar o diretório index e indexação do genoma
```python=
mkdir index
cd index
bowtie2-build -f genome.fa genome
```
Alinhamento
```python=
cd ../
time bowtie2 -x ./refs/genome -q -1 ./raw/SAMPLEA1_R1.fastq -2 ./raw/SAMPLEA1_R2.fastq -S bowtie2-test.sam
```
Conversão de SAM para BAM
```python=
samtools view -b bowtie2-test.sam > bowtie2-test.bam
samtools sort bowtie2-test.bam -o bowtie2-test-sorted.bam
samtools index bowtie2-test-sorted.bam

```

Verificar o alinhamento de uma região
```python=
bowtie2 -c "AGAAAAAC" -x genome -S bowtie2-test.sam
```
:::info
Devem ser copiados os seguintes arquivos para serem avertos no IGV: 
bowtie2-test-sorted.bam
bowtie2-test-sorted.bam.bai
:::

Mapeamento de RNA-seq reads x Genoma referência:
![bowtie2](https://hackmd.io/_uploads/Syv7ack86.png)


### Alinhamento com Burrows-Wheeler Aligner (BWA)
O BWA é um pacote de software para mapear sequências de baixa divergência em relação a um grande genoma de referência, como o genoma humano. 

Possui como algortimos:
* BWA-backtrack: leitura de sequências Illumina de até 100 pb
* BWA-SW e BWA-MEM: leitura de sequências Illumina de 70bp to 1Mbp, com suporte de leituras longas e alinhamento dividido 

**Obs.:** O BWA-MEM é recomendado para consultas de alta qualidade, visto que é mais rápido e preciso. 

```python=
bwa index genome.fa
bwa mem genome.fa ../raw/SAMPLEA1_R1.fastq ../raw/SAMPLEA1_R2.fastq > bwa_SAMPLEA1.sam
```

### Comparar BWA e bowtie2
```python=
samtools sort -n -o bwa_SAMPLEA1_sorted.sam bwa_SAMPLEA1.sam
SAM_nameSorted_to_uniq_count_stats.pl bwa_SAMPLEA1_sorted.sam > bwa_SAMPLEA1_sorted_results.txt
vimdiff test_sorted_results.txt bwa_SAMPLEA1_sorted_results.txt 
```

### TopHat2
O TopHat é um programa que alinha leituras de RNA-Seq a um genoma para identificar *exon-exon splice junctions*. Foi desenvolvido com base no programa de mapeamento ultrarrápido de leituras curtas *Bowtie* e projetado para trabalhar com leituras produzidas pelo *Illumina Genome Analyzer*, embora os usuários tenham conseguido usar o TopHat com leituras de outras tecnologias. 

Para executar o TopHat2 devemos entrar no ambiente do python2
```python=
source /usr/local/bioinfo/anaconda3/etc/profile.d/conda.sh
conda activate python2
```
Utilizar o Samtools para converter um arquivo no formato BAM (Binary Alignment Map) para o formato SAM (Sequence Alignment/Map).

Ao usar o TopHat2, o arquivo SAM é um componente crucial do fluxo de trabalho de análise de RNA-Seq. Ele contém informações detalhadas sobre como as leituras se alinham ao genoma de referência, permitindo uma análise mais aprofundada e uma compreensão mais completa do transcritoma estudado.


```
cd tophat_out/
samtools view accepted_hits.bam > accepted_hits.sam
```

Executar TopHat2 de acordo com os parâmetros desejados utilizar o comando ```tophat2 --help``` quando for definir 
```python=
tophat2  --min-intron-length 500000 --max-intron-length 500000 --no-novel-juncs  genome ../raw/SAMPLEA1_R1.fastq ../raw/SAMPLEA1_R2.fastq
tophat2  --GTF ../ref/genome.gtf --min-intron-length 500000 --max-intron-length 500000 --no-novel-juncs  genome ../raw/SAMPLEA1_R1.fastq ../raw/SAMPLEA1_R2.fastq
tophat2 --b2-very-sensitive --transcriptome-index ./transcriptome --GTF ../ref/genome.gtf --min-intron-length 500000 --max-intron-length 500000 --no-novel-juncs  genome ../raw/SAMPLEA1_R1.fastq ../raw/SAMPLEA1_R2.fastq
tophat2 --b2-very-sensitive --library-type fr-firststrand  --transcriptome-index ./transcriptome --GTF ../ref/genome.gtf --min-intron-length 500000 --max-intron-length 500000 --no-novel-juncs  genome ../raw/SAMPLEA1_R1.fastq ../raw/SAMPLEA1_R2.fastq
```
### Spliced Transcripts Alignment to a Reference (STAR)
STAR é um algoritmo de alinhamento utilizado para mapear leituras de RNA ao genoma de referência. Foi desenvolvido para lidar eficientemente com dados de RNA-Seq e reconhecer junções de *splicing*, identificar transcritos e quantificar a expressão gênica.

Sincronizar o diretório ```salmonella/```  com o diretório ```salmonella_test/``` 
* ```-a``` sincroniza recursivamente os diretórios e mantém as permissões, datas e horários dos arquivos; 
* ```-v```exibe informações detalhadas sobre o processo de sincronização. 

```python=
rsync -av salmonella/ salmonella_test/
```
Criação do índice
```python=
STAR --runThreadN        2 \
     --runMode           genomeGenerate \
     --genomeFastaFiles  ./genome.fa \
     --genomeDir         ./star_index \
     --sjdbGTFfile       ./genome.gtf \
     --sjdbOverhang      149
```
Verificar se funcionou corretamente 
```python=
STAR --runThreadN        2      --runMode           genomeGenerate      --genomeFastaFiles  ./genome.fa      --genomeDir         ./star_index      --sjdbGTFfile       ./genome.gtf      --sjdbOverhang      149
        /usr/lib/rna-star/bin/STAR-avx2 --runThreadN 2 --runMode genomeGenerate --genomeFastaFiles ./genome.fa --genomeDir ./star_index --sjdbGTFfile ./genome.gtf --sjdbOverhang 149
```
Descompactar arquivos fastq do diretório ```../salmonella_test/raw/```
```python=
gzip -d *.fastq.gz
```
Executar o alinhamento 
```python=
STAR --runThreadN             2      --genomeDir      ./ref/star_index      --readFilesIn      ./raw/SAMPLEA1_R1.fastq ./raw/SAMPLEA1_R2.fastq      --sjdbGTFfile       ./ref/genome.gtf      --outFileNamePrefix ./star_out/
        /usr/lib/rna-star/bin/STAR-avx2 --runThreadN 2 --genomeDir ./ref/star_index --readFilesIn ./raw/SAMPLEA1_R1.fastq ./raw/SAMPLEA1_R2.fastq --sjdbGTFfile ./ref/genome.gtf --outFileNamePrefix ./star_out/
```

```
       Mapping speed, Million of reads per hour |       3.12
                          Number of input reads |       866
                      Average input read length |       302
                                    UNIQUE READS:
                   Uniquely mapped reads number |       833
                        Uniquely mapped reads % |       96.19%
                          Average mapped length |       301.09
                       Number of splices: Total |       0
            Number of splices: Annotated (sjdb) |       0
                       Number of splices: GT/AG |       0
                       Number of splices: GC/AG |       0
                       Number of splices: AT/AC |       0
               Number of splices: Non-canonical |       0
                      Mismatch rate per base, % |       0.72%
                         Deletion rate per base |       0.00%
                        Deletion average length |       0.00
                        Insertion rate per base |       0.00%
                       Insertion average length |       1.00
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       1
             % of reads mapped to multiple loci |       0.12%
        Number of reads mapped to too many loci |       0
             % of reads mapped to too many loci |       0.00%
                                  UNMAPPED READS:
  Number of reads unmapped: too many mismatches |       0
       % of reads unmapped: too many mismatches |       0.00%
            Number of reads unmapped: too short |       2
                 % of reads unmapped: too short |       0.23%
                Number of reads unmapped: other |       30
                     % of reads unmapped: other |       3.46%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
```

## Cufflinks
O conjunto de ferramentas Cufflinks pode ser usado para realizar vários tipos diferentes de análises para experimentos de RNA-Seq. O conjunto de Cufflinks inclui vários programas diferentes que trabalham juntos para executar essas análises.
Compreende os programas:
**1.Cuffcompare**: compara os conjuntos de anotações de transcritos de RNA e avalia a sobreposição entre eles; 

**2.Cuffmerge**: combina e unifica as montagnes de transcritos geradas por diferentes amostras individuais provenientes do experimento. O resultado final é um arquivo GTF unificado com a montagem consolidada; 

**3.Cuffquant**: realiza a quantificação da expressão gênica a partir dos dados de RNA-seq;

**4.Cuffdiff**: compara as abundâncias de transcritos e genes entre diferentes condições experimentais e identifica aquelas que são estatisticamente significativas 

**5.Cuffnorm**: normalização dos dados 

Executar o cufflinks
```python=
cufflinks ../index/tophat_out/accepted_hits.bam
cufflinks --GTF-guide ../ref/genome.gtf ../index/tophat_out/accepted_hits.bam
```

Para visualização
```python=
column -t -s$'\t'  genes.fpkm_tracking | less -S
column -t -s$'\t'  isoforms.fpkm_tracking | less -S
```


## Execução do pipeline rnaseq-ref-prok.sh
```python=
./rnaseq-ref-prok.sh star ./input ./output 3 ./ref/genome.gtf ./ref/genome.fa stringtie
```
Obs.: Informar a qualidade do cutoff 

Avaliar as etapas de todo processamento (se foi concluído com êxito)
```python=
cat ./output/star_index/STAR.genomeGenerate.log.*.txt
        /usr/lib/rna-star/bin/STAR-avx2 --runThreadN 3 --runMode genomeGenerate --genomeFastaFiles genome.fa --genomeDir ./ --sjdbGTFfile /srv/bioaat/icampos/salmonella/ref/genome.gtf --genomeSAindexNbases 10 --sjdbOverhang 149
```
Visualizar os valores de abundância
```python=
grep 'yes' gene_exp.diff | column -t -s$'\t'
```
```
MSTRG.1  MSTRG.1  tal        NC_011274.1:7554-8508    CNTRL  TREAT  OK  141809       301112       1.08635    4.39113   5e-05  5e-05  yes
MSTRG.2  MSTRG.2  satP       NC_011274.1:9265-9832    CNTRL  TREAT  OK  1.04437e+06  549715       -0.925875  -3.78783  5e-05  5e-05  yes
MSTRG.3  MSTRG.3  bcfB       NC_011274.1:25000-25687  CNTRL  TREAT  OK  1.03273e+06  409750       -1.33365   -6.43293  5e-05  5e-05  yes
MSTRG.4  MSTRG.4  nhaA       NC_011274.1:46074-47241  CNTRL  TREAT  OK  55149        155724       1.49758    4.44532   5e-05  5e-05  yes
MSTRG.5  MSTRG.5  ileS,lspA  NC_011274.1:53738-57073  CNTRL  TREAT  OK  1.05929e+06  2.25856e+06  1.0923     4.88835   5e-05  5e-05  yes
MSTRG.6  MSTRG.6  citS       NC_011274.1:66929-68268  CNTRL  TREAT  OK  133025       44686.1      -1.5738    -4.59758  5e-05  5e-05  yes
```
Obs.: Devido ao montador utilizado, os genes *ileS* e *lspA* foram agrupados. 


## Avaliação da expressão gênica diferencial 
Realizar a fusão dos dados resultantes do cuffdiff (gene_exp.diff) com descrições do gene obtidas de uma consulta do NCBI
```python=
esearch -db gene -query "550538[Taxonomy ID]" \
		| efetch -format tabular > ./ref/gene_result.txt
```

Selecionar apenas as colunas de interesse (3, 6 e 8),  respectivamente, GeneID, Symbol e description, para um novo arquivo ./refs/gene_info.txt.

```python=
cut -f 3,6,8 ./ref/gene_result.txt > ./ref/gene_info.txt
```
```
GeneID  Symbol  description
10021625  SGt58   tRNA-Pro
10021624  SGt35   tRNA-Val
10021623  SGt06   tRNA-Arg
10021622  SGt04   tRNA-Asp
10021621  SGt71   tRNA-Gly
10021620  SGt75   tRNA-Leu
10021619  SGt70   tRNA-Phe
10021618  SGt68   tRNA-Val
10021617  SGt66   tRNA-Sec
10021616  SGt48   tRNA-Leu
10021615  SGt26   tRNA-Ser
10021614  SGt17   tRNA-Ser
10021613  SGt13   tRNA-Lys
10021612  SGt53   tRNA-Thr
10021611  SGt45   tRNA-Gly
10021610  SGt30   tRNA-Asn
10021609  SGt29   tRNA-Asn
10021608  SGt27   tRNA-Ser
10021607  SGt25   tRNA-Arg
10021606  SGt24   tRNA-Arg
10021605  SGt21   tRNA-Tyr
10021604  SGt16   tRNA-Ser
10021603  SGt07   tRNA-Gln
10021602  SGt74   tRNA-Leu
10021601  SGt62   tRNA-Trp
10021600  SGt67   tRNA-Pro
10021599  SGt34   tRNA-Ala
10021598  SGt33   tRNA-Arg
10021597  SGt32   tRNA-Pro
10021596  SGt31   tRNA-Asn
10021595  SGt28   tRNA-Asn
10021594  SGt22   tRNA-Val
10021593  SGt05   tRNA-Thr
gene_info.txt
```

```python=
 splitteR.R  --x="./output/star_stringtie_cuffdiff/gene_exp.diff"             --col.x="gene"             --by.x=","             --out="./output/star_stringtie_cuffdiff/gene_exp.diff.splitted.txt"
```
```
test_id            gene_id            gene       locus                        sample_1  sample_2  status  value_1      value_2
MSTRG.1            MSTRG.1            tal        NC_011274.1:7554-8508        CNTRL     TREAT     OK      141809       301112
MSTRG.2            MSTRG.2            satP       NC_011274.1:9265-9832        CNTRL     TREAT     OK      1.04437e+06  549715
MSTRG.3            MSTRG.3            bcfB       NC_011274.1:25000-25687      CNTRL     TREAT     OK      1.03273e+06  409750
MSTRG.4            MSTRG.4            nhaA       NC_011274.1:46074-47241      CNTRL     TREAT     OK      55149        155724
MSTRG.5            MSTRG.5            ileS,lspA  NC_011274.1:53738-57073      CNTRL     TREAT     OK      1.05929e+06  2.25856e+06
MSTRG.6            MSTRG.6            citS       NC_011274.1:66929-68268 

log2(fold_change)  test_stat  p_value  q_value  significant
1.08635            4.39113    5e-05    5e-05    yes
-0.925875          -3.78783   5e-05    5e-05    yes
-1.33365           -6.43293   5e-05    5e-05    yes
1.49758            4.44532    5e-05    5e-05    yes
1.0923             4.88835    5e-05    5e-05    yes
-1.5738            -4.59758   5e-05    5e-05    yes
```

Obs.: Note que com o splitteR.R foi possível separar os genes *ileS* e *lspA* que estavam anteriormente agrupados.

Em seguida, o script mergeR.R realiza a fusão de ambas as tabelas geradas anteriormente

```python=
mergeR.R    --x="./output/star_stringtie_cuffdiff/gene_exp.diff.splitted.txt" \
            --all.x \
            --y="./ref/gene_info.txt" \
            --by.x="gene" \
            --by.y="Symbol" \
            --out="./output/star_stringtie_cuffdiff/gene_exp.diff.spplitted.merged.txt"
```

```
gene-SG_RS17105_1       gene-SG_RS17105_1       NC_011274.1:3544753-3545065     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS11470_1       gene-SG_RS11470_1       NC_011274.1:2332347-2337301     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS05715_1       gene-SG_RS05715_1       NC_011274.1:1208464-1209006     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS17135_1       gene-SG_RS17135_1       NC_011274.1:3548052-3549952     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS24495_1       gene-SG_RS24495_1       NC_011274.1:1279198-1280239     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS17145_1       gene-SG_RS17145_1       NC_011274.1:3550142-3550946     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS17150_1       gene-SG_RS17150_1       NC_011274.1:3550979-3551876     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS00030_1       gene-SG_RS00030_1       NC_011274.1:5855-7286   CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS11525_1       gene-SG_RS11525_1       NC_011274.1:2345699-2347343     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS00055_1       gene-SG_RS00055_1       NC_011274.1:10730-11135 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS17175_1       gene-SG_RS17175_1       NC_011274.1:3555999-3558036     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS05775_1       gene-SG_RS05775_1       NC_011274.1:1220310-1223728     CNTRL   TREAT   NOTEST  0       0       0       0       1         1       no      NA      NA
gene-SG_RS00070_1       gene-SG_RS00070_1       NC_011274.1:14905-15853 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS00075_1       gene-SG_RS00075_1       NC_011274.1:15979-16324 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS00080_1       gene-SG_RS00080_1       NC_011274.1:16384-16918 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS00085_1       gene-SG_RS00085_1       NC_011274.1:16934-17378 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA
gene-SG_RS00090_1       gene-SG_RS00090_1       NC_011274.1:17758-19855 CNTRL   TREAT   NOTEST  0       0       0       0       1       1         no      NA      NA

```
Entrar no ambiente R pelo navegador: http://200.145.102.33:8787/

Para a contrução do heatmap foi utilizado o seguinte script e obtenção de tabela com os principais resultados:
```python=
salmonella.df <- read.delim('/data/bioaat/icampos/salmonella/output/star_stringtie_cuffnorm/genes.fpkm_table', sep="\t")
head(salmonella.df)
colnames(salmonella.df)

samples.CNTRL <- colnames(salmonella.df)[ grep("^CNTRL\\_", colnames(salmonella.df)) ]
samples.TREAT <- colnames(salmonella.df)[ grep("^TREAT\\_", colnames(salmonella.df)) ]
salmonella.samps <- c(samples.CNTRL, samples.TREAT)

dim(salmonella.df)


genes_diff.df <-  read.delim('/data/bioaat/icampos/salmonella/output/star_stringtie_cuffdiff/gene_exp.diff',
                             sep='\t')


dim(genes_diff.df)

signif.genes_diff.df <- subset(genes_diff.df, q_value <= 0.05 & 
                                 abs(log2.fold_change.) >= log2(1.5) )

signif.salmonella.df <- merge(x=signif.genes_diff.df,
                              y=salmonella.df,
                              by.x='test_id',
                              by.y='tracking_id')

dim(signif.salmonella.df)

library('xlsx')

write.xlsx( signif.salmonella.df[, c('gene','sample_1','value_1',samples.CNTRL,
                                            'sample_2','value_2',samples.TREAT,
                                            'log2.fold_change.', 'q_value')],
            file="/data/bioaat/icampos/salmonella/signif_genes.xlsx",
            row.names=FALSE)


library('pheatmap')


dfh<-data.frame(group=as.character(salmonella.samps))

dfh$group <- ifelse( unlist(
  lapply( as.character(salmonella.samps) , function(x) { return(regexpr("^TREAT", x)[1]==1); } )
)
, "mutant", "wild-type" )

dfh$group <- as.factor(dfh$group) 
rownames(dfh) <- as.character(salmonella.samps)

rownames(signif.salmonella.df) <- signif.salmonella.df$gene

png(filename ="/data/bioaat/icampos/salmonella/salmonella_heatmap.png",
    res=300, width=3000, height=3000)

ann_colors <- list(
  group=c('mutant'='green', 'wild-type'='blue')
)

pheatmap(as.matrix(signif.salmonella.df[,salmonella.samps]),scale="row", 
         annotation_col = dfh,
         color=colorRampPalette(c("darkmagenta", "white", "gold"))(50),
         annotation_colors = ann_colors
         )

dev.off()
```
![salmonella_heatmap](https://hackmd.io/_uploads/ry1rQ4Jw6.png)


![image](https://hackmd.io/_uploads/B1Y0gcJva.png)

De acordo com os resultados obtidos a bactéria mutante apresentou uma maior expressão dos genes *tal*, *nhaA* e *ileS, lspA* quando comparada a estirpe selvagem, nesta os genes *satP*, *bcfB* e *citS* apresentou maior expressão. Esses resultados obtidos estão de acordo com os valores simulados nos arquivos de **abundance_A.txt** (CNTRL) e **abundance_B.txt** (TREAT).

## Montagem *de novo* de transcritomas
Nesse caso não há sequências que podem ser utilizadas como referências. Esse tipo de sequênciamento necessita de uma montagem das sequências com apenas os dados obtidos do sequênciamento. Envolve um processo de alinhamento entre as sequências geradas, que permite a obtenção de sequências consenso, os alinhamentos são analisados para a reconstrução dos transcritos. 

Diferentemente das montagens com referência, pode haver perda de informação biológica, por isso deve-se priorizar aquelas que possuem referência. 

Primeiramente, deve-se copiar os scripts rnaseq-denovo.sh e rnaseq-trinity-abundance.sh  para o diretório ```./salmonella```
```python=
cd /data/bioaat/icampos/salmonella
cp /data/bioaat/scripts/rnaseq-denovo.sh rnaseq-trinity-abundance.sh ./
```
Execução do script rnaseq-denovo.sh
```python=
./rnaseq-denovo.sh ./output/processed/cleaned/ ./output/ 8 20
```

Iremos indexar o transcritoma montado pelo trinity *vs* o transcrito real onde iniciamos a simulação dentro do diretório ```./raw```

**Obs.:** Antes da indexação, dentro do diretório ```./raw``` deve haver o arquivo transcriptoma.fa

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

No diretório ```./salmonella``` executar o blastn para realizar a comparação
```python=
blastn        -max_hsps 1 \
        -max_target_seqs 1 \
        -num_threads 8 \
        -query ./output/trinity_assembled/Trinity.fasta \
        -task blastn \
        -db ./raw/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
```
Observe que após executar o blastn obtemos os mesmos genes do nosso arquivo ACCS.txt
```
TRINITY_DN0_c0_g1_i1    nhaA    100.000 1148    0       0       1       1148    12      1159    0.0     2071
TRINITY_DN0_c0_g2_i1    satP    100.000 566     0       0       1       566     1       566     0.0     1021
TRINITY_DN0_c0_g3_i1    tal     100.000 952     0       0       1       952     2       953     0.0     1718
TRINITY_DN0_c0_g4_i1    citS    100.000 1318    0       0       1       1318    11      1328    0.0     2378
TRINITY_DN0_c1_g1_i1    bcfB    100.000 686     0       0       1       686     1       686     0.0     1238
TRINITY_DN0_c2_g1_i1    lspA    100.000 500     0       0       1       500     1       500     0.0     902
```

### Verificar a abundância
Copiar os scripts Detonate.sh, Transrate.sh,  run_transdecoder.sh e rnaseq-trinity-abundance.sh para o diretório ```./salmonella``` 

```python=
cp /data/bioaat/scripts/Detonate.sh Transrate.sh run_transdecoder.sh rnaseq-trinity-abundance.sh ./
```
#### Detonate.sh
O script Detonate.sh realiza o alinhamento das leituras utilizando o Bowtie2, verifica se o arquivo BAM resultante já existe antes de realizar o alinhamento, que será executado em paralelo para cada par de arquivos de leitura. Após o alinhamento, os arquivos BAM resultantes são fundidos em um único arquivo "All.aligned.bam". 

A ferramenta "rsem-eval-estimate-transcript-length-distribution" é utilizada para avaliar a distribuição dos tamanhos dos transcritos. 

O score DETONATE é calculado com "rsem-eval-calculate-score", considerando a distribuição de tamanho dos transcritos previamente avaliada. Trata-se de uma métrica importante para avaliar o quão bem a montagem *de novo* representa o transcritoma verdadeiro, auxiliando na validação da qualidade da montagem. 

Execução do Detonate.sh
```python=
./Detonate.sh ./output/processed/cleaned/ ./output/trinity_assembled/ ./raw/transcriptoma.fa ./output/
```
```
PE-detonate.score
Score   -239273.18
BIC_penalty     -29.71
Prior_score_on_contig_lengths   -42.59
Prior_score_on_contig_sequences -7167.14
Data_likelihood_in_log_space_without_correction -232039.22
Correction_term -5.48
Number_of_contigs       6
Expected_number_of_aligned_reads_given_the_data 4557.00
Number_of_contigs_smaller_than_expected_read/fragment_length    0
Number_of_contigs_with_no_read_aligned_to       0
Maximum_data_likelihood_in_log_space    -231971.38
Number_of_alignable_reads       4557
Number_of_alignments_in_total   4557
Transcript_length_distribution_related_factors  -38.50
```
```
PE-detonate.score.genes.results
gene_id transcript_id(s)        length  effective_length        expected_count  CPM     FPKM
TRINITY_DN0_c0_g1       TRINITY_DN0_c0_g1_i1    1148.00 999.00  524.00  78844.37        133477.96
TRINITY_DN0_c0_g2       TRINITY_DN0_c0_g2_i1    566.00  417.00  648.00  176646.23       299049.64
TRINITY_DN0_c0_g3       TRINITY_DN0_c0_g3_i1    952.00  803.00  813.00  144066.68       243894.76
TRINITY_DN0_c0_g4       TRINITY_DN0_c0_g4_i1    1318.00 1169.00 525.00  69853.16        118256.48
TRINITY_DN0_c1_g1       TRINITY_DN0_c1_g1_i1    686.00  537.00  1259.00 294000.99       497722.99
TRINITY_DN0_c2_g1       TRINITY_DN0_c2_g1_i1    500.00  351.00  788.00  236588.57       400527.80
```

#### Transrate.sh
O script Transrate.sh permite a avaliação da qualidade de montagem de transcriptomas.

Para executar o Transrate.sh é necessário um arquivo de proteoma referência; 
Criar arquivo do proteoma em ```./ref```:
```python=
gffread -g genome.fa genome.gtf -y protein.fa
```

Executar o Transrate.sh dentro do diretório ./salmonella:
```python=
Transrate.sh ./output/processed/cleaned/ ./output/trinity_assembled/Trinity.fasta ./ref/protein.fa ./output/
```
Visualizar o resultado o Transrate.sh
```python=
less ./output/Transrate.log.out.txt
```
```
[ INFO] 2023-12-14 17:41:43 : Loading assembly: /srv/bioaat/icampos/salmonella/output/trinity_assembled/Trinity.fasta
[ INFO] 2023-12-14 17:41:43 : Analysing assembly: /srv/bioaat/icampos/salmonella/output/trinity_assembled/Trinity.fasta
[ INFO] 2023-12-14 17:41:43 : Results will be saved in /srv/bioaat/icampos/salmonella/output/Transrate_good/Trinity
[ INFO] 2023-12-14 17:41:43 : Calculating contig metrics...
[ INFO] 2023-12-14 17:41:43 : Contig metrics:
[ INFO] 2023-12-14 17:41:43 : -----------------------------------
[ INFO] 2023-12-14 17:41:43 : n seqs                            6
[ INFO] 2023-12-14 17:41:43 : smallest                        500
[ INFO] 2023-12-14 17:41:43 : largest                        1318
[ INFO] 2023-12-14 17:41:43 : n bases                        5170
[ INFO] 2023-12-14 17:41:43 : mean len                     861.67
[ INFO] 2023-12-14 17:41:43 : n under 200                       0
[ INFO] 2023-12-14 17:41:43 : n over 1k                         2
[ INFO] 2023-12-14 17:41:43 : n over 10k                        0
[ INFO] 2023-12-14 17:41:43 : n with orf                        6
[ INFO] 2023-12-14 17:41:43 : mean orf percent              96.56
[ INFO] 2023-12-14 17:41:43 : n90                             566
[ INFO] 2023-12-14 17:41:43 : n70                             686
[ INFO] 2023-12-14 17:41:43 : n50                             952
[ INFO] 2023-12-14 17:41:43 : n30                            1148
[ INFO] 2023-12-14 17:41:43 : n10                            1318
[ INFO] 2023-12-14 17:41:43 : gc                             0.53
[ INFO] 2023-12-14 17:41:43 : bases n                           0
[ INFO] 2023-12-14 17:41:43 : proportion n                    0.0
[ INFO] 2023-12-14 17:41:43 : Contig metrics done in 0 seconds
[ INFO] 2023-12-14 17:41:43 : Calculating read diagnostics...
[ INFO] 2023-12-14 17:41:46 : Read mapping metrics:
[ INFO] 2023-12-14 17:41:46 : -----------------------------------
[ INFO] 2023-12-14 17:41:46 : fragments                      4860
[ INFO] 2023-12-14 17:41:46 : fragments mapped               4726
[ INFO] 2023-12-14 17:41:46 : p fragments mapped             0.97
[ INFO] 2023-12-14 17:41:46 : good mappings                  4694
[ INFO] 2023-12-14 17:41:46 : p good mapping                 0.97
[ INFO] 2023-12-14 17:41:46 : bad mappings                     32
[ INFO] 2023-12-14 17:41:46 : potential bridges                 0
[ INFO] 2023-12-14 17:41:46 : bases uncovered                   0
[ INFO] 2023-12-14 17:41:46 : p bases uncovered               0.0
[ INFO] 2023-12-14 17:41:46 : contigs uncovbase                 0
[ INFO] 2023-12-14 17:41:46 : p contigs uncovbase             0.0
[ INFO] 2023-12-14 17:41:46 : contigs uncovered                 0
[ INFO] 2023-12-14 17:41:46 : p contigs uncovered             0.0
[ INFO] 2023-12-14 17:41:46 : contigs lowcovered                0
[ INFO] 2023-12-14 17:41:46 : p contigs lowcovered            0.0
[ INFO] 2023-12-14 17:41:46 : contigs segmented                 0
[ INFO] 2023-12-14 17:41:46 : p contigs segmented             0.0
[ INFO] 2023-12-14 17:41:46 : Read metrics done in 3 seconds
[ INFO] 2023-12-14 17:41:46 : Calculating comparative metrics...
[ INFO] 2023-12-14 17:41:48 : Comparative metrics:
[ INFO] 2023-12-14 17:41:48 : -----------------------------------
[ INFO] 2023-12-14 17:41:48 : CRBB hits                         6
[ INFO] 2023-12-14 17:41:48 : n contigs with CRBB               6
[ INFO] 2023-12-14 17:41:48 : p contigs with CRBB             1.0
[ INFO] 2023-12-14 17:41:48 : rbh per reference               0.0
[ INFO] 2023-12-14 17:41:48 : n refs with CRBB                  6
[ INFO] 2023-12-14 17:41:48 : p refs with CRBB                0.0
[ INFO] 2023-12-14 17:41:48 : cov25                             6
[ INFO] 2023-12-14 17:41:48 : p cov25                         0.0
[ INFO] 2023-12-14 17:41:48 : cov50                             6
[ INFO] 2023-12-14 17:41:48 : p cov50                         0.0
[ INFO] 2023-12-14 17:41:48 : cov75                             6
[ INFO] 2023-12-14 17:41:48 : p cov75                         0.0
[ INFO] 2023-12-14 17:41:48 : cov85                             6
[ INFO] 2023-12-14 17:41:48 : p cov85                         0.0
[ INFO] 2023-12-14 17:41:48 : cov95                             6
[ INFO] 2023-12-14 17:41:48 : p cov95                         0.0
[ INFO] 2023-12-14 17:41:48 : reference coverage              0.0
[ INFO] 2023-12-14 17:41:48 : Comparative metrics done in 3 seconds
[ INFO] 2023-12-14 17:41:48 : -----------------------------------
[ INFO] 2023-12-14 17:41:48 : TRANSRATE ASSEMBLY SCORE     0.8735
[ INFO] 2023-12-14 17:41:48 : -----------------------------------
[ INFO] 2023-12-14 17:41:48 : TRANSRATE OPTIMAL SCORE      0.7856
[ INFO] 2023-12-14 17:41:48 : TRANSRATE OPTIMAL CUTOFF     0.8397
[ INFO] 2023-12-14 17:41:48 : good contigs                      6
[ INFO] 2023-12-14 17:41:48 : p good contigs                  1.0
[ INFO] 2023-12-14 17:41:48 : Writing contig metrics for each contig to /srv/bioaat/icampos/salmonella/output/Transrate_good/Trinity/contigs.csv
[ INFO] 2023-12-14 17:41:48 : Writing analysis results to assemblies.csv
```

#### run_transdecoder.sh 
O script run_transdecoder.sh realiza a anotação de ORFs em uma montagem de transcriptoma por meio da ferramenta TransDecoder e, em seguida, realiza a predição de codificação de sequências de proteínas usando informações de domínios Pfam e hits de busca BLASTp. 

Execução do run_transdecoder.sh dentro do diretório ```./salmonella```
```python=
./run_transdecoder.sh
```

#### rnaseq-trinity-abundance.sh
O script rnaseq-trinity-abundance.sh realiza várias etapas para analisar a abundância diferencial de genes em dados de RNA-Seq usando a montagem de transcriptoma obtida pelo Trinity.

Execução do rnaseq-trinity-abundance.sh dentro do diretório ```./salmonella```
```python=
rnaseq-trinity-abundance.sh ./output/renamed/ ./output/trinity_assembled/ ./output/
```
Visualização da matriz gerada
```python=
column -t -s$'\t' output/abundance/gene/DESeq2/abundance.gene.counts.matrix.CNTRL_vs_TREAT.DESeq2.DE_results | less -S
```
```
sampleA            sampleB  baseMeanA  baseMeanB         baseMean          log2FoldChange    lfcSE              stat               pvalue        >
TRINITY_DN0_c1_g1  CNTRL    TREAT      264.591476163054  107.789977079092  186.190726621073  1.29524397223327   0.115661958026464  11.19853056557>
TRINITY_DN0_c2_g1  CNTRL    TREAT      195.458479642819  414.143884665288  304.801182154054  -1.08193359853203  0.102008813066948  -10.6062757324>
TRINITY_DN0_c0_g2  CNTRL    TREAT      226.850640466564  114.504203130992  170.677421798778  0.983639958878017  0.11752979288567   8.369281819758>
TRINITY_DN0_c0_g3  CNTRL    TREAT      35.2010151060934  77.5262392962313  56.3636272011623  -1.13371239610826  0.180998126580625  -6.26366922976>
TRINITY_DN0_c0_g1  CNTRL    TREAT      13.0617402283097  38.0379471854087  25.5498437068592  -1.54454517414085  0.272120976804655  -5.67595042571>
TRINITY_DN0_c0_g4  CNTRL    TREAT      29.7961882026374  11.2406958707237  20.5184420366805  1.40996887586438   0.297719309258592  4.735899997133>
```

## Fluxograma de análises de Transcritomas
![Diagrama2](https://hackmd.io/_uploads/S1M2sEXva.jpg)

> Obs.: Os scripts citados foram abordados ao longo do relatório. 