# 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
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:
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1_R1.atropos_final_fastqc.html
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
```
file:///C:/Users/Usuario/Desktop/Transcriptoma/SAMPLEA1.prinseq_1_fastqc.html
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:
### 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()
```
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
> Obs.: Os scripts citados foram abordados ao longo do relatório.
Sign in with Wallet
Connect another wallet