# Bioinformática II: Análise de Transcriptomas
**<font size="1px" color="green">Atualizado em Nov/2021</font>**
---
URL curta: http://bit.ly/bioinfotranscriptomics
## Sumário
[TOC]
## Simulação
A simulação a seguir, necessita do pacote E-Direct com as ferramentas [E-utilities](https://www.ncbi.nlm.nih.gov/books/NBK25501/). Para maiores informações, consultar a documentação [Entrez Direct E-utilities on the UNIX Command Line](https://www.ncbi.nlm.nih.gov/books/NBK179288/).
### Requisitos
Seguir as instruções para instalação do [e-direct](https://www.ncbi.nlm.nih.gov/books/NBK179288/).
- Caso precise atualizar o conjunto de programas do e-direct, remova primeiro a pasta correspontende
```bash=
cd ~
rm -fr ./edirect
```
- Copiar e colar o código a seguir
```bash=
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
```bash=
echo "export PATH=\$PATH:\$HOME/edirect" >> $HOME/.bash_profile
```
- Para testar:
```bash=
esearch -db nucleotide -query XM_024914431.1 | efetch -format fasta
```
### Construção do transcriptoma para simulação
Teremos que utilizar um transcriptoma com um genoma de referência. O exemplo a seguir utilizará genes de *Trichoderma harzianum* (Agente modelo para uso como biocontrole de fitopatógenos).
::: info
Obtenha o Taxonomy ID para esta espécie: 5544
Dessa forma, para as consultads nos bancos de dados do NCBI, utilizar: **txid5544[Organism:exp]**
:::
A seguir vamos utilizar as ferramentas esearch/efetch para fazer o download e construir o arquivo fasta contendo os transcritos referência para a simulação dos experimentos de investigação do transcriptoma.
Para conhecer os IDs de interesse, consultar a página do banco de dados Gene do NCBI.
e depois de selecionar o gene de interesse:
Identificar o número de acesso do transcrito. No caso de transcritos no banco de dados de referência RefSeq, os IDs iniciam com o prefixo "NM_" ou "NR_", respectivamente para transcritos codificadores de proteínas (mRNAs) e transcritos não codificadores de proteínas (ncRNAS). Os IDs também podem ser respectivamente "XM_" ou "XR_", caso essas entradas tenham somente uma predição computacional, sem ainda terem passado pela etapa de curadoria.
Procurar pelos números de acesso dos transcritos. No caso do gene acima, há 2 transcritos (2 isoformas).
::: success
Números de acesso do gene *M431DRAFT_343923*:
**XM_024914431.1** → XP_024777235.1 hypothetical protein M431DRAFT_343923 [Trichoderma harzianum CBS 226.95]
**XM_024914432.1** → XP_024777236.1 hypothetical protein M431DRAFT_343923 [Trichoderma harzianum CBS 226.95]
:::
```bash=
esearch -db nucleotide -query XM_024914431.1 | efetch \
-format fasta > transcriptoma.fa
esearch -db nucleotide -query XM_024914432.1 | efetch \
-format fasta >> transcriptoma.fa
```
Selecionar mais transcritos de outros genes:
```bash=
esearch -db nucleotide -query XM_024918093.1 | efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_024918802.1 | efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_024914151.1 | efetch \
-format fasta >> transcriptoma.fa
esearch -db nucleotide -query XM_024922757.1 | efetch \
-format fasta >> transcriptoma.fa
```
Com organismos procariotos é necessário copiar e colar manualmente o resultado da busca do NCBI no banco de dados gene com as coordenadas do gene referentes ao transcrito.
Veja exemplo do gene COR42_RS23510 que codifica proteína *DUF3018 family protein* de *Xanthomonas citri* pv. citri strain Xcc29-1 plasmid pXAC64.
```bash=
# COR42_RS23510
efetch -db nucleotide -id NZ_CP024033.1 \
-format fasta \
-chr_start 31467 -chr_stop 31685 -strand minus \
> transcriptoma.fa
```
:::warning
O número de início (*start*) e fim (*stop*) têm que ser ajustados. No caso acima, COR42_RS23510 possui as seguintes coordenadas *31468..31686 complement*, dessa forma, devemos considerar 1 base a menos no início e uma base a menos no fim, ou seja, neste caso 31467 e 31685, respectivamente para início e fim, assim como o fato de estar na fita complementar (*minus*).
:::
E no caso do gene COR42_RS23335 que codifica proteína *recombinase family protein* de *Xanthomonas citri* pv. citri strain Xcc29-1 plasmid pXAC33 ([NCBI Genome](https://www.ncbi.nlm.nih.gov/genome/527?genome_assembly_id=412276)).
```bash=
# COR42_RS23335
efetch -db nucleotide -id NZ_CP024032.1 \
-format fasta \
-chr_start 30895 -chr_stop 31497 -strand plus \
>> transcriptoma.fa
```
:::warning
No caso de COR42_RS23335, o mesmo ajuste deve ser feito.
:::
A fim de simplificar, podemos fazer o mesmo processo utilizando uma estrutura de repetição.
```bash=
rm -f transcriptoma.fa
for acc in XM_024918093.1 XM_024914151.1 XM_024918802.1 XM_024922757.1 XM_024914431.1 XM_024914432.1 ; do
echo "Pegando FASTA para ${acc} ..."
esearch -db nucleotide -query ${acc} | efetch \
-format fasta >> transcriptoma.fa
done
```
Para facilitar ainda mais, podemos criar um arquivo com os números de acesso na primeira coluna e utilizá-lo com o código a seguir.
```bash=
rm -f transcriptoma.fa
for acc in $(cut -f 1 ACCS.txt); do
echo "Pegando FASTA para ${acc} ..."
esearch -db nucleotide -query ${acc} | efetch \
-format fasta >> transcriptoma.fa
done
```
:::warning
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> e não espaços simples**). Isso auxiliará no futuro a identificar quais transcritos relacionam-se com quais genes.
:::
::: info
Arquivo **ACCS.txt**:
<pre>
XM_024918093.1 M431DRAFT_502101
XM_024914151.1 M431DRAFT_300419
XM_024918802.1 M431DRAFT_506422
XM_024922757.1 M431DRAFT_78900
XM_024914431.1 M431DRAFT_343923
XM_024914432.1 M431DRAFT_343923
</pre>
:::
Para procariotos:
::: info
Arquivo **ACCS.txt**:
<pre>
COR42_RS23510 NZ_CP024033.1 31467 31685 minus
COR42_RS23335 NZ_CP024032.1 30895 31497 plus
</pre>
:::
```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
```
::: success
Se executar o comando a seguir e encontrar uma linha para cada transcrito, o processo foi bem sucedido.
```bash=
grep '^>' transcriptoma.fa
```
:::
### Criando os arquivos de abundância
Os arquivos abundance_A.txt e abundance_B.txt devem ser criados manualmente, respectivamente para a condição biológica A e condição biológica B. Os mesmos números de acesso usados na construção do arquivo transcriptoma.fa devem ser utilizados (coluna 1) juntamente com a proporção de cada transcrito (coluna 2) considerando que pode haver diferenças de abundância entre essas duas condições biológicas.
::: danger
**ATENÇÃO:** As duas colunas devem ser **separadas por TAB**, **sem linhas adicionais** e **sem espaços**.
:::
- Para eucariotos
::: info
**abundance_A.txt**
<pre>
XM_024918093.1 0.50
XM_024914151.1 0.05
XM_024918802.1 0.05
XM_024922757.1 0.25
XM_024914431.1 0.15
XM_024914432.1 0
</pre>
**abundance_B.txt**
<pre>
XM_024918093.1 0.05
XM_024914151.1 0.05
XM_024918802.1 0.05
XM_024922757.1 0.25
XM_024914431.1 0.10
XM_024914432.1 0.50
</pre>
:::
- Para procariotos
::: info
**abundance_A.txt**
<pre>
COR42_RS23510 0.8
COR42_RS23335 0.2
</pre>
**abundance_B.txt**
<pre>
COR42_RS23510 0.2
COR42_RS23335 0.8
</pre>
:::
::: success
Execute os comandos abaixo, o resultado da soma dos valores da coluna 2 tem que ser 1 para ambos os arquivos:
```bash=
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
```
:::
#### Genoma reduzido
Caso o genoma seja muito grande, é possível criar uma versão reduzida do genoma somente com os cromossomos dos transcritos selecionados. Para as sequências em fasta:
```bash=
echo -e "NW_020209250\nNW_020209251\nNW_020209249\nNW_020209252" | pullseq -N -i genome.fa > toygenome.fa
```
E para as coordenadas no GFF:
```bash=
grep -P '^(##gff|NW_020209250\t|NW_020209251\t|NW_020209249\t|NW_020209252\t)' genome.gff > toygenome.gff
```
### Gerando os arquivos .fastq
Para exemplificar, abaixo estão apenas os comandos para gerar a primeira réplica da simulação da corrida para a amostra A utilizando o (abundance_A.txt).
1. Gerar "25000" sequências de fragmentos aleatórios em média de "300 pb" com desvio padrão de "30 pb" a partir das sequências do transcriptoma montado anteriormente "transcriptoma.fa", considerando a abundância contida em "abundance_A.txt"
```bash=
generate_fragments.py -r transcriptoma.fa \
-a ./abundance_A.txt \
-o ./tmp.frags.r1 \
-t 25000 \
-i 300 \
-s 30
```
2. Renomear as sequências para "SAMPLEA" juntamente com um número hexadecimal sequencial começando por "0000000001" e posicionar todas as bases em uma única linha (até 1000 pb). A linha de SED irá substituir a descrição adicional.
```bash=
cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \
-if FASTA \
-of FASTA \
-p SAMPLEA1 \
-w 1000 | \
sed 's/^>\(\S\+\).*/>\1/' \
> ./frags_A1.fa
```
3. Fazer a simulação considerando leituras "paired" 2x "151" bp , com parâmetros configurados ([https://www.ebi.ac.uk/goldman-srv/simNGS/runfiles5/151cycPE/s_4_0099.runfile]) usando simNGS:
```bash=
cat ./frags_A1.fa | simNGS -a \
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-p paired \
/usr/local/bioinfo/simNGS/data/s_4_0099.runfile \
-n 151 > ./SAMPLEA1.fastq 2> SAMPLEA1.err.txt
```
4. Desintercalar as leituras em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq no diretório ../raw
```bash=
mkdir -p ./raw
deinterleave_pairs SAMPLEA1.fastq \
-o ./raw/SAMPLEA1_R1.fastq \
./raw/SAMPLEA1_R2.fastq
```
5. Remover arquivos desnecessários:
```bash=
rm -f ./tmp.frags.r1.1.fasta ./frags_A1.fa ./SAMPLEA1.fastq ./SAMPLEA1.err.txt
```
6. Verificar o número de sequências geradas
```bash=
wc -l ./raw/SAMPLEA1_R1.fastq
wc -l ./raw/SAMPLEA1_R2.fastq
```
7. Repetir o mesmo processo para a réplica 2 da condição biológica A (SAMPLEA2) utilizando o arquivo abundance_A.txt
...
8. Repetir o mesmo processo para a réplica 1 e réplica 2 para a condição biológica B (SAMPLEB).
...
### Impacientes
Os impacientes que já tiverem os arquivos **ACCS.txt**, **abundance_A.txt** e **abundance_B.txt** e podem acompanhar a partir deste ponto. Os impacientes que não tiverem os arquivos, manter a calma e ler com atenção as seções anteriores.
:::success
As etapas de simulação realizadas anteriormente podem ser feitas utilizando os scripts [sim.sh](https://github.com/dgpinheiro/bioaat/blob/master/sim.sh) para eucariotos ou [sim-prok.sh](https://github.com/dgpinheiro/bioaat/blob/master/sim-prok.sh) para procariotos.
:::
## Obtenção das referências para análises
### Sequências e Anotações
Sequências no formato FASTA e anotações gênicas no formato GFF.
Página do genoma para *Trichoderma harzianum*. Utilizar a consulta pelo Taxonomy ID ().
- [Sequências genômicas](ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.fna.gz)
- [Anotações gênicas](ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.gff.gz)
```bash=
mkdir ./refs
cd ./refs
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.gff.gz
```
* Utilize gunzip para descompactar
```bash=
gunzip *.gz
```
### Como fazer o download de corridas do SRA:
::: warning
**ATENÇÃO**: NÃO É NECESSÁRIO FAZER O DOWNLOAD DOS DADOS DO SRA CASO ESTEJA PARTINDO DOS DADOS PROVENIENTES DA ETAPA ANTERIOR DE SIMULAÇÃO.
ESTA SEÇÃO É PARA OS QUE DESEJAM APRENDER UMA FORMA DE OBTER OS DADOS DE EXPERIMENTOS REAIS E SEGUIR O MESMO PROTOCOLO. OBVIAMENTE SERÃO NECESSÁRIAS ADAPTAÇÕES NOS NOMES DAS CONDIÇÕES BIOLÓGICAS E DOS ARQUIVOS.
APÓS ESTA SEÇÃO, É POSSÍVEL SEGUIR O TUTORIAL COM OS DADOS DA SIMULAÇÃO.
:::
::: danger
**ATENÇÃO**: NÃO FAZER O DOWNLOAD DOS ARQUIVOS DO NCBI SRA NO SERVIDOR HAMMER!
Os comandos a seguir podem ser executados em um outro computador com mais espaço.
:::
- Criar diretório apropriado
```bash=
cd ../
mkdir ./raw
cd ./raw
```
- Consultar por organismo *Trichoderma harzianum*:
:::info
txid5544[Organism:noexp]
:::
- Identificar o ID da corrida (prefixo SRR) a partir do experimento (prefixo SRX):
- Utilizar o comando *prefetch* do NCBI Toolkit:
```bash=
prefetch SRR8707026
```
- O arquivo .sra será baixado para o diretório ~/ncbi/public/sra/
```bash=
ls -lh ~/ncbi/public/sra/SRR8707026.sra
```
- Obter os arquivos .fastq referentes a esta corrida. O parâmetro --split-3 indica que se deseja separar as sequências em arquivos _1.fastq e _2.fastq. O parâmetro --origfmt indica que se deseja utilizar o nome original das reads.
```bash=
fastq-dump \
--split-3 \
--origfmt SRR8707026
```
:::info
Read 22339605 spots for SRR8707026
Written 22339605 spots for SRR8707026
:::
## Limpeza dos arquivos referência
Voltar para diretório ./refs:
```bash=
cd ../refs/
```
### Limpeza dos cabeçalhos FASTA
Criar arquivo com o script (cleanfasta.sh) abaixo:
```bash=
nano cleanfasta.sh
```
Ctrl+c & Ctrl+v
```bash=
#!/bin/bash
infile=$1
if [ ! ${infile} ]; then
echo "Missing input fasta file"
exit
fi
if [ ! -e ${infile} ]; then
echo "Not found input fasta file (${infile})"
exit
fi
sed 's/^>\([^\.]\+\).*/>\1/' ${infile}
```
Dar permissão:
```bash=
chmod a+x cleanfasta.sh
```
Executar:
```bash=
./cleanfasta.sh GCF_003025095.1_Triha_v1.0_genomic.fna > genome.fa
```
### Ajustes do arquivo GFF:
Executar o script [fixNCBIgff.sh](https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixNCBIgff.sh) do repositório [Github](https://github.com/dgpinheiro/bioinfoutilities):
```bash=
fixNCBIgff.sh GCF_003025095.1_Triha_v1.0_genomic.gff genome.gff
```
::: warning
O script [**fixNCBIgff.sh**](https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixNCBIgff.sh) possui dependências:
- [fixGFFgenestruct.pl](https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixGFFgenestruct.pl)
- [fixGFFdupID.pl](https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixGFFdupID.pl)
- [sortGFF.pl](https://github.com/dgpinheiro/bioinfoutilities/blob/master/sortGFF.pl)
:::
Veja descrição das colunas do arquivo [GFF](https://www.ensembl.org/info/website/upload/gff3.html) e [GFF2/GTF](https://www.ensembl.org/info/website/upload/gff.html).
#### Testes com GFF
- Para a conversão de GFF versão 3 para um arquivo GTF:
```bash=
gffread genome.gff -g genome.fa -T -o genome.gtf
```
- Para obtenção de um arquivo com a sequência FASTA do trasncriptoma:
```bash=
gffread genome.gff -g genome.fa -w transcriptome.fa
```
- Para obtenção de um arquivo com a sequência FASTA do proteoma:
```bash=
gffread genome.gff -g genome.fa -y proteome.fa
```
### Testes de alinhamento com Bowtie2
:::warning
NÃO É NECESSÁRIO REALIZAR ALINHAMENTOS COM BOWTIE2 PARA AS ANÁLISES USANDO O GENOMA REFERÊNCIA. É APENAS **TESTE**.
:::
Baixar o [Integrative Genomics Viewer](https://software.broadinstitute.org/software/igv/) (IGV) para visualização dos alinhamentos.
- Indexação do genoma
```bash=
cd ./refs
bowtie2-build -f genome.fa genome
```
- Alinhamento
```bash=
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
```bash=
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
```
:::info
Devem ser copiados os seguintes arquivos para serem abertos no [IGV](https://software.broadinstitute.org/software/igv/):
**bowtie2-test-sorted.bam**
**bowtie2-test-sorted.bam.bai**
:::
Mapeamento de RNA-Seq reads x Genoma:
:::warning
Observar que há falta de cobertura nas regiões correspondentes a introns.
:::
### Testes de alinhamento com TopHat
::: danger
**ATENÇÃO:** O TopHat tem problemas de compatibilidade com versões mais recentes do bowtie2. Verifique a compatibilidade antes de utilizá-los.
:::
```bash=
tophat2 --num-threads 2 \
--GTF ./refs/genome.gtf \
--microexon-search \
--coverage-search \
--read-edit-dist 2 \
--b2-very-sensitive \
./refs/genome \
./raw/SAMPLEA1_R1.fastq ./raw/SAMPLEA1_R2.fastq
```
No caso do TopHat basta indexar o arquivo .bam resultante para incorporá-lo ao IGV.
```bash=
samtools index ./tophat_out/accepted_hits.bam
```
::: success
Notar o alinhamento das leituras que atravessam os introns e alinham concordantemente com a estrutura gênica contida no arquivo GTF.
:::
### Testes de alinhamento com STAR
- Criação do índice
- Entrar no diretório
```bash=
cd ./refs
```
- Criar diretório
```bash=
mkdir ./star_index
```
- Criação do índice
```bash=
STAR --runThreadN 2 \
--runMode genomeGenerate \
--genomeFastaFiles ./genome.fa \
--genomeDir ./star_index \
--sjdbGTFfile ./genome.gtf \
--sjdbOverhang 149
```
Alinhamento:
```bash=
cd ../
STAR --runThreadN 2 \
--genomeDir ./refs/star_index \
--readFilesIn ./raw/SAMPLEA1_R1.fastq ./raw/SAMPLEA1_R2.fastq \
--sjdbGTFfile ./refs/genome.gtf \
--outFileNamePrefix ./star_out/
```
O resultado no formato .sam deve ser convertido para .bam para que possa ser utilizado no [IGV](https://software.broadinstitute.org/software/igv/).
```
samtools view -b ./star_out/Aligned.out.sam > ./star_out/Aligned.out.bam
samtools sort ./star_out/Aligned.out.bam -o ./star_out/Aligned.out-sorted.bam
samtools index ./star_out/Aligned.out-sorted.bam
```
Procura por: XM_024918093.1
::: success
Notar o alinhamento das leituras que atravessam os introns e alinham concordantemente com a estrutura gênica contida no arquivo GTF.
:::
### Testes com fastqc (pre)
```bash=
mkdir -p output/processed/fastqc/pre
fastqc -t 2 \
./raw/SAMPLEA1_R1.fastq \
-o ./output/processed/fastqc/pre/
fastqc -t 2 \
./raw/SAMPLEA1_R2.fastq \
-o ./output/processed/fastqc/pre/
```
### Testes com prinseq-lite.pl
Poda de qualidade baseada em janela deslizante.
::: spoiler
<pre>
<b>AATAGATACGATGAGAGACCAGTAGACGA</b>TCA
||||||
|||||10
||||20
|||25
||33
|30
30
mean(c(25, 20, 10)) = 18.33
mean(c(33, 25, 20)) = 26
mean(c(30, 33, 25)) = 29.33
mean(c(30, 30, 33)) = 31
-trim_qual_window 3
-trim_qual_step 1
-trim_qual_right 30
-trim_qual_rule lt
-trim_qual_type mean
<b>AATAGATACGATGAGAGACCAGTAGACGA</b>
</pre>
:::
```bash=
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
```
### Testes com fastqc (pos)
```bash=
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/
```
### Criando script com loop
Utilizar nano ou vim para criar o arquivo **./scriptfor.sh** com conteúdo abaixo:
```bash=
#!/bin/bash
mkdir -p ./output/processed/fastqc/pre
for r1 in `ls raw/*_R1.fastq`; do
r2=`echo ${r1} | sed 's/_R1.fastq/_R2.fastq/'`
name=`basename ${r1} | sed 's/_R1.fastq//'`
echo "FastQC evaluation using sample ${name}: ${r1} & ${r2} ..."
fastqc -t 2 \
${r1} \
-o ./output/processed/fastqc/pre/ \
> ./output/processed/fastqc/pre/${name}_R1.log.out.txt \
2> ./output/processed/fastqc/pre/${name}_R1.log.err.txt
fastqc -t 2 \
${r2} \
-o ./output/processed/fastqc/pre/ \
> ./output/processed/fastqc/pre/${name}_R2.log.out.txt \
2> ./output/processed/fastqc/pre/${name}_R2.log.err.txt
done
```
Salvar e sair do editor. Depois, dar permissão de execução e, por fim, executar:
```bash=
chmod a+x ./scriptfor.sh
./scriptfor.sh
```
#### Script FastQC(pre)+PrinSeq+FastQC(pos)
O script [preprocess1.sh](https://github.com/dgpinheiro/bioaat/blob/master/preprocess1.sh) foi desenvolvido para uma análise simples considerando que há arquivos *_R1.fastq e *_R2.fastq no diretório ./raw. O script deve ser gravado no diretório que contém o subdiretório ./raw :
```bash=
wget https://raw.githubusercontent.com/dgpinheiro/bioaat/master/preprocess1.sh
chmod a+x preprocess1.sh
```
... e executado da seguinte forma:
```bash=
./preprocess1.sh
```
O script [preprocess2.sh](https://github.com/dgpinheiro/bioaat/blob/master/preprocess2.sh) tem a mesma funcionalidade que o anterior, porém é muito mais generalista, pois permite a execução a partir de qualquer local, desde que sejam informados o **diretório de entrada** contendo os arquivos *_R1.fastq e *_R2.fastq (1º argumento) e o **diretório de saída** (2º argumento), onde serão gravados os diretórios referentes às análises com FastQC e PrinSeq.
```bash=
wget https://raw.githubusercontent.com/dgpinheiro/bioaat/master/preprocess2.sh
chmod a+x preprocess2.sh
```
```bash=
mkdir -p ./output
./preprocess2.sh ./raw ./output
```
]
### Testes com atropos (*adapter trimming*)
Para chamar o comando atropos no servidor hammer é necessário ativar o ambiente Python compatível com o programa.
```bash=
pyenv activate atropos
```
Com o ambiente carregado é possível invocar o comando e visualizar as opções disponíveis para a poda (trimagem).
```bash=
atropos trim -h
```
Para desativar o ambiente e retornar a o ambiente padrão:
```bash=
pyenv deactivate
```
:::warning
**ATENÇÃO**: Não desativar o comando até terminar a execução do atropos
:::
A primeira execução será no modo **insert**:
```bash=
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
```
Consulte o [manual](https://atropos.readthedocs.io/en/1.1/guide.html) para descrever o uso de cada um desses parâmetros usados na linha de comando. Também é possível consultar por meio da chamada ao help:
```bash=
atropos trim --help
```
A saída do atropos no modo **insert** das *reads* que não foram podadas serão as entradas do atropos no modo **adapter**:
```bash=
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:
```bash=
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
```
Para depois remover os arquivos já concatenados no arquivo com sufixo **.atropos_final.fastq**.
```bash=
rm -f ./output/processed/atropos/SAMPLEA1_R1.atropos_insert.fastq \
./output/processed/atropos/SAMPLEA1_R1.atropos_adapter.fastq \
./output/processed/atropos/SAMPLEA1_R2.atropos_insert.fastq \
./output/processed/atropos/SAMPLEA1_R2.atropos_adapter.fastq
```
Os adaptadores utilizados nesses comandos referem-se à seguinte construção:
E às sequências:
-g <font color="#48E471">AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT</font>
-G <font color="#E46C0A
">GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT</font>
-a <font color="#FF00FF">AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC</font>
-A <font color="#31859C">AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT</font>
## Montando o script de pré-processamento final
A ordem de execução será a seguinte:
**FastQC** (pré) -> **ATROPOS** (insert) -> ATROPOS (adater) -> resultados concatanados ATROPOS (insert) & ATROPOS (adapter) -> **PrinSeq** -> **FastQC** (pós)
Dessa forma, o script [preprocess2.sh](https://github.com/dgpinheiro/bioaat/blob/master/preprocess2.sh) foi atualizado para incorporar as etapas de pré-processamento com ATROPOS no script [preprocess3.sh](https://github.com/dgpinheiro/bioaat/blob/master/preprocess3.sh).
## Montando o script de análise RNA-Seq com referência
O script [rnaseq-ref.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-ref.sh) foi criado utilizando o script [preprocess3.sh](https://github.com/dgpinheiro/bioaat/blob/master/preprocess3.sh) para fazer o pré-processamento e posteriormente executar o STAR com as leituras *Paired-End* e as leituras *Single-End* (*singletons*), além de combinar os resultados, ordenar e executar o script [SAM_nameSorted_to_uniq_count_stats.pl](https://github.com/trinityrnaseq/trinityrnaseq/blob/master/util/misc/SAM_nameSorted_to_uniq_count_stats.pl) para coletar as estatísticas de alinhamento.
```bash=
wget https://raw.githubusercontent.com/dgpinheiro/bioaat/master/preprocess3.sh
wget https://raw.githubusercontent.com/dgpinheiro/bioaat/master/rnaseq-ref.sh
```
Para verificar que o script foi executado com sucesso, confira com a listagem dos arquivos resultantes dos alinhamentos e que foram posteriormente ordenados com samtools.
```bash=
ls -lh ./output/star_out_final/*/Aligned.out.sorted.bam
```
Verifique a porcentagem geral de alinhamento em pares de modo adequado (*proper pairs*):
```bash=
find ./output/star_out_final/ -name 'Aligned.stats.txt' -exec bash -c 'echo $(basename $(dirname $0)) ":" $(grep "^Overall" $0)' {} \;
```
### Execução do pipeline
O script [rnaseq-ref.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-ref.sh) deve ser executado da seguinte forma:
```bash=
./rnaseq-ref.sh ./raw ./output ./refs/genome.gtf ./refs/genome.fa
```
Caso tenha problemas de execução, fazer as seguintes conferências:
#### *CHECKLIST* do sucesso:
:smiley:
- Os seguintes diretórios e arquivos de entrada existem: **./raw**, **./output**, **./refs/genome.gtf** e **./refs/genome.fa**. No caso dos arquivos (genome.gtf e genome.fa), eles devem possuir conteúdo, representado na linha abaixo em número de Bytes, KBytes ou MegaBytes:
```bash=
ls -dlh ./raw ./output ./refs/genome.gtf ./refs/genome.fa
```
- Os arquivos contidos em **./raw** possuem nomenclatura adequada e devem estar em pares R1 e R2:
```bash=
ls -dlh ./raw/*
```
::: success
**NOMEDAAMOSTRA_R[12].fastq**
Exemplo:
SAMPLEB1_R1.fastq e SAMPLEB2_R2.fastq
:::
- Todos os arquivos **./raw** possuem conteúdo, seja por número de Bytes, número de linhas ou quantidade de sequências, respectivamente:
```bash=
ls -dlh ./raw/*
wc -l raw/*.fastq
find ./raw -name '*_R[12].fastq' -exec bash -c 'echo -e "$0\t$(echo "$(cat $0 | wc -l)/4" | bc)"' {} \;
```
... usando um outro código :neutral_face: shell para contar sequências em arquivo .fastq:
```bash=
find ./raw -name '*_R[12].fastq' -exec bash -c 'echo -e "$0\t$(cat $0 | paste - - - - | cut -f 1 | wc -l)"' {} \;
```
### Avaliando a presença de contaminantes
Supondo que possa haver contaminantes na amostra, fazemos um mapeamento das leituras com sequências do contaminante.
No diretório ./refs obter e formatar a(s) sequência(s) do contaminante. Neste caso [*Bacillus subtilis*](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=1423&lvl=3&lin=f&keep=1&srchmode=1&unlock):
```bash=
cd ./refs
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
```
Descompactar:
```bash=
gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz
```
Construir o índice:
```bash=
bowtie2-build GCF_000009045.1_ASM904v1_genomic.fna bsubtillis
```
Alinhar:
```bash
cd ../
bowtie2 -x ./refs/bsubtillis \
-1 ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq_1.fastq \
-2 ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq_2.fastq \
--very-sensitive \
--no-mixed \
--no-discordant \
| samtools view -S -f 2 - | cut -f 1 | nsort -u > ./refs/bsubtillis.txt
pullseq -e \
-n ./refs/bsubtillis.txt \
-i ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq_1.fastq \
> ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq.cleaned_1.fastq
pullseq -e \
-n ./refs/bsubtillis.txt \
-i ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq_2.fastq \
> ./output/processed/prinseq/SAMPLEA1.atropos_final.prinseq.cleaned_2.fastq
```
#### Redução do genoma
Caso estejam limitados com relação à memória e processamento do servidor e tenham que reduzir o tamanho do genoma, devem fazer isso com o arquivo FASTA das sequências genômicas e também o arquivo GTF/GFF.
```bash=
cd ./refs
```
O primeiro passo é criar um arquivo contendo uma lista com os nomes das sequências desejadas.
Vamos selecionar somente os cromossomos referentes aos transcritos selecionados na simulação. Para começar a filtrar as sequências dos cromossomos, vamos montar uma lista não redundante com os números de acesso das sequências de RNA.
```bash=
cut -f 1 abundance_A.txt abundance_B.txt | sort -u > rnas.txt
```
Depois utilizá-la para filtrar as sequências dos cromossomos:
```bash=
grep -w -f ./rnas.txt genome.gff | cut -f 1 | sort -u > selected.txt
```
Para o arquivo com as sequências no formato FASTA:
```bash=
pullseq -n selected.txt -i genome.fa > smallgenome.fa
```
E para o arquivo GFF/GTF:
```bash=
grep -w -f selected.txt genome.gtf > smallgenome.gtf
```
### Montagem dos transcritomas
Para seguir, é necessário ter as sequências das leituras alinhadas com o genoma referência. Neste caso, com o alinhador STAR.
```bash=
tree ./output/star_out_final/
```
:::spoiler
./output/star_out_final/
├── SAMPLEA1
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
├── SAMPLEA2
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
├── SAMPLEB1
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
└── SAMPLEB2
├── Aligned.out.bam
├── Aligned.out.sorted.bam
├── Aligned.stats.txt
├── samtools.merge.log.err.txt
├── samtools.merge.log.out.txt
├── samtools.sort.log.err.txt
└── samtools.sort.log.out.txt
4 directories, 28 files
</pre>
:::
Execução do cufflinks para montagem dos transcritos baseada em alinhamento de leituras (*reads*) com referência genômica.
```bash=
mkdir -p ./output/cufflinks/SAMPLEA1/
cufflinks --output-dir ./output/cufflinks/SAMPLEA1/ \
--num-threads 2 \
--GTF-guide ./refs/genome.gtf \
--frag-bias-correct ./refs/genome.fa \
--multi-read-correct \
--library-type fr-unstranded \
--frag-len-mean 300 \
--frag-len-std-dev 50 \
--total-hits-norm \
--max-frag-multihits 20 \
--min-isoform-fraction 0.20 \
--max-intron-length 10000 \
--min-intron-length 100 \
--overhang-tolerance 4 \
--max-bundle-frags 999999 \
--max-multiread-fraction 0.45 \
--overlap-radius 10 \
--3-overhang-tolerance 300 \
./output/star_out_final/SAMPLEA1/Aligned.out.sorted.bam
```
::: warning
O comando a seguir é um exemplo com a execução para SAMPLEA1.
Para a execução considerando todas as amostras, obter o script
[rnaseq-ref.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-ref.sh) do GitHub.
:::
O script auxiliar [checkMinMaxIntronSize.sh](https://github.com/dgpinheiro/bioinfoutilities/blob/master/checkMinMaxIntronSize.sh) baseado no script [introntab.pl]() o qual foi modificado do script original ([introntab.pl](http://wfleabase.org/genome-summaries/gene-structure/introntab.pl)) obtido de [wFleaBase](http://wfleabase.org/).
```bash=
checkMinMaxIntronSize.sh refs/genome.gff ./output/ 0.9
```
Este número **0.9** que está na linha de comando representa o 90% percentil da distribuição de tamanhos de introns. O resultado deste script são dois valores separados por vírgula, onde o primeiro valor refere-se ao tamanho mínimo e o segundo ao tamanho máximo do intron no arquivo GFF referência. O resultado da linha de comando acima é **127,9697**. Os valores obtidos podem ser aproximados, como foi feito no exemplo acima, 127 para 100 e 9697 para 10000.
::: warning
Veja este [link](https://www.biostars.org/p/344264/) sobre a orientação das sequências em protocolo RNA-Seq fita-específica.
:::
Com a execução do cufflinks, a expectativa é que os arquivos **transcripts.gtf** tenham sido gerados com sucesso.
```bash=
tree ./output/cufflinks/
```
:::spoiler
<pre>
./output/cufflinks/
├── SAMPLEA1
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
├── SAMPLEA2
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
├── SAMPLEB1
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
└── SAMPLEB2
├── genes.fpkm_tracking
├── isoforms.fpkm_tracking
├── skipped.gtf
└── transcripts.gtf
4 directories, 16 files
</pre>
:::
É importante verificar se os arquivos têm conteúdo.
```bash=
ls -lh ./output/cufflinks/*/transcripts.gtf
```
Para visualizar o resultado da quantificação de expressão em FPKM para todos os genes:
```bash=
cat ./output/cufflinks/SAMPLEA1/genes.fpkm_tracking \
| column -t -s$'\t' \
| less -S
```
... ou somente para os genes da simulação:
Com a linha de comando abaixo conseguimos acesso aos IDs dos genes referentes aos RNAs selecionados na simulação.
```bash=
cut -f 1 abundance_A.txt abundance_B.txt \
| sort -u > ./refs/rnas.txt
grep -w -f ./refs/rnas.txt ./refs/genome.gff \
| grep -w -P '\tm?RNA\t' \
| sed 's/^.*Parent=//' \
| sed 's/;.*$//' \
| sort -u > ./refs/genes.txt
```
Vamos utilizar essa lista de genes para filtrar o transcripts.gtf gerado com cufflinks.
```bash=
grep -w -f ./refs/genes.txt ./output/cufflinks/SAMPLEA1/genes.fpkm_tracking \
| column -t -s$'\t' \
| less -S
```
Podemos montar uma outra lista, agora com números de acesso RefSeq dos RNAs selecionados anteriormente. Este ID é próprio do GFF e não corresponde neste caso aos números de acesso RefSeq.
```bash=
grep -w -f ./refs/rnas.txt ./refs/genome.gff \
| grep -w -P '\tm?RNA\t' \
| sed 's/^.*ID=//' \
| sed 's/;.*$//' \
| sort -u > ./refs/rnaIDs.txt
```
... usando essa lista de IDs de RNAs (do GFF/GTF) para obter as coordenadas desejadas.
```bash=
grep -w -f ./refs/rnaIDs.txt ./output/cufflinks/SAMPLEA1/isoforms.fpkm_tracking \
| column -t -s$'\t' \
| less -S
```
### Montagem do transcritoma referência
O transcritoma referência será constituído por uma fusão dos transcritomas individuais gerados com cufflinks. O programa cuffmerge é responsável por esta fusão.
```bash=
find ./output/cufflinks/ -name 'transcripts.gtf' > ./output/cuffmerge/assembly_GTF_list.txt
cuffmerge -o ./output/cuffmerge/ \
--ref-gtf ./refs/genome.gtf \
--ref-sequence ./refs/genome.fa \
--min-isoform-fraction 0.20 \
--num-threads 2 \
./output/cuffmerge/assembly_GTF_list.txt \
> ./output/cuffmerge/cuffmerge.log.out.txt \
2> ./output/cuffmerge/cuffmerge.log.err.txt
```
Para obter uma contagem de ocorrências de montagem de acordo com o [código de classe](http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/index.html) do cufflinks, encontrado na documentação do programa cuffcompare:
```bash=
perl -F"\t" -lane 'my ($transcript_id)=$F[8]=~/transcript_id "([^"]+)"/; my ($class_code)=$F[8]=~/class_code "([^"]+)"/; print join("\t", $transcript_id, $class_code);' \
./output/cuffmerge/merged.gtf \
| nsort -u \
| cut -f 2 \
| nsort \
| uniq -c \
| sed 's/^ *//' \
| awk 'BEGIN{OFS="\t";}{ print $2,$1;} ' \
| sort -t$'\t' -k 2nr,2nr
```
### Quantificação baseada na nova referência (merged.gtf)
Para a quantificação de cada uma das réplicas executamos
o programa cuffquant. Os arquivos são obtidos a partir do
alinhamento realizado anteriormente com STAR (* verificar utilização de ) ou TopHat e :
```bash=
mkdir -p ./output/cuffquant/SAMPLEA1
cuffquant --output-dir ./output/cuffquant/SAMPLEA1 \
--frag-bias-correct ./refs/genome.fa \
--multi-read-correct \
--num-threads 2 \
--library-type fr-unstranded \
--frag-len-mean 300 \
--frag-len-std-dev 50 \
--max-bundle-frags 9999999 \
--max-frag-multihits 20 \
./output/cuffmerge/merged.gtf \
./output/star_out_final/SAMPLEA1/Aligned.out.sorted.bam \
> ./output/cuffquant/SAMPLEA1/cuffquant.log.out.txt \
2> ./output/cuffquant/SAMPLEA1/cuffquant.log.err.txt
```
Ele resultará em arquivos .cxb que poderão ser repassados
ao programa cuffdiff, para o cálculo da expressão gênica diferencial,
ou ao programa cuffnorm, para a obtenção de tabelas de contagens em formato texto.
### Obtenção das matrizes de contagem
O programa cuffnorm será executado para a obtenção das matrizes de contagens de reads por gene/isoforma/cds/tss e normalizadas.
```bash=
mkdir -p ./output/cuffnorm
cuffnorm --output-dir ./output/cuffnorm \
--labels SAMPLEA,SAMPLEB \
--num-threads 2 \
--library-type fr-unstranded \
--library-norm-method geometric \
--output-format simple-table \
./output/cuffmerge/merged.gtf \
./output/cuffquant/SAMPLEA1/abundances.cxb,./output/cuffquant/SAMPLEA2/abundances.cxb \
./output/cuffquant/SAMPLEB1/abundances.cxb,./output/cuffquant/SAMPLEB2/abundances.cxb \
> ./output/cuffnorm/cuffnorm.log.out.txt \
2> ./output/cuffnorm/cuffnorm.log.err.txt
```
Para obter os dados de contagem sem a normalização realizada pelo cuffnorm, é necessário utilizar o script [de-normalize-cuffnorm.R](https://github.com/dgpinheiro/bioinfoutilities/blob/master/de-normalize-cuffnorm.R).
```bash=
de-normalize-cuffnorm.R \
--in=./output/cuffnorm/genes.count_table \
--st=./output/cuffnorm/samples.table \
--out=./output/cuffnorm/genes.raw.count_table
```
### Análises de Expressão Gênica Diferencial
A expressão diferencial será realizada por gene/isoforma/cds/tss e normalizadas a partir do resultado do cuffquant com o programa cuffdiff.
```bash=
mkdir -p ./output/cuffdiff
cuffdiff --output-dir ./output/cuffdiff \
--labels SAMPLEA,SAMPLEB \
--frag-bias-correct ./refs/genome.fa \
--multi-read-correct \
--num-threads 2 \
--library-type fr-unstranded \
--frag-len-mean 300 \
--frag-len-std-dev 50 \
--max-bundle-frags 9999999 \
--max-frag-multihits 20 \
--total-hits-norm \
--min-reps-for-js-test 2 \
--library-norm-method geometric \
--dispersion-method per-condition \
--min-alignment-count 10 \
./output/cuffmerge/merged.gtf \
./output/cuffquant/SAMPLEA1/abundances.cxb,./output/cuffquant/SAMPLEA2/abundances.cxb \
./output/cuffquant/SAMPLEB1/abundances.cxb,./output/cuffquant/SAMPLEB2/abundances.cxb \
> ./output/cuffdiff/cuffdiff.log.out.txt \
2> ./output/cuffdiff/cuffdiff.log.err.txt
```
### Acrescentando informações aos resultados de expressão gênica diferencial
A intenção aqui é mostrar um método para acrescentar informações, ou fundir dados, baseado em coluna comum (usando um **ID**entificador comum) entre os dois arquivos formatados como tabelas, ou matrizes, ou, utilizando uma terminologia de R, **data.frame**s.
Neste exemplo, a seguir, vamos fundir os dados resultantes do cuffdiff (**gene_exp.diff**) com descrições do gene obtidas de uma consulta do NCBI.
Para a consulta ao NCBI, vamos utilizar diretamente via linha de comando,
as ferramentas [E-utilities](https://www.ncbi.nlm.nih.gov/books/NBK25501/) do NCBI.
A seguir, será realizada a busca pelos genes de **Trichoderma harzianum CBS 226.95** utilizando
o Taxonomy ID específico desta linhagem (983964):
```bash=
esearch -db gene -query "983964[Taxonomy ID]" \
| efetch -format tabular > ./refs/gene_result.txt
```
Nem todas as colunas são necessárias, podemos selecionar somente as de interesse. Neste caso, selecionaremos somente as colunas 3, 6 e 8, respectivamente, *GeneID*, *Symbol* e *description*, para um novo arquivo ./refs/**gene_info.txt**.
```bash=
cut -f 3,6,8 ./refs/gene_result.txt > ./refs/gene_info.txt
```
:::warning
O cufflinks/cuffmerge, durante o seu processamento, pode fundir genes considerados inicialmente distintos, ou seja, que possuem IDs distintos (exemplo: M431DRAFT_12071 e M431DRAFT_502167). Isso pode ser observado no arquivo **gene_exp.diff**. Neste caso, o script [mergeR.R](https://github.com/dgpinheiro/bioinfoutilities/blob/master/mergeR.R), criado para realizar a tarefa final de fundir os dados, não irá funcionar, afinal, na linha correspondente, os identificadores estarão unidos por vírgula (exemplo: M431DRAFT_12071,M431DRAFT_502167) impossibilitando o relacionamento. Essa impossibilidade se dá pelo fato de que no arquivo das descrições os dois identificadores estarão cada um em uma linha distinta no arquivo ./refs/**gene_info.txt**.
:::
```bash=
splitteR.R --x="./output/cuffdiff/gene_exp.diff" \
--col.x="gene" \
--by.x="," \
--out="./output/cuffdiff/gene_exp.diff.splitted.txt"
```
Após gerar o novo arquivo, executar a tarefa de fusão dos dados tabulares. Note que utilizaremos o parâmetro **--all.x** para que todas as linhas do arquivo ./output/cuffdiff/**gene_exp.diff.spplitted.txt** estejam no resultado final, no arquivo ./output/cuffdiff/**gene_exp.diff.spplitted.merged.txt**, mesmo que não tenha uma linha correspondente em ./refs/**gene_info.txt**. Essa última situação pode acontecer quando o cufflinks cria um novo mapeamento gênico, que não corresponde a nenhum outro contido no GFF/GTF original, usado como guia para alinhamento e montagem do transcritoma.
```bash=
mergeR.R --x="./output/cuffdiff/gene_exp.diff.spplitted.txt" \
--all.x \
--y="./refs/gene_info.txt" \
--by.x="gene" \
--by.y="Symbol" \
--out="./output/cuffdiff/gene_exp.diff.spplitted.merged.txt"
```
## Montagem *de novo* de transcriptomas
[AULA](https://www.dropbox.com/s/93l12efb021axh1/UNESP-2019-Montagem.pdf?dl=0)
Nesta prática de montagens utilizaremos o programa [Trinity](https://github.com/trinityrnaseq/trinityrnaseq/wiki).
Os dados de entrada para este programa serão os mesmos utilizados com o protocolo baseado em referência usando Cufflinks (Tuxedo). Esses arquivos processados estão no diretório **./output/processed/prinseq**.
Para preparar os dados para montagem com Trinity além de construir a linha de comando, foi desenvolvido o script [rnaseq-denovo.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-denovo.sh) para realizar a montagem *de novo*.
```bash=
./rnaseq-denovo.sh ./output/processed/prinseq ./output
```
::: info
O script [rnaseq-denovo.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-denovo.sh) vai realizar um processo inicial de renomear as leituras para que os cabeçalhos contenham o sufixo **/1** e **/2**, respectivamente para as reads R1 e R2.
Tanto as sequências pareadas quanto as que perderam um dos membros do par (*singletons*) serão renomeadas no script.
:::
::: danger
O script [rnaseq-denovo.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-denovo.sh) assume que as sequências foram previamente processadas com prinseq, e contêm os sufixos **_1.fastq**/**_2.fastq** e **_1_singletons.fastq**/**_2_singletons.fastq**
:::
::: danger
O pipeline de montagem com Trinity funciona nas versões atuais com python 3.6, portanto, é necessário carregar o ambiente com esta versão de python, antes de chamar o pipeline [rnaseq-denovo.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-denovo.sh).
```bash=
pyenv activate trinity
```
:::
O Trinity também tem uma opção para executar a montagem utilizando como guia os alinhamentos das leituras X genoma referência. Neste caso, as entradas são os alinhamentos obtidos com STAR, nas etapas anteriores quando executadas os procedimentos seguindo o protocolo baseado em referência (Tuxedo). Para executar este modo de montagem, foi criado o script [rnaseq-ref-trinity.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-ref-trinity.sh). Neste script, o diretório que serve como entrada é o **./output/star_out_final**, onde estão os arquivos com nome **Aligned.out.sorted.bam**. A primeira etapa deste script antes de realizar a montagem é fundir todos esses alinhamentos em um único arquivo (**./output/All.sorted.bam**).
```bash=
./rnaseq-ref-trinity.sh ./output/star_out_final ./output
```
## Para avaliação das montagens
Indexação do transcriptoma utilizado para a simulação:
```bash=
makeblastdb -title "RefTrans" \
-dbtype nucl \
-out RefTrans \
-in transcriptoma.fa \
> makebleastdb.log.out.txt \
2> makeblastdb.log.err.txt
```
:::warning
Para conferir se esta indexação foi bem sucedida, neste caso, aparecerão os seguintes arquivos:
- RefTrans.nhr
- RefTrans.nin
- RefTrans.nsq
:::
Para executar a comparação vamos utilizar o **blastn**:
Primeiro, para a montagem com Trinity *de novo*:
```bash=
blastn -max_hsps 1 \
-max_target_seqs 1 \
-num_threads 8 \
-query ./output/trinity_assembled/Trinity.fasta \
-task blastn \
-db ./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
```
Visualizar os nomes das sequências *query* (qseqid) e *subject* (sseqid) HSPs de todos HITs obtidos:
```bash=
cut -f 1,2 ./Trinity_x_RefTrans.txt
```
Contar quantas ocorrências de sequências da base de dados (*subject*) foram encontradas:
```bash=
cut -f 2 ./Trinity_x_RefTrans.txt | sort | uniq -c
```
Depois, para a montagem com Trinity *Genome Guided*:
```bash=
blastn -max_hsps 1 \
-max_target_seqs 1 \
-num_threads 8 \
-query ./output/trinity_GG_assembled/Trinity-GG.fasta \
-task blastn \
-db ./RefTrans \
-out ./Trinity-GG_x_RefTrans.txt \
-evalue 1e-5 \
-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
> ./Trinity-GG_x_RefTrans.log.out.txt \
2> ./Trinity-GG_x_RefTrans.log.err.txt
```
Visualizar os nomes das sequências *query* (qseqid) e *subject* (sseqid) HSPs de todos HITs obtidos:
```bash=
cut -f 1,2 ./Trinity-GG_x_RefTrans.txt
```
Contar quantas ocorrências de sequências da base de dados (*subject*) foram encontradas:
```bash=
cut -f 2 ./Trinity-GG_x_RefTrans.txt | sort | uniq -c
```
::: warning
**ATENÇÃO:** Se o número de transcritos montados for menor do que o esperado pela simulação, reveja os parâmetros, por exemplo, o tamanho mínimo de contigs gerados por padrão no script é de 600 (--min_contig_length 600), e se o tamanho do transcrito esperado for menor que este tamanho, ele será omitido do resultado final.
:::
## Para estimar a abundância e encontrar os genes/isoformas diferencialmente expressos
A seguir, vamos obter os valores de *Log Fold-Change* e os p-valores ajustados referentes às comparações pareadas entre **SAMPLEA** e **SAMPLEB**.
Para isso utilizaremos scripts auxiliares do programa Trinity.
Para começar a análise, a variável ambiente do linux ${TRINITY_HOME} deve apontar para o caminho do diretório base do programa Trinity.
Para checar:
```bash=
echo ${TRINITY_HOME}
```
...e para testar se os scripts auxiliares align_and_estimate_abundance.pl e abundance_estimates_to_matrix.pl funcionam e verificarmos as opções de parâmetros:
Primeiro este:
```bash=
${TRINITY_HOME}/util/align_and_estimate_abundance.pl
```
... depois este:
```bash=
${TRINITY_HOME}/util/abundance_estimates_to_matrix.pl
```
Além disso utilizaremos o script [run-DESeq2.R](https://github.com/dgpinheiro/bioinfoutilities/blob/master/run-DESeq2.R) para executar uma análise comparativa pareada de expressão gênica diferencial.
Verificar se ele está disponível para execução:
```bash=
run-DESeq2.R
```
Para executar esse protocolo de análise de expressão gênica diferencial de uma forma automatizada, foi criado o script [rnaseq-trinity-abundance.sh](https://github.com/dgpinheiro/bioaat/blob/master/rnaseq-trinity-abundance.sh).
```bash=
rnaseq-trinity-abundance.sh ./output/renamed/ \
./output/trinity_assembled/ \
./output
```
Para visualizar o resultado da anális com DESeq2 de modo alinhado:
```bash=
cat ./output/abundance/DEG/SAMPLEA-SAMPLEB.txt | column -s$'\t' -t | less -S
```
### Acessando o RStudio para algumas análises dos resultados
Para acessar a interface do RStudio no servidor utilizar este [link](
http://hammer.fcav.unesp.br:8787/).
Na interface do RStudio inserir o script para selecionar os genes diferencialmente expressos:
```R=
# Carregando a biblioteca xlsx para lidar com arquivos do Excel
library("xlsx")
# Carregando a biblioteca da função para gerar o HeatMap - heatmap.2()
library("gplots")
# Carregando a biblioteca com a função para criar a paleta de cores
library("RColorBrewer")
# Equivalente ao comando "pwd" no linux
getwd()
# Equivalente ao comando "cd" no linux
setwd("~/")
getwd()
# Carrega arquivo texto separado por TAB ("\t") com cabeçalho e sem transformar
# strings em fatores
deg.df <- read.delim(file="./output/abundance/DEG/SAMPLEA-SAMPLEB.txt",
sep="\t",
header=TRUE,
stringsAsFactors = FALSE)
# exibir somente as 2 primeiras linhas para checagem
head(deg.df,2)
# Dimensão do data.frame, onde o primeiro valor refere-se à quantidade de linhas e
# o segundo valor à quantidade de colunas
dim(deg.df)
# selecionar um subconjunto de linhas do data.frame deg.df em um outro data.frame
# ("subset.deg.df") contendo somente os genes que possuírem logFC >= 2 ou logFC <= -2,
# além de possuírem um p-valor corrigido ("FDR") <= 0.05
subset.deg.df <- subset(deg.df, (((logFC >= 2) | (logFC <= -2) ) & (FDR <=0.05)) )
dim(subset.deg.df)
head(subset.deg.df,2)
# criando um vetor de nomes de colunas selecionadas,
# ou seja, colunas que não são "X", "logFC", "PValue" ou "FDR",
# as quais, sabemos previamente que contém os valores de expressão.
# Os nomes das colunas serão ordenados depois de secionados com a
# função "setdiff", a qual obtém
# um conjunto de valores (vetor) coma a diferença entre dois conjuntos de
# valores: o de todas as colunas ("colnames(subset.deg.df)") menos os nomes
# das colunas indesejadas.
sel_columns <- sort(
setdiff( colnames(subset.deg.df),
c("X", "logFC", "PValue", "FDR")
)
)
# criando uma nova matriz "expression_data" que contém somente as colunas
# selecionadas
expression_data <- as.matrix( subset.deg.df[,sel_columns] )
# Atribuindo para os nomes das linhas da matriz o conteúdo da coluna "X",
# onde sabemos previamente que contém o identificador dos genes
rownames(expression_data) <- subset.deg.df$X
head(expression_data,2)
# Função para gravar o data.frame em um arquivo Excel, na planilha
# de nome "DEGS_SAMPLEA-SAMPLEB"
write.xlsx(subset.deg.df,
file="./output/abundance/DEG/SAMPLEA-SAMPLEB.xlsx",
sheetName="DEGS_SAMPLEA-SAMPLEB"
)
# cria uma paleta personalizada de 299 cores do vermelho ao verde,
# passando pelo amarelo
my_palette <- colorRampPalette(c("red", "yellow", "green"))(n = 299)
expression_data_to_plot <- log2(expression_data+1)
retval <- expression_data_to_plot
retval$rowMeans <- rm <- rowMeans(expression_data_to_plot, na.rm = TRUE)
expression_data_to_plot <- sweep(expression_data_to_plot, 1, rm)
retval$rowSDs <- sx <- apply(expression_data_to_plot, 1, sd, na.rm = TRUE)
scaled_expression_data_to_plot <- sweep(expression_data_to_plot, 1, sx, "/")
# define as quebras das cores manualmente para transição de cores
col_breaks = c(seq(-1,-0.5,length=100), # for red
seq(-0.49,0.5,length=100), # for yellow
seq(0.51,1,length=100)) # for green
# criação de uma imagem de tamanho 5 x 5 polegadas
png("./output/abundance/DEG/heatmap_DEGS_SAMPLEA-SAMPLEB.png", # cria arquivo do tipo PNG for the heat map
width = 5*300, # 5 x 300 pixels de largura
height = 5*300, # 5 x 300 pixels de altura
res = 300, # 300 pixels por polegada
pointsize = 8) # tamanho da fonte
heatmap.2(scaled_expression_data_to_plot, # a matriz com os valores de expressão
main = "Correlation distance (z-score)", # Título do HeatMap
density.info="none", # desabilita o gráfico de densidade dentro da legenda
trace="none", # desabilita as linhas dentro do HeatMap
margins =c(12,12), # definiação das margens no entorno do gráfico
col=my_palette, # nome do objeto contendo a paleta de cores criada anteriormente
breaks=col_breaks, # pontos de quebra para a transição de cores
symbreaks=TRUE,
dendrogram="both", # desenhar dendrograma para linhas e colunas
distfun = function(x) as.dist(1-cor(t(x))), # distância baseada em correlação
hclustfun = function(x) hclust(x, method="centroid") # método de ligação pelo centróide
)
dev.off() # fecha o arquivo da imagem PNG
## OR ##
# criação de uma imagem de tamanho 5 x 5 polegadas
png("./output/abundance/DEG/heatmap_DEGS_SAMPLEA-SAMPLEB.png", # cria arquivo do tipo PNG for the heat map
width = 5*300, # 5 x 300 pixels de largura
height = 5*300, # 5 x 300 pixels de altura
res = 300, # 300 pixels por polegada
pointsize = 8) # tamanho da fonte
heatmap.2(expression_data_to_plot, # a matriz com os valores de expressão
main = "Correlation distance", # Título do HeatMap
density.info="none", # desabilita o gráfico de densidade dentro da legenda
trace="none", # desabilita as linhas dentro do HeatMap
margins =c(12,12), # definiação das margens no entorno do gráfico
col=my_palette, # nome do objeto contendo a paleta de cores criada anteriormente
scale="row", # pontos de quebra para a transição de cores
symbreaks=TRUE,
dendrogram="both", # desenhar dendrograma para linhas e colunas
distfun = function(x) as.dist(1-cor(t(x))), # distância baseada em correlação
hclustfun = function(x) hclust(x, method="centroid") # método de ligação pelo centróide
)
dev.off() # fecha o arquivo da imagem PNG
```
::: danger
Troque o caminho da pasta **~/** pelo nome da sua pasta contendo as análises.
:::
Há duas versões utilizando os dados padronizados por linha (Z-score). Na primeira versão, acontece primeiro a padronização manual dos valores da estimativa de abundância em log2. Em seguida, há a definição dos valores onde haverá transições de cores.
A segunda versão é mais automatizada, e, portanto, não é preciso especificar as transições. Dessa forma, ela é mais rápida de escrever, porém não permite uma especificação mais precisa das quebras para a transição de cores.
Sign in with Wallet
Connect another wallet