# Bioinformática Aplicada II: Análise de Transcritomas
> [name=João Victor dos Anjos Almeida]
> [name=Docente: Professor Doutor Daniel Guariz Pinheiro]
---
* [Artigo utilizado](#Artigo-utilizado)
* [Contextualização](#Contextualização)
* [Metodologia seguida no artigo para análise dos transcritos](#Metodologia-seguida-no-artigo-para-análise-dos-transcritos)
* [Metodologia seguida pelo discente](#Metodologia-seguida-pelo-discente)
* [Organização de diretórios e obtenção dos dados](#Organização-de-diretórios-e-obtenção-dos-dados)
* [Pré-processamento](#Pré-processamento)
+ [FastQC](#FastQC)
+ [Atropos](#Atropos)
+ [PrinSeq](#PrinSeq)
* [Processamento](#Processamento)
+ [Alinhamento](#Alinhamento)
+ [Montagem](#Montagem)
+ [Análise dos dados de RNA-Seq](#Análise-de-dados-de-RNA-Seq)
* [Plot em R](#Plot-em-R)
- [Comparações com os Resultados obtidos](#Comparações-com-os-Resultados-obtidos)
<h2 id="artigo-utilizado">Artigo utilizado</h2>
WANG, Weizhe *et al*. Unraveling the mechanism of raffinose utilization in *Ligilactobacillus salivarius* Ren by transcriptomic analysis. 3 **Biotech**, v. 12, n. 9, p. 229, 2022. Doi: [10.1007/s13205-022-03280-6 ](https://doi.org/10.1007/s13205-022-03280-6)
<h2 id="contextualizacao">Contextualização</h2>
<div style="text-align: justify;">
O artigo em questão utiliza a espécie bacteriana *Ligilactobacillus salivarius* como organismo experimental, sendo esta amplamente reconhecida na literatura como um potente probiótico presente tanto no trato gastrointestinal quanto na cavidade oral de mamíferos. Essa espécie tem se destacado por suas propriedades benéficas à saúde, como a modulação da microbiota intestinal, o fortalecimento do sistema imunológico, a redução de inflamações gastrointestinais e a produção de substâncias antibacterianas, as bacteriocinas. O crescimento de *L. salivarius* pode ser estimulado por oligossacarídeos, como os galacto-oligossacarídeos e a inulina.
De acordo com o artigo e suas fontes, dentre os carboidratos dietéticos, a rafinose é um trissacarídeo de difícil digestão, mas que, em quantidades moderadas, exerce uma função prebiótica, favorecendo o crescimento de microrganismos benéficos. Bactérias ácido-láticas, como *L. plantarum*, *L. reuteri* e *Bifidobacterium lactis*, são capazes de degradar a rafinose por meio de duas principais vias: a levansacarase, que quebra a rafinose em frutose e melibiose, e a α-galactosidase, que converte a rafinose em sacarose e galactose. A capacidade de metabolizar a rafinose, no entanto, varia entre diferentes cepas de uma mesma espécie.
Nesse contexto, o estudo investigou o crescimento da cepa *L. salivarius* Ren em presença de rafinose e outros carboidratos, utilizando análise transcritômica (RNA-Seq) para identificar os genes envolvidos na captação e no metabolismo da rafinose.
</div>
<h2 id="Metodologiaart">Metodologia seguida no artigo para análise dos transcritos</h2>
<div style="text-align: justify;">
Para a extração de RNA, a cepa *L. salivarius* Ren foi inoculada (1% v/v) em 50 mL de mMRS, um caldo MRS modificado contendo os seguintes componentes: triptona (10 g/L), extrato de carne (10 g/L), extrato de levedura (5 g/L), Tween 80 (1 g/L), citrato de amônio (2 g/L), fosfato de potássio monobásico (3 g/L), fosfato de potássio dibásico (3 g/L), acetato de sódio (5 g/L), sulfato de magnésio (0,1 g/L) e sulfato de manganês (0,05 g/L), com a adição de 2% (p/v) de rafinose. Como controle, a linhagem foi cultivada em mMRS com 2% (p/v) de glicose ao invés de rafinose. Três replicados biológicos independentes foram realizados neste estudo. Os procedimentos, kits e equipamentos utilizados para o sequenciamento estão detalhados no artigo.
O sequenciamento foi realizado utilizando a plataforma Illumina Hiseq 2000 - como mencionado no texto do artigo; gerando leituras paired-end. Foram preferencialmente selecionados fragmentos de cDNA com comprimento entre 150 e 200 pb. Todas as leituras limpas foram alinhadas ao genoma de *L. salivarius* Ren (NZ_CP011403.1-NZ_CP011405.1) utilizando o [Bowtie2](https://github.com/BenLangmead/bowtie2) v2.0.6. A análise de genes diferencialmente expressos foi realizada com o software [DESeq](https://www.bioconductor.org/packages/2.11/bioc/html/DESeq.html) v1.10.1. Os valores de P resultantes foram ajustados pela abordagem de Benjamini e Hochberg para controle da taxa de falso positivo. Genes com valores de P ajustados inferiores a 0,05 foram considerados diferencialmente expressos.
Os dados brutos podem ser acessados no GEO com o número de acesso [GSE174336](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE184336).
</div>
## Metodologia seguida pelo discente
### Organização de diretórios e obtenção dos dados
<div style="text-align: justify;">
Primeiramente, foi preparado o diretório para a disciplina, criando-se os subdiretórios conforme a estrutura padrão apresentada ao longo das aulas. Essa organização visa facilitar o armazenamento e o acesso aos arquivos, garantindo um fluxo de trabalho eficiente.
```bash
cd /data/bioaat2024/joao/
mkdir L_salivarius
cd L_salivarius
mkdir input
mkdir output
mkdir refs
mkdir raw
`````
Os dados brutos foram obtidos a partir dos números de acesso das corridas depositadas no banco SRA do NCBI. Para isso, foi seguido o fluxo de informações disponível no site, começando pela identificação do BioProject do estudo ([PRJNA729477](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA729477/)) e os respectivos BioSamples para as três réplicas tratadas com suplementação de rafinose no meio (LSRR) - [SAMN19117227](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117227), [SAMN19117226](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117226) e [SAMN19117225](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117225) - e para o controle (LSRG) com glicose - [SAMN19117231](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117231), [SAMN19117230](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117230) e [SAMN19117229](https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN19117229). Os identificadores das corridas foram os seguintes:
</div>
```diff!
SRR14510976 LSRR1
SRR14510977 LSRR2
SRR14510978 LSRR3
SRR14510973 LSRG1
SRR14510974 LSRG2
SRR14510975 LSRG3
```````
<div style="background-color: #ffcccc; border-left: 5px solid red; padding: 10px; color: black; font-weight: bold;">
Apesar do texto do artigo indicar o uso da plataforma Illumina Hiseq 2000, tanto no GEO como no repositório SRA do NCBI, é indicada a plataforma Illumina Hiseq 2500.
</div>
Para obter os dados brutos, utilizou-se o seguinte comando, baseado em uma lista contendo os identificadores de cada arquivo a ser baixado:
```bash
cd /data/bioaat2024/joao/L_salivarius/raw
nano lista_sra.txt
for i in `cat list_sra.txt`; do fastq-dump --split-3 $i; done
`````
<div style="text-align: justify;">
Para facilitar a manipulação dos dados brutos gerados e garantir maior segurança no processo de utilização das leituras, foram criados links simbólicos no diretório de entrada (input). Esses links, direcionados aos arquivos correspondentes. A criação dos links simbólicos permite o acesso rápido aos arquivos sem duplicação de dados, o que é especialmente importante em cenários reais de processamento de grandes volumes de dados.
</div>
```bash
cd /data/bioaat2024/joao/L_salivarius/input
nano mklink.sh
#### mklink.sh
#!/bin/bash
ln -s ../raw/SRR14510976_1.fastq LSRR1_R1.fastq
ln -s ../raw/SRR14510976_2.fastq LSRR1_R2.fastq
ln -s ../raw/SRR14510977_1.fastq LSRR2_R1.fastq
ln -s ../raw/SRR14510977_2.fastq LSRR2_R2.fastq
ln -s ../raw/SRR14510978_1.fastq LSRR3_R1.fastq
ln -s ../raw/SRR14510978_2.fastq LSRR3_R2.fastq
ln -s ../raw/SRR14510973_1.fastq LSRG1_R1.fastq
ln -s ../raw/SRR14510973_2.fastq LSRG1_R2.fastq
ln -s ../raw/SRR14510974_1.fastq LSRG2_R1.fastq
ln -s ../raw/SRR14510974_2.fastq LSRG2_R2.fastq
ln -s ../raw/SRR14510975_1.fastq LSRG3_R1.fastq
ln -s ../raw/SRR14510975_2.fastq LSRG3_R2.fastq
####
chmod a+x mklink.sh
./mklink.sh
```
<div style="text-align: justify;">
Para obter os dados de referência da espécie utilizada no estudo, foi consultado o repositório FTP RefSeq do NCBI, com o identificador da montagem de *L. salivarius Ren* [(GCF_001011095.1)](https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/011/095/GCF_001011095.1_ASM101109v1/). Dessa forma, foi possível baixar os arquivos contendo a sequência genômica (genomic.fna), a anotação (genomic.gtf) e as características (feature_table.txt), que foram essenciais para a descrição e análise dos transcritos montados no estudo. Esses dados de referência são fundamentais para a correta interpretação e comparação das leituras obtidas com os transcritos da cepa analisada.
</div>
```bash
cd /data/bioaat2024/joao/L_salivarius/refs
#sequência genômica
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/011/095/GCF_001011095.1_ASM101109v1/GCF_001011095.1_ASM101109v1_genomic.fna.gz
#anotação
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/011/095/GCF_001011095.1_ASM101109v1/GCF_001011095.1_ASM101109v1_genomic.gtf.gz
#características
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/011/095/GCF_001011095.1_ASM101109v1/GCF_001011095.1_ASM101109v1_feature_count.txt
```
<div style="text-align: justify;">
<h3 id="pre-processamento">Pré-processamento</h3>
Através desses arquivos, foi possível realizar a execução inicial do primeiro script disponibilizado pelo docente e adaptado em aula, com o objetivo de pré-processar os dados brutos obtidos.
O script começa com a definição de argumentos e a validação desses parâmetros, garantindo que os comandos sejam executados corretamente. Além disso, cria diretórios específicos para armazenar as saídas dos processos realizados. Os programas utilizados nesta etapa incluem com as respectivas versões durante o período de execução da disciplina: [FastQC](https://github.com/s-andrews/FastQC) v0.12.1, [Atropos](https://github.com/jdidion/atropos) v1.1.32 e [Prinseq](https://github.com/uwb-linux/prinseq) v0.20.4.
#### FastQC
O primeiro programa a ser executado, para cada par de arquivos R1 e R2, é o FastQC, que realiza uma análise de qualidade inicial das leituras. O script gera arquivos de log e relatórios detalhados sobre a qualidade dos dados brutos. Este procedimento é realizado tanto antes quanto após o pré-processamento dos dados, permitindo comparar a qualidade das leituras antes e depois das etapas de filtragem e correção.
```bash
# RECORTE DO SCRIPT COM A ETAPA DE EXECUÇÃO DO FASTQC
name=`basename ${r1} | sed 's/_R1.fastq\(\.gz\)\?//'`
echo -e "Processing library ${name} ..."
if [ ! -e "${outdir}/processed/fastqc/pre/${name}_R1_fastqc.html" ] ||
[ ! -s "${outdir}/processed/fastqc/pre/${name}_R1_fastqc.html" ]; then
echo -e "\tFastQC evaluation using sample ${name}: ${r1} ...\n"
fastqc -t ${THREADS} \
${r1} \
-o ${outdir}/processed/fastqc/pre/ \
1> ${outdir}/processed/fastqc/pre/${name}_R1.log.out.txt \
2> ${outdir}/processed/fastqc/pre/${name}_R1.log.err.txt
fi
if [ ! -e "${outdir}/processed/fastqc/pre/${name}_R2_fastqc.html" ] ||
[ ! -s "${outdir}/processed/fastqc/pre/${name}_R2_fastqc.html" ]; then
echo -e "\tFastQC evaluation using sample ${name}: ${r2} ...\n"
fastqc -t ${THREADS} \
${r2} \
-o ${outdir}/processed/fastqc/pre/ \
1> ${outdir}/processed/fastqc/pre/${name}_R2.log.out.txt \
2> ${outdir}/processed/fastqc/pre/${name}_R2.log.err.txt
fi
```
Ao decorrer desse processo, os arquivos de saída gerados são submetidos ao programa [MultiQC](https://github.com/MultiQC/MultiQC) v1.25.1, que realiza uma visualização integrada das métricas de qualidade. O MultiQC compila os relatórios do FastQC e outros programas em uma única interface gráfica, facilitando a análise comparativa e permitindo uma avaliação geral da qualidade dos dados antes e após o pré-processamento.
```bash
if [ ! -e ${outdir}/processed/fastqc/pre/multiqc_report.html ]; then
echo "Running multiqc (pre) ..."
multiqc --force \
--outdir ${outdir}/processed/fastqc/pre/ \
${outdir}/processed/fastqc/pre/ \
> ${outdir}/processed/fastqc/pre/multiqc.log.out.txt \
2> ${outdir}/processed/fastqc/pre/multiqc.log.err.txt
fi
if [ ! -e ${outdir}/processed/fastqc/pos/multiqc_report.html ]; then
echo "Running multiqc (pos) ..."
multiqc --force \
--outdir ${outdir}/processed/fastqc/pos/ \
${outdir}/processed/fastqc/pos/ \
> ${outdir}/processed/fastqc/pos/multiqc.log.out.txt \
2> ${outdir}/processed/fastqc/pos/multiqc.log.err.txt
fi
```
#### Atropos
Em seguida a checagem inicial, o script prossegue com a remoção de adaptadores utilizando o Atropos. Este programa aplica uma série de filtros nos arquivos R1 e R2, como taxa de erro, comprimento mínimo das leituras, entre outros critérios. O objetivo é garantir que os dados processados sejam mais limpos e precisos, removendo sequências indesejadas ou de baixa qualidade antes da análise posterior.
```bash
# RECORTE DO SCRIPT COM A ETAPA DE EXECUÇÃO DO ATROPOS
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 AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTGAAAAA \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTNNNNNNNNNNNNNNNNNNNNGTGGTCGCCGTATCATTAAAAAA \
-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 \
2> ${outdir}/processed/atropos/${name}.atropos.log.err.txt
echo -e "\tRunning atropos (adapter) for adapter trimming using sample ${name}: ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq & ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq ...\n"
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 AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTGAAAAA \
-g TTTTTTAATGATACGGCGACCACNNNNNNNNNNNNNNNNNNNNACACTCTTTCCCTACACGACGCTCTTCCGATCT \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTNNNNNNNNNNNNNNNNNNNNGTGGTCGCCGTATCATTAAAAAA \
-G TTTTTCAAGCAGAAGACGGCATACGAGATNNNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT \
-T ${THREADS} \
-o ${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
-p ${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}.atropos_adapter.log.out.txt \
2> ${outdir}/processed/atropos/${name}.atropos_adapter.log.err.txt
echo -e "\tMerging atropos adapter trimming results using sample ${name}: ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq and ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq + ${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq and ${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq ...\n"
cat ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq
cat ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq
echo -e "\tRemoving useless atropos results ...\n"
rm -f ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq
fi
```
#### PrinSeq
Por fim, o script utiliza o PrinSeq para limpar as leituras que passaram pelo processo de Atropos, removendo adaptadores e outras sequências de baixa qualidade. Essas etapas têm como objetivo aprimorar ainda mais a qualidade das leituras, garantindo que os dados estejam adequados para análises subsequentes, com maior precisão e confiabilidade.
```bash
# RECORTE DO SCRIPT COM A ETAPA DE EXECUÇÃO DO PRINSEQ
if [ ! -e "${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1.fastq" ] ||
[ ! -s "${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1.fastq" ] ; then
echo -e "\tPrinSeq processing: ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq & ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq ...\n"
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 50 \
-trim_tail_right 5 \
-trim_tail_left 5 \
-ns_max_p 80 \
-min_qual_mean 30 \
> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.out.log \
2> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.err.log
fi
```
Todo o script então foi executado com:
```bash
cd /data/bioaat2024/joao/L_salivarius
./rnaseq-preprocess.sh input output 10
```
<div style="background-color: #ffcccc; border-left: 5px solid yellow; padding: 10px; color: black; font-weight: bold;">
O script também possibilita a utilização um arquivo de contaminantes, quando oferecido o script usa o Bowtie2 para alinhar as sequências contra os contaminantes e excluir as leituras que correspondem a esses contaminantes. Contudo, essa opção não foi utilizada em nosso roteiro.
</div>
A ferramenta [bamPEFragmentSize](https://deeptools.readthedocs.io/en/develop/content/installation.html) v3.5.5 foi utilizada para verificar os tamanhos dos fragmentos após a filtragem a partir dos arquivos BAM de sequenciamento paired-end gerados.
```bash
cd /data/bioaat2024/joao/L_salivarius/output/align_out_final/
bamPEFragmentSize -hist fragmentSize.png -b LSRG1/Aligned.sorted.bam LSRG2/Aligned.sorted.bam LSRG3/Aligned.sorted.bam LSRR1/Aligned.sorted.bam LSRR2/Aligned.sorted.bam LSRR3/Aligned.sorted.bam --samplesLabel LSRG1 LSRG2 LSRG3 LSRR1 LSRR2 LSRR3
```
### Processamento
Com os dados pré-processados, foi iniciado o processamento das sequências limpas por meio de outro script discutido em aula, específico para genomas de procariotos. Assim como no script anterior, há uma etapa de preparação dos diretórios e manipulação dos dados para a execução do processo. O script verifica se os diretórios de entrada e saída foram corretamente passados como parâmetros na linha de comando. Além disso, caso o procedimento anterior não tenha sido executado de forma independente, há a possibilidade de fornecer o caminho para o script anterior, permitindo que ambos os processos sejam executados de maneira conjunta.
#### Alinhamento
O script, para o desenvolvimento dessa etapa,permite a utilização de três alinhadores diferentes: [Tophat2](https://ccb.jhu.edu/software/tophat/index.shtml) v2.1.1, [HISAT2](https://daehwankimlab.github.io/hisat2/) v2.2.1 e [STAR](https://github.com/alexdobin/STAR) v1.5.3. Para a execução do alinhamento nesta atividade, foi escolhido o HISAT2 (Hierarchical Indexing for Spliced Transcript Alignment) sendo um alinhador de alta performance desenvolvido para mapear sequências de RNA-Seq em genomas de referência, e já utilizado durante as aulas.
```bash
# RECORTE DO SCRIPT COM A ETAPA DE PROCESSOS LIGADOS A ESCOLHA DO HISAT2
elif [ "${aligner}" == "hisat" ]; then
DTAPARAM=""
if [ "${assembler}" == "cufflinks" ]; then
DTAPARAM="--dta-cufflinks"
else
DTAPARAM="--dta"
fi
if [ ! -e "${outdir}/hisat_out_pe/${name}/accepted_hits.bam" ]; then
mkdir -p ${outdir}/hisat_out_pe/${name}
echo -e "\tHISAT2 alignment (${name}) paired-end reads X genome ..."
hisat2 ${DTAPARAM} \
--no-spliced-alignment \
--all \
--minins 300 \
--maxins 800 \
--threads ${THREADS} \
--rna-strandness RF \
-S ${outdir}/hisat_out_pe/${name}/accepted_hits.sam \
-x ${outdir}/hisat_index/genome \
-1 ${r1} \
-2 ${r2} \
> ${outdir}/hisat_out_pe/${name}.log.out.txt \
2> ${outdir}/hisat_out_pe/${name}.log.err.txt
samtools view -b --threads ${THREADS} --output ${outdir}/hisat_out_pe/${name}/accepted_hits.bam ${outdir}/hisat_out_pe/${name}/accepted_hits.sam
rm -f ${outdir}/hisat_out_pe/${name}/accepted_hits.sam
else
echo -e "\tFound HISAT2 output for PE (${name})..."
fi
if [ ! -e "${outdir}/hisat_out_se/${name}/accepted_hits.bam" ]; then
mkdir -p ${outdir}/hisat_out_se/${name}
if [ -s "${singletons}" ]; then
echo -e "\tHISAT2 alignment (${name}) singleton reads X genome ..."
hisat2 ${DTAPARAM} \
--no-spliced-alignment \
--all \
--no-spliced-alignment \
--threads ${THREADS} \
-S ${outdir}/hisat_out_se/${name}/accepted_hits.sam \
-x ${outdir}/hisat_index/genome \
-U ${singletons} \
> ${outdir}/hisat_out_se/${name}.log.out.txt \
2> ${outdir}/hisat_out_se/${name}.log.err.txt
samtools view -b --threads ${THREADS} --output ${outdir}/hisat_out_se/${name}/accepted_hits.bam ${outdir}/hisat_out_se/${name}/accepted_hits.sam
rm -f ${outdir}/hisat_out_se/${name}/accepted_hits.sam
fi
else
echo -e "\tFound HISAT2 output for SE (${name})..."
fi
if [ ! -e "${outdir}/hisat_out_final/${name}/accepted_hits.bam" ]; then
mkdir -p ${outdir}/hisat_out_final/${name}
if [ -s "${outdir}/hisat_out_pe/${name}/accepted_hits.bam" ]; then
if [ -s "${outdir}/hisat_out_se/${name}/accepted_hits.bam" ]; then
echo -e "\tMerging HISAT2 results ..."
samtools view -H ${outdir}/hisat_out_pe/${name}/accepted_hits.bam > ${outdir}/hisat_out_final/${name}/Header.txt
samtools merge -n --threads ${THREADS} \
-h ${outdir}/hisat_out_final/${name}/Header.txt \
${outdir}/hisat_out_final/${name}/accepted_hits.bam \
${outdir}/hisat_out_pe/${name}/accepted_hits.bam \
${outdir}/hisat_out_se/${name}/accepted_hits.bam \
> ${outdir}/hisat_out_final/${name}.log.out.txt \
2> ${outdir}/hisat_out_final/${name}.log.err.txt
else
pe_result_abs_path=$(readlink -f ${outdir}/hisat_out_pe/${name}/accepted_hits.bam)
cd ${outdir}/hisat_out_final/${name}/
ln -s ${pe_result_abs_path} accepted_hits.bam
cd ${curdir}
fi
else
if [ -s "${outdir}/hisat_out_se/${name}/accepted_hits.bam" ]; then
se_result_abs_path=$(readlink -f ${outdir}/hisat_out_se/${name}/accepted_hits.bam)
cd ${outdir}/hisat_out_final/${name}/
ln -s ${se_result_abs_path} accepted_hits.bam
cd ${curdir}
else
echo -e "[ERROR] Not found any alignment for PE or SE reads." 1>&2
fi
fi
else
echo -e "\tFound HISAT2 output final (${name})..."
fi
if [ ! -e "${outdir}/hisat_out_final/${name}/Aligned.sorted.bam" ]; then
echo -e "\tSorting alignments (${name})..."
samtools sort --threads ${THREADS} \
-o ${outdir}/hisat_out_final/${name}/Aligned.sorted.bam \
${outdir}/hisat_out_final/${name}/accepted_hits.bam \
> ${outdir}/hisat_out_final/${name}/Aligned.sorted.out.txt \
2> ${outdir}/hisat_out_final/${name}/Aligned.sorted.err.txt
fi
# SEMPRE VAMOS REMOVER O LINK SIMBÓLICO PARA QUE AO ESCOLHER UM OUTRO
# ALINHADOR ELE SEJA SUBSTITUÍDO
#if [ ! -e "${outdir}/align_out_final/${name}/Aligned.out.bam" ]; then
rm -f ${outdir}/align_out_final/${name}/Aligned.out.bam
rm -f ${outdir}/align_out_final/${name}/Aligned.sorted.bam
if [ -e "${outdir}/hisat_out_final/${name}/accepted_hits.bam" ]; then
align_final_out=`readlink -f ${outdir}/hisat_out_final/${name}/accepted_hits.bam`
align_sorted_out=`readlink -f ${outdir}/hisat_out_final/${name}/Aligned.sorted.bam`
cd ${outdir}/align_out_final/${name}
ln -s ${align_final_out} Aligned.out.bam
ln -s ${align_sorted_out} Aligned.sorted.bam
cd ${curdir}
else
echo "[ERROR] Not found HISAT2 final output (${outdir}/hisat_out_final/${name}/accepted_hits.bam)" 2>&1
exit
fi
```
#### Montagem
O script oferece duas opções para a montagem dos transcritos: [Cufflinks](https://cole-trapnell-lab.github.io/cufflinks/) v2.2.1 ou [StringTie](https://github.com/gpertea/stringtie) v2.2.3. Para a execução desta atividade, foi escolhido o programa Cufflinks.
Cufflinks é uma ferramenta amplamente utilizada na análise de dados de RNA-Seq para a montagem de transcritos e quantificação de expressão gênica. Como na opção pelo HISAT2, o Cufflinks também foi um dos programas explorados durante a disciplina, propiciando a sua para execução. A primeira ferramenta ligada ao conjunto utilizada após a seleção do alinhador como Cufflinks foi o pacote, responsável por combinar múltiplos arquivos de montagem de transcritos em uma única montagem unificada.
```bash
# RECORTE DO SCRIPT COM A ETAPA DE PROCESSOS LIGADOS A ESCOLHA DO CUFFLINKS
if [ "${assembler}" == "cufflinks" ]; then
if [ ! -e "${outdir}/${aligner}_cufflinks/${name}/transcripts.gtf" ]; then
echo -e "\tAssembly transcriptome (${name}) [cufflinks]\n"
cufflinks --output-dir ${outdir}/${aligner}_cufflinks/${name} \
--num-threads ${THREADS} \
--GTF ${refgtf} \
--frag-bias-correct ${refseq} \
--multi-read-correct \
--library-type fr-firststrand \
--frag-len-mean 500 \
--frag-len-std-dev 200 \
--total-hits-norm \
--min-isoform-fraction 0.5 \
--pre-mrna-fraction 0.5 \
--min-frags-per-transfrag 10 \
--junc-alpha 0.0000009 \
--small-anchor-fraction 0.5 \
--overhang-tolerance 0 \
--min-intron-length 5000000 \
--max-intron-length 1 \
--trim-3-avgcov-thresh 0.05 \
--trim-3-dropoff-frac 0.01 \
--max-multiread-fraction 0.75 \
--overlap-radius 1 \
--3-overhang-tolerance 0 \
--intron-overhang-tolerance 0 \
${outdir}/align_out_final/${name}/Aligned.sorted.bam \
> ${outdir}/${aligner}_cufflinks/${name}/cufflinks.out.txt \
2> ${outdir}/${aligner}_cufflinks/${name}/cufflinks.err.txt
fi
...
if [ "${assembler}" == "cufflinks" ]; then
if [ ! -e "${outdir}/${aligner}_cuffmerge/merged.gtf" ]; then
echo -e "\tMerging transcriptomes (${outdir}/assembly_GTF_list.txt) in a transcriptome reference [cuffmerge]\n"
cuffmerge -o ${outdir}/${aligner}_cuffmerge \
--ref-gtf ${refgtf} \
--ref-sequence ${refseq} \
--min-isoform-fraction 0.5 \
--num-threads ${THREADS} \
${outdir}/assembly_GTF_list.txt \
> ${outdir}/${aligner}_cuffmerge/cuffmerge.out.txt \
2> ${outdir}/${aligner}_cuffmerge/cuffmerge.err.txt
fi
```
#### Análise de dados de RNA-Seq
Com a escolha do Cufflinks como montador, o processo de análise envolveu uma série de ferramentas adicionais para interpretação dos dados gerados. Após a montagem inicial dos transcritos com o Cufflinks, é possível usar outras ferramentas da [Cufflinks suite](https://cole-trapnell-lab.github.io/cufflinks/manual/) para aprofundar a análise de expressão gênica e detectar diferenças entre condições experimentais.
O primeiro passo após a montagem, seguindo o script fornecido é o uso do Cuffcompare, que compara os transcritos montados com as sequências anotadas no genoma de referência. Essa comparação permite validar a precisão da montagem.
Depois de realizar essa comparação, o próximo passo é a quantificação da expressão dos transcritos, que é feita com o Cuffquant. Essa ferramenta calcula a abundância relativa de cada transcrito nas diferentes amostras, gerando um arquivo com a quantificação das expressões.
Com as quantificações de expressão em mãos, o Cuffnorm entra em ação para normalizar os dados. A normalização ajuda a ajustar as quantificações para que as diferenças observadas entre as amostras sejam de fato biológicas, e não causadas por artefatos técnicos.
O último passo do processo é a análise de expressão diferencial, que é realizada com o Cuffdiff. Essa ferramenta compara os níveis de expressão dos transcritos entre diferentes condições experimentais e identifica quais genes ou isoformas são diferencialmente expressos. O Cuffdiff gera um relatório com valores de p ajustados, o que permite identificar com confiança os genes que estão sendo regulados de forma significativa nas condições testadas.
```bash
#inicio cuffcompare
if [ ! -e "${outdir}/${aligner}_${assembler}_cuffcompare/cuffcmp.combined.gtf" ]; then
echo -e "\tRunning cuffcompare with ${aligner} & ${assembler} transcriptome reference (${transcriptomeref})..."
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
fi
#fim cuffcompare
#inicio cuffquant
# LISTA DE VALORES NÃO REDUNDANTES (NOME DO GRUPO BIOLÓGICO)
# Ex.: (CONTROL TEST)
biogroup_label=()
for bamfile in `find ${outdir}/align_out_final -name Aligned.sorted.bam`; do
name=`basename $(dirname ${bamfile})`
if [ ! -e "${outdir}/${aligner}_${assembler}_cuffquant/${name}/abundances.cxb" ]; then
echo -e "\tRunning cuffquant using sample ${name} as using ${aligner} & ${assembler} (${transcriptomeref}) ..."
mkdir -p ${outdir}/${aligner}_${assembler}_cuffquant/${name}
cuffquant --output-dir ${outdir}/${aligner}_${assembler}_cuffquant/${name} \
--frag-bias-correct ${refseq} \
--multi-read-correct \
--num-threads ${THREADS} \
--library-type fr-firststrand \
--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
fi
groupname=`echo ${name} | sed 's/[0-9]\+$//'`
biogroup_label=($(printf "%s\n" ${biogroup_label[@]} ${groupname} | sort -u ))
done
#fim cuffquant
biogroup_files=()
echo -e "\tCollecting Expression Data from cuffquant output (*.cxb) ..."
for label in ${biogroup_label[@]}; do
echo -e "\t\tCollecting .cxb files for ${label} ..."
group=()
for cxbfile in `ls ${outdir}/${aligner}_${assembler}_cuffquant/${label}*/abundances.cxb`; do
echo -e "\t\t\tFound ${cxbfile}"
group=(${group[@]} "${cxbfile}")
done
biogroup_files=(${biogroup_files[@]} $(IFS=, ; echo "${group[*]}") )
done
echo -e "Starting Gene Expression Analysis ..."
echo -e "\t\tLabels.: " $(IFS=, ; echo "${biogroup_label[*]}")
echo -e "\t\tFiles..: " ${biogroup_files[*]}
if [ ! -e "${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.count_table" ]; then
echo -e "\t\t\tGenerating abundance matrices (cuffnorm) ..."
cuffnorm --output-dir ${outdir}/${aligner}_${assembler}_cuffnorm \
--labels $(IFS=, ; echo "${biogroup_label[*]}") \
--num-threads ${THREADS} \
--library-type fr-firststrand \
--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
fi
if [ ! -e "${outdir}/${aligner}_${assembler}_cuffnorm/isoforms.raw_count_table.txt" ]; then
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.out.txt \
2> ${outdir}/${aligner}_${assembler}_cuffnorm/de-normalize-cuffnorm.err.txt
fi
if [ ! -e "${outdir}/${aligner}_${assembler}_cuffdiff/isoform_exp.diff" ]; then
echo -e "\t\t\tAnalysing differential expression (cuffdiff) ..."
```
Todo o script foi executado com:
```bash
cd /data/bioaat2024/joao/L_salivarius
./rnaseq-ref-prok.sh hisat input/ output/ 10 refs/genome_fixed.gtf refs/GCF_001011095.1_ASM101109v1_genomic.fna cufflinks
```
<div style="background-color: #ffcccc; border-left: 5px solid yellow; padding: 10px; color: black; font-weight: bold;">
Como no pré-processamento, nesse caso também há a possibilidade de se utilizar um arquivo de contaminantes, onde recebendo um arquivo de contaminantes, o script cria um índice para o Bowtie2 usando esse arquivo. Contudo, essa opção não foi utilizada em nosso roteiro.
</div>
### Plot em R
Com os resultados dos cálculos de expressão diferencial, foram utilizados scripts em R para gerar um heatmap que compare a expressão gênica entre os dois grupos analisados na pesquisa, auxiliando na visualização dos dados obtidos com a realização da atividade proposta.
O arquivo description.txt foi construído com base nos Ids genéricos e associados presentes nos arquivos resultantes da montagem com Cufflinks e processamento com demais programas do suite, e com as descrições presentes no arquivo de caracterśiticas baixados anteriormente para o genoma de referência.
```R
# Carregando a biblioteca da função para gerar o HeatMap - heatmap.2()
library("gplots")
library("RColorBrewer")
# Caminho para o diretório
setwd("/data/bioaat2024/joao/L_salivarius")
# Leitura dos dados
genes_fpkm <- read.delim("./output/hisat_cufflinks_cuffnorm/genes.fpkm_table", sep="\t")
head(genes_fpkm)
colnames_samples <- colnames(genes_fpkm)
samples_only <- colnames_samples[grep("^(LSRG|LSRR)_", colnames_samples)]
gene_diff <- read.delim("./output/hisat_cufflinks_cuffdiff/gene_exp_teste.diff", sep="\t")
head(gene_diff)
description <- read.delim("./description.txt", sep="\t")
head(description)
# Merge dos dados com anotação
df.genes_fpkm.diff <- merge(
x = genes_fpkm,
y = gene_diff,
by.x = 'tracking_id',
by.y = 'test_id'
)
df.genes_fpkm.diff.annot <- merge(
x = df.genes_fpkm.diff,
y = description,
by.x = 'tracking_id',
by.y = 'id'
)
df.diff <- subset(df.genes_fpkm.diff.annot, q_value<0.05 & abs(log2.fold_change.) > 2 )
png(filename = "/tmp/golub.diffs.png",
bg = "white",
res = 300,
width = 3000, height = 3000)
hv1 <- heatmap.2(
as.matrix(log2(df.diff[, samples_only])),
scale = "row",
col = greenred,
labRow = df.diff[, "description"],
key = TRUE,
symkey = TRUE,
density.info = "none",
trace = "none",
Rowv = TRUE,
Colv = TRUE,
cexRow = 0.8,
cexCol = 1,
keysize = 1,
margins = c(10, 25),
dendrogram = c("both"),
main = "HeatMap"
)
graphics.off()
##########################
# Volcano Plot
# Carregando bibliotecas necessárias
library(tidyverse)
library(RColorBrewer)
library(ggrepel)
# Criando uma coluna para expressão diferencial
df.genes_fpkm.diff.annot$diffexpressed <- "Not significant"
df.genes_fpkm.diff.annot$diffexpressed[df.genes_fpkm.diff.annot$log2.fold_change. > 0 & df.genes_fpkm.diff.annot$q_value < 0.05] <- "Upregulated"
df.genes_fpkm.diff.annot$diffexpressed[df.genes_fpkm.diff.annot$log2.fold_change. < 0 & df.genes_fpkm.diff.annot$q_value < 0.05] <- "Downregulated"
# Criando o gráfico
p <- ggplot(data = df.genes_fpkm.diff.annot, aes(x = log2.fold_change., y = -log10(q_value), col = diffexpressed, fill = diffexpressed)) +
geom_vline(xintercept = c(-2, 2), col = "gray", linetype = 'dashed') +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = 'dashed') +
geom_point(size = 2.5, shape = 21, stroke = 0.3, color = "black") +
# position = position_jitter(width = 0, height = 0.1)) + # Ajuda a ver pontos muito próximos
scale_color_manual(
values = c("blue", "grey", "red"),
labels = c("Downregulated", "Not significant", "Upregulated")
) +
scale_fill_manual(
values = c("blue", "grey", "red")
) +
theme_minimal() +
theme(
axis.text = element_text(size = 12, family = "Arial"),
axis.title = element_text(size = 14, family = "Arial"),
legend.title = element_blank(),
text = element_text(family = "Arial")
)
```
</div>
## Comparações com os Resultados obtidos
<div style="text-align: justify;">
Como descrito pelos autores, a capacidade de *L. salivarius* Ren em utilizar diferentes carboidratos foi medida no estudo escolhido, onde, a espécie cresceu de forma limitada em mMRS contendo xilose, celobiose e arabinose. Contudo, em rafinose o crescimento foi, significativamente superior ao valor em glicose. Além disso, a taxa de crescimento de L. salivarius Ren em rafinose foi 0,91 ± 0,03 h−1, superior àquela observada em glicose (0,83 ± 0,01 h−1). E a análise de RNA-seq foi utilizada para entender melhor essa resposta dessa cepa à rafinose.
O número de leituras limpas obtidas de cada amostra e as informações sobre as leituras mapeadas para o genoma de referência também estão apresentadas no texto suplementar do artigo. As tabelas seguintes mostram uma comparação com os resultados obtidos tanto no trabalho, quanto durante a atividade da disciplina
</div>
Tabela 1: Número de leituras antes e após o uso do script de pré-processamento.
| Nome das amostras | Leituras brutas | Leituras limpas no estudo | Leituras limpas duraten a atividade |
|:-----------------:|:---------------:|:-------------------------:|:-----------------------------------:|
| LSRG_1 | 9.321.970 | 9.228.612 | 9.062.236 |
| LSRG_2 | 10.132.006 | 10.031.726 | 9.895.742 |
| LSRG_3 | 8.033.212 | 7.960.520 | 7.848.328 |
| LSRR_1 | 8.192.438 | 8.124.940 | 7.922.904 |
| LSRR_2 | 10.460.828 | 10.379.636 | 10.211.942 |
| LSRR_3 | 8.760.588 | 8.681.970 | 8.529.290 |
<div style="text-align: justify;">
A comparação entre os gráficos gerados pelo MultiQC podem ser verificados logo abaixo, tanto para:
</div>
Antes do processamento;
Quanto para após o processamento.
Tabela 2: Informações das leituras mapeadas para o genoma de referência.
| Nome das amostras | LSRG_1 | LSRG_2 | LSRG_3 | LSRR_1 | LSRR_2 | LSRR_3 |
|:--------------------------:|:----------------:|:----------------:|:----------------:|:----------------:|:-----------------:|:----------------:|
| Total de leituras mapeadas no estudo | 9.118.664, 98.8% | 9.933.683, 99.0% | 7.872.606, 98.9% | 7.980.887, 98.2% | 10.297.765, 99.2% | 8484489, 97.7% |
| Total de leituras mapeadas na atividade | 9.043.232, 99,8% | 9.866.558, 99,7% | 7.818.162, 99,6% | 7.879.784, 99,4% | 10.195.898, 99,8% | 8.400.362, 98,4% |
<div style="text-align: justify;">
A distribuição dos fragmentos formados podem ser verificados logo abaixo para cada uma das amostras, com valor médio maior de 300 pares de bases.
<\div>
<div style="text-align: justify;">
A análise de expressão diferencial de genes diferencialmente expressos (DEGs) foi conduzida, no artigo como menciondo, com valor de P ajustado < 0,05 e considerando para discussão um Log2 Fold Change > 0,8. Para melhor visualizar os genes mais diferencialmente expressos, o valor de LogFC foi ajustado para 2, sendo verificados genes envolvidos no transporte de carboidratos, metabolismo da rafinose, metabolismo da galactose e metabolismo do piruvato, como também verificados tanto no heatmap gerado nessa atividade, quanto na [tabela 1](https://link.springer.com/article/10.1007/s13205-022-03280-6/tables/1) do artigo. Além disso, foi gerado o gráfico Volcano identificar a dispersão de genes diferencialmente expressos, combinando a mudança na expressão (eixo x, log2 do fold change), com o valor negativo do logaritmo base 10 do valor p (eixo y). Está destacado no gráfico também linhas de corte, tanto verticais para fold change utilizado no heatmap, quanto horizontal para significância estatística também utilizada.
heatmap
Volcano
O heatmap ilustra a expressão diferencial de genes nas duas condições diferentes amostradas. As cores indicam o nível de expressão em relação à média, com vermelho representando maior expressão, verde indicando menor expressão e preto simbolizando valores próximos à média. A análise dos dendrogramas revela padrões de agrupamento tanto entre os genes quanto entre as condições, sugerindo que genes com perfis de expressão semelhantes podem estar funcionalmente relacionados ou co-regulados, enquanto condições experimentais próximas podem compartilhar respostas moleculares semelhantes.
Os resultados indicam que alguns genes, como aqueles relacionados ao transporte de açúcares (como subunidades do sistema de transporte PTS) e enzimas metabólicas (como UDP-glucose-4-epimerase Gale), apresentam expressão aumentada em condições específicas, em LSRR_1 e LSRR_2, sugerindo uma possível ativação metabólica sob as maiores concentrações de rafinose. Por outro lado, genes com expressão reduzida em LSRQ_1 e LSRQ_0 podem estar sendo reprimidos ou não necessários frente as sinalizações que não indicam o enriquecimento do meio com rafinose
Os resultados então alcançados pela execução da análise de transcritos com uma abordagem diferente a executada pelos autores corroborá com os achados e apontamentos feitos pelos autores para *L. salivarius* Ren onde seu perfil de transcritoma por RNA-Seq indica um aumentou a expressão de enzimas envolvidas no transporte de rafinose, hidrólise de rafinose, utilização de galactose e metabolismo de sacarose para utilizar rafinose.
<div style="background-color: #ADD8E6; border-left: 5px solid blue; padding: 10px; color: black; font-weight: bold;">
Os scripts de pré-processamento e processamento dos dados podem ser encontrados em suas versões anteriores a atualizada durante a disciplina no GitHub do Prof. Dr.
<a href="https://github.com/bioinfo-fcav/bioapprna" target="_blank">Daniel Guariz Pinheiro</a>
</div>
Sign in with Wallet
Connect another wallet