# Relatório final da disciplina Bioinformática Aplicada II: Análise de Transcritomas
Universidade Estadual Paulista
Faculdade de Ciências Agrárias e Veterinárias (FCAV)
Disciplina: Bioinformártica Aplicada II - Análise de transcritomas
Discente: Bruna Marques Moreno
Mestranda em Genética e Melhoramento de Plantas
Contato: bruna.marques.moreno@usp.br
O organismo modelo para análise escolhido foi o feijoeiro (*Phaseolus vulgaris*). Seu taxonomy ID é 3885 na plataforma do NCBI.
Para facilitar o processamento durante a aula, foram selecionados 6 genes de um mesmo cromossomo, além da redução do genoma de referência e arquivo de extensão .gff.
# 1. Obtenção dos dados para análise
Para iniciar as análises foram contruídos os diretórios, em que cada um recebeu um tipo de extensão de dados e resultados de análises
```
mkdir /data/bioaat/bmoreno/pvulgaris
mkdir /data/bioaat/bmoreno/pvulgaris/input
mkdir /data/bioaat/bmoreno/pvulgaris/output
mkdir /data/bioaat/bmoreno/pvulgaris/ref
mkdir /data/bioaat/bmoreno/pvulgaris/raw
```
* /pvulgaris = é o diretório central da análise deste mini projeto de análise de RNA-seq;
* /pvulgaris/input = contém os dados de sequenciamento que serão analisados nomeados de maneira apropriada;
* /pvulgaris/output = contém todos os dados pós processamento, ou seja, resultado da análise com o programa FastQC, o s arquivos resultantes da trimagem, alinhamento, montagem e por fim a tabela que indica os genes diferentemente expressos;
* /pvulgaris/ref = contém os arquivos de referência como o genoma de referência reduzido (com apenas o cromossomo que tem os 6 genes analisados) e o arquivo gff (também reduzido);
* /pvulgaris/raw = contém os dados brutos de sequenciamento (que no exemplo da aula foram gerados com simulação e com os arquivos de abundância de cada gene para cada condição (A e B)).
Esta estrutura está representada no esquema abaixo:
Os diferentes tons de verde representam os diretórios dentro de outros diretórios. O mais escuro (bmoreno) engloba os outros mais claros.
## Preenchendo a pasta **ref**
### Baixando o genoma de referência e arquivos GFF e GTF
Dentro do diretório ref, foi baixado o genoma de referência do *P. vulgaris* e seu arquivo gtf do site NCBI, utilizando o protocolo ftp.(https://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/Phaseolus_vulgaris/latest_assembly_versions/GCF_000499845.1_PhaVulg1_0/)
```
cd ref
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/Phaseolus_vulgaris/latest_assembly_versions/GCF_000499845.1_PhaVulg1_0/GCF_000499845.1_PhaVulg1_0_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/Phaseolus_vulgaris/latest_assembly_versions/GCF_000499845.1_PhaVulg1_0/GCF_000499845.1_PhaVulg1_0_genomic.gtf.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/Phaseolus_vulgaris/latest_assembly_versions/GCF_000499845.1_PhaVulg1_0/GCF_000499845.1_PhaVulg1_0_genomic.gff.gz
```
Para descompactar eles, usou-se o comando "gunzip".
```
gunzip GCF_000499845.1_PhaVulg1_0_genomic.fna.gz
gunzip GCF_000499845.1_PhaVulg1_0_genomic.gtf.gz
gunzip GCF_000499845.1_PhaVulg1_0_genomic.gff.gz
```
Para limpar o arquivo fasta de genoma
```
nano cleanfasta.sh #Criação de um script para isso
```
```
#!/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}
```
```
chmod a+x cleanfasta.sh #Dar permissão
```
```
./cleanfasta.sh GCF_000499845.1_PhaVulg1_0_genomic.fna > genoma.fa #Para limpar e redirecionar os resultados para o arquivo genome.fa
```
Para selecionar apenas o cromossomo 10 do arquivo de genoma já limpo.
```
echo -e "NC_023750" | pullseq -N -i genoma.fa > genome.fa
```
Para arrumar o arquivo gff, utilizou-se o script [fixNCBIgff.sh](https://github.com/dgpinheiro/bioinfoutilities/blob/master/fixNCBIgff.sh) (que está no Github)
```
fixNCBIgff.sh GCF_000499845.1_PhaVulg1_0_genomic.gff genome.gff
```
Para trabalhar apenas com a porção do gff referente ao cromossomo 10:
```
grep -P '^(##gff|NC_023750.1\t)' genoma.gff > genome.gff
```
Para criar um arquivo gtf a partir de um gff
```
gffread genome.gff -g genome.fa -T -o genome.gtf
```
Para criar um transcriptoma de referência baseado em arquivo gff
```
gffread genome.gff -g genome.fa -w transcriptome.fa
```
Para criar um proteoma de referência baseado em arquivo gff
```
gffread genome.gff -g genome.fa -y proteome.fa
```
## Preenchendo a pasta **raw**
Os 6 genes selecionados do cromossomo 10 (RefSeq NC_023750.1) a partir do próprio site do NCBI:
*gene (número de acesso do NCBI da sequência)*
* PHAVU_010G142600g (XM_007135533.1)
* PHAVU_010G142900g (XM_007135538.1)
* PHAVU_010G142900g (XP_007135599.1)
* PHAVU_010G143200g (XM_007135547.1)
* PHAVU_010G143600g (XM_007135551.1)
* PHAVU_010G143900g (XM_007135555.1)
> Um detalhe que passou despercebido foi que um dos transcritos escolhidos na verdade não era um mRNA, então ele acabou não sendo analisado posteriormente.
Depois de selecionados criou-se um arquivo txt.
```
cd raw
nano ACCS.txt
```
E, dentro deste txt, as seguintes informações separadas por **tabs**.
> XM_007135533.1 PHAVU_010G142600g
> XM_007135538.1 PHAVU_010G142900g
> XP_007135599.1 PHAVU_010G142900g
> XM_007135547.1 PHAVU_010G143200g
> XM_007135551.1 PHAVU_010G143600g
> XM_007135555.1 PHAVU_010G143900g
Para coletar as sequências fasta de cada gene, utilizou-se o script "pegargene.sh"
```
nano pegargene.sh
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
```
```
chmod a+x pegargene.sh
```
```
./pegargene.sh
```
Depois para criar abundâncias para cada gene simulando um sequenciamento da vida real para cada condição criou-se dois arquivos txt: abundance_A e abundance_B. **(Sempre separado po TABS)**
```
nano abundance_A.txt
XM_007135533.1 0.02
XM_007135538.1 0.08
XM_007135547.1 0.50
XM_007135551.1 0.3
XM_007135555.1 0.1
nano abundance_B.txt
XM_007135533.1 0.4
XM_007135538.1 0.02
XM_007135547.1 0.5
XM_007135551.1 0.07
XM_007135555.1 0.01
```
Para simular dados reais, foi utilizado o script simLib.pl que já estava no PATH. Foram simuladas 1000 sequências para cada tratamento em três repetições, paired-end com tamanho de 151 nt.
```
simLib.pl -a abundance_A.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA1
simLib.pl -a abundance_A.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA2
simLib.pl -a abundance_A.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEA3
simLib.pl -a abundance_B.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB1
simLib.pl -a abundance_B.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB2
simLib.pl -a abundance_B.txt -i transcriptoma.fa -rf /usr/local/bioinfo/simNGS/data/s_4_0099.runfile -n 1000 -p SAMPLEB3
```
## Preenchendo a pasta **input**
No nosso exemplo de aula, para simular os dados que serão utilizados nas análises para a montagem do transcritoma e análise dos genes diferentemente expressos, os dados do diretório raw foram copiados para o diretório input, com alteração dos nomes para TRAT (tratamento) e CNTL (controle).
```
cd ../raw
#Para compactar os dados da mesma forma que vamos receber da facility
gzip *.fastq
#Arquivo utilizado para alterar os nomes dos arquivos de forma automatizada
nano data.txt
SAMPLEA1_R1.fastq.gz CNTRL_B1_R1.fastq.gz
SAMPLEA1_R2.fastq.gz CNTRL_B1_R2.fastq.gz
SAMPLEA2_R1.fastq.gz CNTRL_B2_R1.fastq.gz
SAMPLEA2_R2.fastq.gz CNTRL_B2_R2.fastq.gz
SAMPLEA3_R1.fastq.gz CNTRL_B3_R1.fastq.gz
SAMPLEA3_R2.fastq.gz CNTRL_B3_R2.fastq.gz
SAMPLEB1_R1.fastq.gz TREAT_B1_R1.fastq.gz
SAMPLEB1_R2.fastq.gz TREAT_B1_R2.fastq.gz
SAMPLEB2_R1.fastq.gz TREAT_B2_R1.fastq.gz
SAMPLEB2_R2.fastq.gz TREAT_B2_R2.fastq.gz
SAMPLEB3_R1.fastq.gz TREAT_B3_R1.fastq.gz
SAMPLEB3_R2.fastq.gz TREAT_B3_R2.fastq.gz
```
Para cria o link simbólico com os nomes da maneira que pensamos melhor em trabalhar
```
createSymLink.sh ./input/data.txt ./raw/ ./input/
```
# 2.Pré-processamento dos dados ##
Durante as aulas foi utilizado um script preprocess5.sh disponível em /data/bioaat/scripts. Nesta etapa foram utilizados três programas: 1.FastQC - para a avaliação da qualidade das sequências, antes e depois do tratamento, 2. Atropos - para a trimagem de insertos provenientes da biblioteca da illumina no processo de sequenciamento, 3.Prinseq - para a trimagem e controle de qualidade de reads (mais comum para análises de metagenoma), que também consegue retirar contaminantes.
O script tem a porção inicial de criação dos diretórios que serão usados como local do output das anáises feitas durante o processamento do script.
## 2.1.Trimagem
(para o atropos funcionar de maneira adequada, o comando conda activate python2 tem que estar ativo).
Primeiramente o script vai analisar os arquivos .fastq utilizando o FastQC a fim de verificar como estão os dados e se é necessário processo de trimagem.
O resultado em formato html está representado abaixo de um dos tratamentos (SAMPLEA1) paras as reads foward e reverse, respectivamente.
R1 (foward) - pré trimagem
R2 (reverse) - pré trimagem
```
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
```
> Acima é um trecho do script que avalia a porção pre do FastQC para as reads R2. O argumento THREADS indica a quantidade de threads destinamos para a análise na linha de comando. o r2, verifica se existe um arquivo R2 nos dados. Em -o é o caminho do output do programa.
Posteriormente, o uso do atropos é pra a trimagem de sequências de baixas qualidade e de bases de baixa qualidade (abaixo de phred 20 - faixa amarela no FastQC).
```
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
```
```
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
```
> Ponto importante de mencionar é que o atropos gera três arquivos *insert*, *adapter* e *untrimmed*. No final o *insert* e *adapter* são concatenados em um chamado *final*.
Um ponto importante é que nas reads R1 tem a presença da cauda poli-A, por conta do enriquecimento feito para a captura de mRNA na fase de montagem da biblioteca para o sequenciamento.
```
prinseq-lite.pl -fastq ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq \
-fastq2 ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq \
-out_format 3 \
-trim_qual_window 3 \
-trim_qual_step 1 \
-trim_qual_type mean \
-trim_qual_rule lt \
-trim_qual_right 30 \
-out_good ${outdir}/processed/prinseq/${name}.atropos_final.prinseq \
-out_bad null \
-lc_method dust \
-lc_threshold 50 \
-min_len 20 \
-trim_tail_right 5 \
-trim_tail_left 5 \
-ns_max_p 80 \
-min_qual_mean 28 \
> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.out.log \
2> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.err.log
```
Após a trimagem:
> É perceptível a diferença da qualidade dos dados, uma vez que a porção final foi trimada por ter valores de phred score abaixo do valor indicado como bom (28, faixa verde da imagem). Apesar de ter melhor qualidade, sequências foram perdidas (1039 antes - 927 depois).
# 3. Alinhamento de reads com genoma de referência
Para as próximas etapas, foi utilizado o script rnaseq-ref.sh (direcionado para eucariotos) em que contém códigos para dois alinhadores: 1.STAR, 2.TopHat2. A decisão de qual usar é feita na linha de comando. Alpem disso têm dois montadores: 1. Cufflinks, 2. Stringtie.
No início do script há a definição das variáveis na linha de comando, como qual a posição do input, output, # de threads a serem usadas, criação de diretórios conforme o alinhador escolhido.
As informações a serem colocadas na linha de comando:
> rnaseq-ref.sh alinhador diretorio-de-input diretório-de-output threads arquivo-gtf genoma-de-referência montador arquivo-de-contaminantes
O primeiro comando utilizado/analisado foi:
```
rnaseq-ref.sh star ./input ./output 3 ./ref/genome.gtf ./ref/genome.fa stringtie
```
Para o STAR funcionar adequadamente, ele precisa de um index. No script é feito um index e alocado no diretório **../output/star_index**
Feito isso, o script passa para esta seguinte parte:
```
STAR --runThreadN ${THREADS} \
--genomeDir ${outdir}/star_index/ \
--readFilesIn ${r1} ${r2} \
--outSAMstrandField intronMotif \
--outFilterIntronMotifs RemoveNoncanonical \
--sjdbGTFfile ${refgtf} \
--outFilterMultimapNmax 20 \
--outFileNamePrefix ${outdir}/star_out_pe/${name}/ \
--outSAMtype BAM Unsorted \
--outFilterType BySJout \
--outSJfilterReads Unique \
--alignSJoverhangMin 8 \
--alignSJDBoverhangMin 1 \
--outFilterMismatchNmax 999 \
--outFilterMismatchNoverReadLmax 0.04 \
--alignIntronMin 85 \ #alteração feita aqui
--alignIntronMax 1200 \ #alteração feita aqui
--alignMatesGapMax 200 \
> ${outdir}/star_out_pe/${name}.log.out.txt \
2> ${outdir}/star_out_pe/${name}.log.err.txt
```
As alterações feitas nesta porção são resultantes do script checkMinMaxIntronSize.sh , em que a linha de comando usada com ele foi:
```
checkMinMaxIntronSize.sh refs/genome.gff ./output/ 0.9
```
Ele procura um valor mínimo e máximo de íntron conforme o arquivo gff que otimiza o processo do alinhamento.
Depois do processamento, se criaram três diretórios: star_out_pe (com alinhamento das reads que permaneceram paired), star_out_se (alinhamento das reads singletons), star_out_final (união dos alinhamentos anteriores).
Dentro da star_out_final:
Dentro deles há o arquivo **.bam** que armazena os dados do alinhamento (sua vizualização é problemática porque é um arquivo binário (binary alinment map). Se transformar em **.sam**, conseguimos ler.
# 4. Montagem do transcriptoma
Para a montagem, o StringTie. Apenas o valor de -g foi alterado para o default.
```
stringtie --merge \
-G ${refgtf} \
-o ${outdir}/${aligner}_stringmerge/merged.gtf \
-m 200 \
-c 1 \
-F 4 \
-T 4 \
-f 0.25 \
-g 50 \ #alterado para o valor default
${outdir}/assembly_GTF_list.txt \
> ${outdir}/${aligner}_stringmerge/stringmerge.out.txt \
2> ${outdir}/${aligner}_stringmerge/stringmerge.err.txt
```
Dentro do diretório tem um diretório para cada condição e repetição, como pode ser visto abaixo.
Além disso, dentro de cada diretório, existem 10 arquivos, dentre eles arquivos de abundância, cobertura. Abaixo é um exemplo dentro do diretório da condição TRAT, da repetição 3.
Após o alinhamento de cada condição é feito uma unição (*merge*) de todos, como um transcriptoma único de referência, que está localizado no diretório ./output/star_stringmerge
Dentro do merged.gtf
```
more merged.gtf
```
# 5. Comparação com o transcriptoma de referência
A próxima etapa foi avaliar a expressão diferencial dos genes. Para tanto, a ferramenta utilizada foi o Cufflinks, mais especificamente o cuffcompare, cuffquant, cuffnorm e cuffdiff.
Com o cuffcompare, inicia-se a comparação entre as condições, no caso da aula CNTLxTRAT.
```
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
```
cuffcmp.combined.gtf:
Com o cuffquant, ele verifica com base na quantidade de reads alinhadas um valor para comparação.
```
cuffquant --output-dir ${outdir}/${aligner}_${assembler}_cuffquant/${name} \
--frag-bias-correct ${refseq} \
--multi-read-correct \
--num-threads ${THREADS} \
--library-type fr-firststrand \ #sequenciamento fita específica
--frag-len-mean 400 \ #mudança pelo nosso exemplo de 2x151 pb
--frag-len-std-dev 50 \
--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
```
Ao final, é gerado um arquivo de extensão .cxb para cada condição e repetição com sua respectiva abundância. No exemplo abaixo o diretório da condição Controle, repetição 1.
# 6. Análise de expressão diferencial
Como os genes apresentam expressão diferente, a quantidade de reads que temos de cada um é diferente. Para poder comparar os genes com quantidades de reads diferentes é necessário normalizar os dados. A etapa do cuffnorm realiza isso.
```
cuffnorm --output-dir ${outdir}/${aligner}_${assembler}_cuffnorm \
--labels $(IFS=, ; echo "${biogroup_label[*]}") \
--num-threads ${THREADS} \
--library-type fr-firststrand \ #mudança por ser seq. fita específica
--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
```
Por exemplo, o arquivo genes.fpkm_table apresenta os dados de cada gene normalizado:
Além disso, há o arquivo com a contagem de reads sem a normalização. Isso porque os programas em R (DESeq2, edgeR) necessitam de dados não normalizados, porém o outpur do Cufflink é normalizado. Assim, uma denormalização é feita para a análise no R, representada abaixo.
Para quantificar a expressão diferencial:
```
cuffdiff --output-dir ${outdir}/${aligner}_${assembler}_cuffdiff \
--labels $(IFS=, ; echo "${biogroup_label[*]}") \
--frag-bias-correct ${refseq} \
--multi-read-correct \
--num-threads ${THREADS} \
--library-type fr-firststrand \ #seq fita específica
--frag-len-mean 400 \ #tamanho médio considerando o sequeciamento 2x151 pb
--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 \
${transcriptomeref} \
${biogroup_files[*]} \
> ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.out.txt \
2> ${outdir}/${aligner}_${assembler}_cuffdiff/cuffdiff.log.err.txt
```
> Para ficar melhor a visualização:
```
grep 'yes' gene_exp.diff | column -t -s$'\t'
```
Além dos genes diferentemente expressos, também tem a informação das isoformas, tss (*transcription starting sites*) e cds.
# 7. Heatmap
Tendo em mãos a expressão dos genes, realizou-se um gráfico do tipo *heatmap* para ilustrar a expressão de todos os genes de uma só vez.
O gráfico foi feito usando o R no servidor. Para acessá-lo em um browser:
```
http://200.145.102.33:8787/
```
O login é o mesmo do servidor.
Nesta etapa usamos a planilha com os dados de FPKM, log2FoldChange. Para acessar a tabela com os dados de FPKM das diferentes condições e suas repetições:
```
#dados de FPKM, tirados diretamente do servidor
pvulgaris.df <- read.delim(file = "/data/bioaat/bmoreno/pvulgaris_relatorio/output/star_stringtie_cuffnorm/genes.fpkm_table",
header = TRUE,
sep= "\t",
dec =".",
stringsAsFactors = FALSE)
```
```
samples.CNTRL <- colnames(pvulgaris.df)[ grep("^CNTRL\\_", colnames(pvulgaris.df)) ] # pegando os nomes dos controles a partir do data.frame
samples.TREAT <- colnames(pvulgaris.df)[ grep("^TREAT\\_", colnames(pvulgaris.df)) ] # pegando os nomes dos tratamentos a partir do data.frame
pvulgaris.samps <- c(samples.CNTRL, samples.TREAT)
dim(pvulgaris.df) #união dos nomes em um só data.frame.
```
Para acessar os dados de log2foldchange das diferentes condições:
```
#dados de log2foldchange do cuffdiff
genes_diff.df <- read.delim('/data/bioaat/bmoreno/pvulgaris_relatorio/output/star_stringtie_cuffdiff/gene_exp.diff',
sep='\t')
dim(genes_diff.df)
```
```
signif.genes_diff.df <- subset(genes_diff.df, q_value <= 0.05) #seleção de apenas os genes diferentemente expressos com q-valor<0,05
signif.pvulgaris.df <- merge(x=signif.genes_diff.df,
y=pvulgaris.df,
by.x='test_id',
by.y='tracking_id') #União dos valores dos genes diferentemente expressos com os nomes de controle e tratamento
dim(signif.pvulgaris.df)
```
```
#solução low-throughput (hard coding), mas é o que consigo pensar agora para arrumar os nomes dos genes
genes <- c("PHAVU_010G142600g(XM_007135533.1)", "PHAVU_010G143200g(XM_007135547.1)", "PHAVU_010G143600g(XM_007135551.1)")
signif.pvulgaris.df ["gene"] <- genes
```
Para exportar este dada frame em forma de .xlsx:
```
library('xlsx')
#Exportar uma planilha com os dados separados anteriormente
write.xlsx( signif.pvulgaris.df[, c('gene','sample_1','value_1',samples.CNTRL,
'sample_2','value_2',samples.TREAT,
'log2.fold_change.', 'q_value')],
file="/data/bioaat/bmoreno/pvulgaris_relatorio/signif_genes.xlsx",
row.names=FALSE)
```
Planilha .xlsx:
```
library('pheatmap')
dfh<-data.frame(group=as.character(pvulgaris.samps))
dfh$group <- ifelse( unlist(
lapply( as.character(pvulgaris.samps) , function(x) { return(regexpr("^TREAT", x)[1]==1); } )
)
, "Tratamento", "Controle" )
dfh$group <- as.factor(dfh$group)
rownames(dfh) <- as.character(pvulgaris.samps)
rownames(signif.pvulgaris.df) <- signif.pvulgaris.df$gene
png(filename ="/data/bioaat/bmoreno/pvulgaris_relatorio/pvulgaris_heatmap.png",
res=300, width=3000, height=3000)
ann_colors <- list(
group=c('Tratamento'='darkgreen', 'Controle'='orange')
)
pheatmap(as.matrix(signif.pvulgaris.df[,pvulgaris.samps]),scale="row",
annotation_col = dfh,
color=colorRampPalette(c("darkblue", "white", "red"))(50),
annotation_colors = ann_colors,
annotation_names_row = TRUE,
show_rownames = TRUE,
cluster_rows = FALSE
)
dev.off()
```
Por fim, o heatmap final da análise feita em aula.
---
# 8. Outras opções
Há outra opção de alinhador que é utilizar desde o começo o Cufflinks. Para isso seus parâmetros estão descritos abaixo:
```
cufflinks --output-dir ${outdir}/${aligner}_cufflinks/${name} \
--num-threads ${THREADS} \
--GTF-guide ${refgtf} \
--frag-bias-correct ${refseq} \
--multi-read-correct \
--library-type fr-firststrand \ #sequenciamento fita específica
--frag-len-mean 400 \ #tamanho médio das reads deste conjunto de dados
--frag-len-std-dev 50 \
--total-hits-norm \
--min-isoform-fraction 0.25 \
--pre-mrna-fraction 0.15 \
--min-frags-per-transfrag 10 \
--junc-alpha 0.001 \
--small-anchor-fraction 0.08 \
--overhang-tolerance 0 \
--min-intron-length 85 \ #alterado com valores do script checkminmaxintronsize.sh
--max-intron-length 1200 \#alterado com valores do script checkminmaxintronsize.sh
--trim-3-avgcov-thresh 5 \
--trim-3-dropoff-frac 0.1 \
--max-multiread-fraction 0.75 \
--overlap-radius 50 \
--3-overhang-tolerance 600 \
--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
```
Para as reads paired-end
```
tophat2 --min-anchor 8 \
--min-intron-length 85 \ # alteração do valor baseado no resultado
do script checkMinMaxIntronSize.sh
--max-intron-length 1200 \ ## alteração do valor baseado no resultado
do script checkMinMaxIntronSize.sh
--max-multihits 20 \
--transcriptome-max-hits 10 \
--prefilter-multihits \
--num-threads ${THREADS} \
--GTF ${refgtf} \
--transcriptome-index ${outdir}/tophat_index/transcriptome \
--mate-inner-dist 100 \ #mais flexível a distância entre as reads para im mesmo transcrito
--mate-std-dev 50 \
--coverage-search \
--microexon-search \
--b2-very-sensitive \
--library-type fr-firststrand \ #alterado por ser fita específica
--output-dir ${outdir}/tophat_out_pe/${name} \
--no-sort-bam \
${outdir}/tophat_index/genome \
${r1} \
${r2} > ${outdir}/tophat_out_pe/${name}.log.out.txt \
2> ${outdir}/tophat_out_pe/${name}.log.err.txt
```
Para as reads singletons:
> Singletons são reads que perderam seu par durante o processo de limpeza e trimagem, mas mesmo assim são utilizadas na análise.
```
tophat2 --min-anchor 8 \
--min-intron-length 85 \
--max-intron-length 1200 \
--max-multihits 20 \
--transcriptome-max-hits 10 \
--prefilter-multihits \
--num-threads ${THREADS} \
--GTF ${refgtf} \
--transcriptome-index ${outdir}/tophat_index/transcriptome \
--coverage-search \
--microexon-search \
--b2-very-sensitive \
--library-type fr-unstranded \ #aqui tem que deixar assim, por ser uma R1 ou R2 que ficou sem par no processo de trimagem.
--output-dir ${outdir}/tophat_out_se/${name} \
--no-sort-bam \
${outdir}/tophat_index/genome \
${outdir}/tophat_out_se/${name}/singletons.fastq \
> ${outdir}/tophat_out_se/${name}.log.out.txt \
2> ${outdir}/tophat_out_se/${name}.log.err.txt
```
O tophat resume o resultado do alinhamento de maneira escrita como a seguinte:
Na imagem acima é um exemplo do tratamento B3, todas as informações sobre o alinhamento de suas reads.
## Diferenças entre os programas
## A. Alinhamento feito com o Tophat2 vs STAR
Utilizando o TopHat2, cerca de 16% (15,98%) em média das reads pareadas alinharam de forma apropriada.
O documento com essas informações está localizado no diretório ./output/align_out_info. Cada condição têm suas informações resumidas nos arquivos de extensão .tophat.out.txt.
Um exemplo de saída com essas informações está representado abaixo:
O STAR, por sua vez, apresentou uma média superior, cerca de 63% de reads alinhadas.
Os arquivos com essas informações estão, também localizados no diretório ./output/align_out_info:
A saída das informações apresentam a mesma estrutura da anterior:
Assim, podemos verificar que o alinhador STAR tem uma performance superior com os dados trabalhados em aula. Este comportamento deve ser observado na definição mais ao final dos genes diferentemente expressos, os resultados que utilizaram o STAR como alinhador, devem ser mais precisos na determinação dos DEGs.
## B. Diferença entre genes expressos
Levando em consideração as abundâncias determinadas inicialmente, era para apenas o gene PHAVU_010G142900g não ser classificado como DEG (gene diferentemente expresso), por apresentar abundância 0,5 em ambas condições.
## Arquivo de saída STAR>StringTie
O arquivo bruto de saída dos genes diferentemente expressos:
Com essa configuração de programas houve apenas um falso positivo. O transcrito terminado em 47.1 não era para ser apontado como DEG, uma vez que sua abundância é igual no controle e tratamento.
## Arquivo de saída STAR>Cufflinks
O arquivo bruto de saída dos genes diferentemente expressos:
O mesmo comportamento é visto utilizando o Cufflinks, o gene PHAVU_010G143200g (com transcrito finalizado em 47.1) é um falso positivo.
## Arquivo de saída TopHat>Stringtie
O arquivo bruto de saída dos genes diferentemente expressos:
O mesmo comportamento é visto anteriormente, o gene PHAVU_010G143200g (com transcrito finalizado em 47.1) é um falso positivo. Porém, com esta combinação apenas três genes foram testados estatisticamente, em que dois deles foram considerados DEGs. Apesar de ter sido testado estatisticamente, o transcrito 47.1 é um falso positivo.
## Arquivo de saída TopHat>Cufflinks
O arquivo bruto de saída dos genes diferentemente expressos:
Neste caso também três transcritos foram testados estatisticamente, e apenas dois foram considerados DEGs. Esta combinação de programas foi a única que não identificou o transcrito 47.1 como DEG. Contudo, outros genes que seriam considerados DEGs não foram identificados.
--
De forma geral os programas apresentam resultados similares, ainda que eles classificam falsos positivos.
## Arquivos dos **output**
Grande parte dos arquivos foram comendados acima em cada etapa. Esta imagem é apenas um resumo de todos os diretórios que ao final ficam no diretório output.
# Detalhes importantes
<details open>
Foi utilizado o programa IGV para ver o alinhamento dos genes com o genoma de referência depois do alinhamento utilizando o programa bowtie2.
* Vimos também brevemente a importância de testar o dm5, para ver se os dados recebidos da facility não estão corrompidos.
* O nome das pastas e arquivos tem a sequência dos programas utilizados. Dessa forma fica mais fácil de manter a sequência de cada programa utilizado.
* O valor de p-adj ou q-value é uma representação dos múltiplos testes. A fim de deixar o estudo com um erro global de alpha = 0.05, cada cada gene quando testado hipótese se ele é diferentemente expresso ou não não pode ter um valor de significância alpha = 0.05, ele tem que ser mais baixo. Para fazer essa correção de múltiplos testes, usamos o FDR (*False Discovery Rate*) de 0,05, que indica que eventualmente haverá falsos positivos. A ideia de maineira bem simplificada é corrigir os p-valores d efrma que apenas 5% das vezes pode ser que haja falsos positivos identificados.
* Todos os comentários (#) dentro do scritp bash estão apenas nestas anotações. Diferentemente do R, os comntários são lidos e dão problemas no seguimento da leitura do script;
* No script colocamos alguns parâmetros que não existiam antes. Em sua grande maioria é auto explicativo, por isso se optou usar o nome mais descritivo do parâmetro, ao invés da opção de um ou duas letras para descrição do parâmetro.
Sign in with Wallet
Connect another wallet