# **Tutorial de Bioinformatica: Análise de trascriptomas**
A seguir será apresentado o relatório de atividades de simulação, referentes à análise de transcriptomas, utilizando o sistema operacional Linux/GNU, para o uso das ferramentas e programas associados a atividade foi feita através de conexão remota (SSH) com o servidor Hammer, localizado no Departamento de Tecnologia da Faculdade de Ciências Agrárias e Veterinárias da UNESP.

Todos os scripts utilizados na etapas seguintes foram copiados do diretório /data/tmp/bioaat para a pasta criada para o projeto no diretório: /state/partition1/<usuário>/

### **Parte 1 - Criando e ajustando as referências**
Primeiro foi selecionado um organismo para ser feita a simulação, no meu caso escolhi uma bactéria causadora de enfermidade em peixes a *Flavobacterium columnare*, o primeiro passo foi baixar o arquivo fasta do genoma referência e o arquivo "gff", este arquivo contém informações referentes a anotação e posição dos genes no genoma referência, para baixar estes arquivos primeiro foi criada uma pasta no diretório /state/partition1/darferreira com o nome do projeto no meu caso criei uma pasta denominada FlavobacteriumRNAseq, o comando para a criação de pasta pode ser feito através de caminho absoluto ou relativo no meu caso eu entrei na pasta do meu usuario e utilizei o seguinte comando.

Entrar na pasta:

cd /state/partition1/darferreira

mkdir ./FlavobacteriumRNAseq


Dentro do diretório criado para o projeto foi criada pasta chamada "Refs" onde foi depositado os arquivos referência.

mkdir ./Refs

No diretório Refs `cd ./Refs` foi utilizado o comando wget para baixar os arquivos referência, para isto é necessário o link do protocolo (fttp:\\)) do site da NCBI, os passos a seguir foram realizados.

1- Entrar no site do GenBank através do endereço na barra de navegação (Foto), escolhendo a opção Genome e digitando na caixa de busca o organismo interesse:


2- No canto superior esquerdo do conteúdo da busca, clicar com o botão direito do *mouse* sobre a opção "genome" e selecionar a opção copiar link (o link copiado já será do protocolo fttp), segue foto:


3- Com o link do genoma copiado ir até o terminador e digitar na linha de comando o comando wget (clicar com o lado esquerdo do *mouse* que funcina como a opção colar), executar o comando na pasta Refs, deste modo será baixado o arquivo fasta do genoma do organismo.

4- Repetir os passos 2 e 3 clicando na opção "gff" ao invés de "genome"

Agora será nescessário descompactar os arquivos baixados pelo comando wget, os arquivos estarão em formato .gzip ou seja estão compactados, para descompactar utilizar o comando gunzip e o nome do arquivo, automaticamente o arquivo será descompactado na pasta que foi executado o comando.

Um arquivo com o codigo de acesso da referência aparecera na pasta um com a extensão .fna (genoma) e outro com a extensão .gff (genome features), para facilitar o trabalho os arquivos foram renomeados com o nome genome.fna e genome.gff, para isto foram utilizados os comandos abaixo:

mv NOMEORIGINALdoARQUIVO.fna ./genome.fna
mv NOMEORIGINALdoARQUIVO.gff ./genome.gff

Após isto foi feita a limpeza do cabeçalho do arquivo do genoma, abaixo segue a foto do exemplo que utilizei.


Para fazer a limpeza do baceçalho e deixar apenas o numero de assesso na linha 1 e abaixo a sequencia foi utilizado o comando "sed" conforme abaixo:

sed 's/^>\([^\.]\+\).*/>\1/' genome.fna > genome.fa

Deste modo o conteúdo do arquivo gerado deve ser este:


Agora será necessário "ajustar" o arquivo gff baixado do GenBank para isto será utilizado o script "[fixNCBIgff.sh](https://https://github.com/dgpinheiro/bioinfoutilities)" que pode ser baixado através do repositório Github do docente Daniel Guariz Pinheiro, lembrando que este script possui outros três scripts dependentes que deverão ser baixados para a execução do script fixador:

Scripts necessários para execução.
- fixGFFgenestruct.pl
- fixGFFdupID.pl
- sortGFF.pl

Não será necessário baixar estes scripts pois no servidor hammer eles já estão presentes e estão com o caminho registrado na variável $(PATH).

Executar então o script fixNCBIgff.sh e direcionar como "in file" o arquivo genome.gff

fixNCBIgff.sh genome.gff

Após isto o arquivo receberá adição do nome genome_fixed.gff, deste modo o arquivo original ja pode ser removido, e o arquivo trabalho pode ser renomeado.

rm genome.gff
mv genome_fixed.gff ./genome.gff

Utilar o comando gffread para a transformação de arquivo gff para arquivo ."gtf" (*genome transcriptional features*), este formato trará detalhadamente as CDS do genoma sem trazer as anotações, (**CURIOSIDADE ADICIONAL**) este comando também constroi os arquivos trascriptome.fa e proteome.fa estes arquivos são o transcriptoma e proteoma em formato fasta para isto executar os comamndos a seguir:

gffread genome.gff -g genome.fa -T -o genome.gtf
gffread genome.gff -g genome.fa -w transcriptome.fa
gffread genome.gff -g genome.fa -y proteome.fa

### Parte 2 - Criando arquivos de abundância e seleção dos genes para a simulação.

Com os arquivos referência já disponiveis, acessar o arquivo gff através do comando less

less genome.gff

Visualizar o seguinte aquivo, ((Nota: todas as anotações estão separadas por tab))


As informações assinaladas serão utilizadas para compor um novo arquivo, este modelo será o utilizado para procariotos, para eucariotos cada gene completo sem introns tem uma referência, portanto será utilizada a referência interna do gene, como o organismo alvo é uma bactéria, as descrições da foto acima foram anotados para a criação de um arquivo conforme abaixo:

nano ACCS.txt 

Para facilitar este processo o comando screen é bem útil pois é possível selecionar os elementos que deseja-se copiar para o arquivo texto em uma janela e na outra colocar o aquivo texto, este arquivo texto gerado será utilizado para fazer a simulação, para isto foram selecionados 6 genes (para não sobrecarregar o servidor) no caso da realização com eucariotos será necessessário além de selecionar os genes criando arquivo texto com os numeros de acesso dos genes, fazer a limpeza do genoma e do arquivo gff, selecionando apenas o cromossomo correspondente ao gene e criando um arquivo chamado toygenome.fa e toygenome.gff.

Abaixo segue exemplo do arquivo contendo o codigo de acesso do gene, codigo de acesso do genoma, posição iniciadora, posição terminadora e se é da fita senso ou antisenso referendo-se como plus fita senso e minus fita antisenso.


~~**~~***NOTA: DEVERÁ SER SUBTRAÍDO UMA BASE DO CÓDON DE INICIAÇÃO E TERMINAÇÃO IDENTIFICADOS NO ARQUIVO GFF E TODAS AS COLUNAS DEVEM SER SEPARADAS POR TABULAÇÃO***~~**~~ instruções sedidas pelo script fixador do NCBI

Criando estimativas de abundancia que devem ser nomeados de abundance_A.txt e abundance_B.txt, que serão referentes ao organismo na condição A e B (**NOTA: OS ARQUIVOS SEMPRE DEVEM SEGUIR EXATAMENTE O MESMO NOME POIS NA HORA DE RODAR OS SCRIPTS JÁ ESTA REGISTRADO ESTES NOMES**)



Os arquivos gerados com estimativas de abundância para cada gene, deverão ter diferenças hipotéticas; Por exemplo pelo menos um gene que esteja altamente expresso na condição A e menos na condição B, a soma dos valores de abundância de A e a soma dos valores de abundancia de B devem ser 1 respectivamente.

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

Simulando a geração de transcritos através do script sim-prok.sh para procariotos e sim.sh para eucariotos (https://github.com/dgpinheiro/bioaat/blob/master/sim-prok.sh) por padrão este script gerará 2500 fragmentos aleatórios com média de 300pb (+/- 30pb):


O *script* criará uma pasta chamada raw e depositará os transcritos simulados, com 3 réplicas para cada condição (A e B), caso queira modificar basta alterar os parametros no script.

### Parte 3 - Pré processamento dos transcritos gerados (corte de adaptadores e trimagem das sequencias)

Ajustar os parâmetros do script preprocess5.sh, primeiramente deve-se copiar o script da pasta biaat para a pasta do projeto:

cp /data/tmp/bioaat/preprocess5.sh /state/partition1/darferreira/FlavobacteriumRNAseq

Para procariotos pode-se alterar o parametro de trimagem pelo programa prinseq encontrando o ponto do script conforme abaixo através do comando nano ./preprocess5.sh:


Os parâmetros trim_tail_right 5 e trim_tail_left 5, os valores podem ser alterados para 0, este parâmetro se refere a trimagem de poli-A tanto no transcrito R1 como em R2, portanto não se aplica ao organismo escolhido. Depois de finalizados os ajustes procurar por demais ajustes e apertar CTRL X para sair do editor, o programa perguntará se quer salvar as alterações, precione y e ENTER
Caso queira adicionar ou alterar algum parâmetro consultar o manual do [PrinSeq](https://https://vcru.wisc.edu/simonlab/bioinformatics/programs/prinseq/prinseq-lite.help.txt) e [Atropos](https://https://atropos.readthedocs.io/en/latest/guide.html).
Agora iremos criar os diretórios de input e output que serão utilizados no pipeline, ir para o diretório do projeto. `cd /state/patition1/daferreira/FlavobacteriumRNAseq` criar os diretórios input e output:

mkdir ./input \
mkdir ./output

Criar links simbólicos dos transcritos gerados (diretório raw) no diretório input, isto pode ser feito executando o script mklink.sh que pode ser copiado do diretório bioaat localizado no caminho /data/tmp/bioaat neste caminho encontram-se todos os scripts utilizados no curso, de todo modo aconselha-se a execução do script mklink.sh para certificar da exatidão dos nomes dados aos links simbólicos.


**NOTA: EXEUTAR O SRIPT mklink.sh NA PASTA DO PROJETO ONDE TEM O DIRETÓRIO raw E AS DUAS PASTAS CRIADAS input E output**.

Deste modo o conteúdo da pasta input ficara da seguinte maneira:


Para conferir se os links foram criados corretamente pode-se utilizar o comando:

ls -lh /state/partition1/darferreira/FlavobacteriumRNAseq/input

O conteúdo deve aparecer da seguinte maneira:


Agora sim, podemos executar o script preprocess5.sh (lembrando que este scrip ja esta embutino no script rnaseq-ref-prok.sh) portanto ao executar o script ja utilizara os parametros setados a seguir.

cd /state/partition1/darferreira/FlavobacteriumRNAseq

preprocess5.sh ./input ./output

Será criado o diretório processed dentro da pasta output, tendo como conteúdo os diretórios cleaned, atropos, fastQC e prinseq todos contendo a pasta pre e pós gravando os resultados das análises com os respectivos programas.

**Observação**: Para vizualizar os arquivos fastq pelo programa FASTQC, realizar o download de um programa de gerenciador de pastas por conxão ssh, eu utilizei o programa WinSCP, nele vc pode conectar ao hammer e tranferir arquivos do servidor para o computador pessoal, o conteúdo da pasta ./output/processed/fastqc pode ser transferido para ser visualizado em opção html:// conforme a baixo:

Exemplo de arquivo do pré processamento:


Exemplo de arquivo do pós processamento:


**obs. mesma sequencia nos dois exemplos pode-se observar instantaneamente a melhoria da qualidade das reads.**

A pasta cleaned no caminho /state/partition1/darferreira/FlavobacteriumRNAseq/output/processed contém links simbólicos do resultado do programa prinseq que faz a trimagem portanto a ultima etapa do preprocessamento, é com estes arquivos que seguiremos para as etapas seguintes.

**IMPORTANTE: O SCRIPT ACIMA ESTA CONCATENADO AO SCRIP rna-seq-prok.sh PORTANTO QUANDO FOR EXECUTADO O SCRIPT ELE EXECUTARA O SCRIPT DE PRÉ PROCESSAMENTO BASTA ALTERAR OS PARÂMETROS CONFORME SUA ESCOLHA E DEPOIS SALVAR**
### Parte 4 - Alinhamento dos transcritos utilizando os softwares Star e TopHat2 e montagem do transcriptoma com Cufflinks

Foi utilizado o script [rnaseq-ref-prok.sh](https://github.com/dgpinheiro/bioinfoutilities) para realizar o alinhamento (STAR e TOPHAT2) e cálculos de expressão diferenciada e montagem do transcriptoma (CUFFLINKS) em forma de pipeline, automaticamente este script executará o script preprocess5.sh, portanto basta auterar os parametros do script conforme a sua necessidade, mas antes vamos conferir os parametros que estão embutidos no scrip rnaseq-ref-prok.sh para isto abrir o scrip com editor de texto "vim" ou "nano" e conferir o bloco onde os parametros são adicionados.

Os argumentos da linha de comando para a execução do script são os seguintes:
O **argumento 1** do script será o alinhador que deseja-se utilizar, o **argumento 2** será diretorio input que no caso é a pasta input que contém os links simbólicos gerados anteriormente, o **argumento 3** diretório de saida que pode ser ./output, **argumento 4** numero de trheads que serão utilizadas, neste caso será no máximo 6 para o servidor que estamos utilizando recomenda-se sempre colocar uma thread a menos para permitir que funções básicas da maquina continuem operando, **argumento 5** arquivo GTF (gerado pelo gffread), **argumento 6** arquivo FASTA do genoma (cabeçalho limpo), **argumento 7** que é qual montador deseja utilizar no caso cufflinks, **argumento 8** que é o arquivo de contaminantes caso não tenha indicar a opção **NA**, **argumento 9** arquivo do protocolo fttp da NCBI que contém as anotações referentes a cada gene de cada oranismo depositado, este arquivo tra uma coluna chamada Tax_id esta coluna contem um indentificador por espécie no qual todos os genes da espécie terão este identificador na linha (este arquivo tem varias versões detalhes serão explicados posteriormente *see*.Parte 8) e por ultimo **argumento 10** que deve ser oo numero do taxonomy id da espécie que se esta trabalhando (*see*.Parte 8)conforme exemplo abaixo:

##### Alterando parâmetros

Para realizar a alteração dos parâmetros no script, basta buscar o ponto do documento no qual se aplicam os parametros, sendo fundamental o uso dos respectivos manuais para se identificar quais são os parametros *defalt*, e o quais são as funções.
Vamos conferir os parametros que estão embutidos no scrip para isto abrir o scrip com editor de texto "vim" ou "nano" e conferir o bloco onde os parametros são adicionados.

**Importante: Caso muitas reads não sejam mapeadas mas tem correspondencia no genoma, pelo alinhador TopHat recomenda-se um processamento extra para estas reads, recomenda-se a leitura do artigo** https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1058-x.

### Parte 5 - Gerando Status do alinhamento para comparação de eficiência

Para a melhor vizualização do resultado final do pré processamento e alinhamento dos transcritos, foi criado pelo docente **Daniel Guariz Pinheiro** o script [stats_preprocess.sh](https://github.com/dgpinheiro/bioinfoutilities). Para a execução deste script o **primeiro argumento** é o diretório **input**, que será o conteúdo do diretório **tophat_out_final** ou **star_out_final** que serão os resultados do alinhamento feito com os dois programas e o segundo argumento o nome do arquivo destino, abaixo segue exemplo do diretório output após a realização do alinhamento dos transcritos com os dois programas:


Após a realização do alinhamento, pode-se utilizar o comando cat para comparar os resultados das duas análises, abaixo segue o **exemplo 1** que foi gerado com o resultado dos dois alinhadores TopHat2 e Star.

**exemplo 1**

Para conferir quais foram os parâmetros utilizados na análise que se quer comparar basta ir ao diretório output:

**<alinhador>_out_final/<amostra>/Header txt**

Exemplo de arquivo Header extraído da réplica um condição controle, mas vale para todos pois a análise sempre é feita igualmente para todoas as réplicas e condições


**^^O arquivo acima apresenta status dos resultados dos parâmetros utilizados na primeira amostragem utilizando o programa STAR**

Para o segundo teste utilizando o alinhador Star, só que alterando o parâmetro <outFilterMismatchNoverReadLmax 0.01> para 0.04, lembrando que os parâmetros devem ser alterados no script rna-seq-prok.sh, abaixo segue o resultado dos status do alinhamento após a alteração do parâmetro.


Comparando ao resultado anterior **exemplo 1** pode se observar maior proporção de reads mapeadas, alcançando um maior aproveitamento das reads sequenciadas.

Para escolha dos parâmetros ler o manual do alinhador que deseja-se utilizar [TopHat2](https://https://ccb.jhu.edu/software/tophat/manual.shtml) ou [Star](https://https://physiology.med.cornell.edu/faculty/skrabanek/lab/angsd/lecture_notes/STARmanual.pdf)

Recomenda-se sempre anotar os parametros que foram alterados inclusive no script rna-seq-prok.sh alguns valores foram alterados para uma maior eficiência no alinhamento de sequencias de mRNA de procariotos, sempre ter em mão quais parametros foram alterados da condição *default*, de todo modo aconselha-se rodar os programas com os parâmetros base e observar os resultados, e ir alterando comforme nescessidade.
O manual para o uso do cufflinks se encontra disponível em http://cole-trapnell-lab.github.io/cufflinks/cufflinks/.

### Parte 6 - Visualizando os resultados do alinhamento utilizando o software IGV

Para a vizualização dos resultados utilizando interface gráfica,foi utilizado o [IGV (Integrate Genome Viwer)](https://https://software.broadinstitute.org/software/igv/download), a seguir será mostrado teste para a vizualização dos resultados, para a amostra controle replica 1.

Primeiramente será necessário gerar um arquivo indice para o resultado **aligned.sorted.bam** o arquivo ganhará a extenção .bai, para gerar o indice ir até o diretório que contém os resultados para esta amostra.

cd /state/partition1/darferreira/FlavobacteriumRNAseq/output/tophat_out_final/CTR01

samtools index ./Aligned.sorted.bam

Depois de gerar o indice, será necessário tranferir para o computador pessoal os arquivos genome.fa (cabeçalho limpo), genome.gtf (gerado pelo comando gffread), o arquivo Aligned.sorted.bam e o arquivo Aligned.sorted.bam.bai (gerado pelo samtools, passo acima), para o computador pessoal, eu utilizei o programa WinSCP.

Após tranferir os arquivos abrir o programa IGV, ir na aba Genomes e selecionar o arquivo genome.fa, depois na aba "File" carregar o arquivo genome.gtf e Aligned.sorted.bam, o indice apenas prescisa estar na mesma pasta automaticamento o programa reconhecerá


### Parte 7 - Vizualização dos valores normalizados (cuffnorm) e valores brutos

O script rna-seq-ref-prok.sh contém nos blocos finais do script, a execução do script em "R" chamado [de-nomarlise-cuffnorm.R](https://https://github.com/dgpinheiro/bioinfoutilities/blob/master/de-normalize-cuffnorm.R) este script irá utilizar o arquivo de saída do programa cuffnorm "isoforms.count_table" para obter os dados brutos (quantas reads mapeadas por locus) e organizar o resultado no arquivo "isoforms.raw_count_table.txt" no mesmo diretório.


O diretório output/<alinhador>_<montador>_cuffcompare/ vai conter o arquivo TCONS-nearstref.txt nele contém a associação TCONS (id do processo) e nearst ref que é o gene mais proximo ou onde o transcrito alinhou, essa informação foi obtida pelo programa através dos arquivos genome.gff e transcriptome.fa.

exemplo do arquivo que contém as associações TCONS (interno) nearest_ref identificado do gene depositado na NCBI, este arquivo tembém é utilizado como auxilo para a fusão de tabelas e fusão da anotação (see Parte8 e 9) pelos scripts seguintes.

#### Erro identificados no processo de alinhamento

Na foto é possível observar que os genes HUE53_RS00120 e HUE53_RS00125 ficaram com os mesmos valores de expressão é possivel observar este fato vizualizando o arquivo gene_exp.diff para facilitar pode-ser utilizado o comando "grep", isto aconteceu devido a estes genes compartilharem sequencia codificadora. Sobretudo no resultado final da expressão contido no arquivo isoforms_exp.diff pode se notar que apenas o gene correto está presente.

### Parte 8 - Vizualização dos valores de expressão diferencial (Cuffdiff)

O último bloco do script rnaseq-prok-ref.sh utiliza dois scripts em "R", que organizarão os resultados em um arquivo único, facilitando a leitura conforme foto abaixo:

O primeiro script [splitteR.R](https://https://github.com/dgpinheiro/bioinfoutilities/blob/master/splitteR.R) irá utilizar o arquivo de saída do programa cuffdiff, localizado na pasta output/**<alinhador>_<montador>**_cuffdiff/gene_exp.diff, este arquivo contém os valores de expressão normalizados para cada tratamento e traz também os genes expressos de acordo com informações do arquivo "gtf", entretanto quando mais de um gene é assimilado ao resultado de expressão devido ao fato da proximidade em alguns loccus gênicos especialmente em procariotos, o "symbol" dos genes relacionados aparecem na mesma coluna separados por "," então o script apenas duplicará a linha colocando cada gene relacionado aquele valor de expressão separados e duplicará as informações contidas no restante das colunas segue foto indicativa do processo.

Arquivo inicial "gene_exp.diff", os genes testados (simulados) traram em sua linha na coluna teste a descrição "OK", então para facilitar a vizualização foi utilizado o comando grep 'OK' gene_exp.diff para mostrar apenas os genes que foram simulados

**NOTA: genes na marcação foram agrupados na mesma linha separados por ","**

arquivo gerado pelo script splitteR.R "gene_exp.diff.splitted.txt"
**NOTA: cada gene separado por linha**

A última estapa irá executar o script [mergeR.R](https://https://github.com/dgpinheiro/bioinfoutilities/blob/master/mergeR.R) que irá unir as informações do arquivo "gene_exp.diff.splitted.txt" e as desrições dos genes e anotação contidos no arquivo "gene_Id" que traz todos os genes da maioria dos taxons (procurar a versão do NCBI que contém o taxon de interesse, caso o organismo que esteja trabalhando não esteja no arquivo baixado pelo protocolo FTTP da NCBI) para baixar o conteúdo das referencias do NCBI e anotação utilizar o comando wget e colar o link do protocolo "fttp" que contenha o taxon de interesse, para encontrar as informações desejadas o script utilizará os dois ultimos argumentos da linha de comando utilizada para executar o script rnaseq-ref-prok.sh que serão o local do arquivo contendo as informações das anotações dos genes com uma coluna relacionada ao "tax_id" no meu caso eu busquei o taxonomy id do organismo que utilizei para simulação atravéz da busca pelo nome da espécie no site taxonomy da NCBI.
Ao referênciar o arquivo contendo as informações seguido do "taxonomy id" o script ira filtrar as anotações referentes aquele taxon do arquivo que no caso é o "gene_info" e criará ou arquivo na pasta output chamado gene_info <taxonomyID> com as informações filtradas ou seja todos os genes e anotações da espécie que esta trabalhando.
Abaixo segue informações de como baixar o arquivo gene info através do [protocolo fttp da NCBI](https://https://ftp.ncbi.nlm.nih.gov/genbank).
**NOTA: Entrar em um dos links no endereço eletrônico apresentado na foto.**

Desta forma o script mergeR.R gerará o arquivo gene_exp.diff.annot.txt que trará as informações de expressão compiladas a anotações contidas no arquivo gbk, o conteúdo deste arquivo é separado por tabulação portanto pode ser tranferido para o computador pessoal e aberto pelo exel que disponibilizará as informações nas células, estas informações compiladas contém os valores de expressão normalizados o gene, locus gênico, valor de "Q" e "P" teste de significancia dentre outros conforme imagem abaixo, caso ja queira colocar o arquivo texto em formato de planilha utilizar o script tsv2xls.sh que fará a conversão depois só tranferir para o computador pessoal.

**Passo 1** após transferir o arquivo "gene_exp.diff.annot.txt" para o computador e abrir no excell, perceberá que a tabela trás informação de todos os genes testados e não testados, para facilitar a vizualização ativar a opção filtro e ir até a coluna "status" e selecionar "OK".
Deste modo a tabela ficará apenas com os genes testados (simulados)
Para aumentar as informações da tabela é interessante trazer os valores brutos de expressão então neste caso é só abrir o arquivo "gene_exp.diff.splitted.txt" no computador pessoal da mesma maneira feito anteriormente e organizar os dados em uma tabela conforme abaixo:

Ainda é possível utilizar o script mergeR.R para unir informações a tabela de trabalho deste modo o arquivo resultante pode ser utilizado para analises posteriores (o quanto mais completo melhor), pode ser usado os seguintes argumentos neste caso eu uni as informações de dados brutos de expressão e expressão normalizada unindo os arquivos cuffnorm/gene_raw.count.table e cuffdiff/isoforms_exp.diff.annot.txt (testes de significancia e log2FOLDCHANGE) no arquivo raw_exp.diff.txt utilizando a linha de comando abaixo:


Sempre lembrar de direcionar o indentificador comum com o titulo do cabeçalho onde ele se encontra no caso para "y" foi "test_id" e "x" foi "trakking_id" eles contém o identificador comun "TCONS" e indicar a opção --print.out.label pois ela preservará os cabeçalhos das tabelas originais.

Para trazer também a anotação dos genes eu utilizei o arquivo cuffdiff/genes_exp.diff.annot.txt e executei o script mergeR.R novamente com a seguinte linha de comando:

**NOTA: PARA ESTE SCRIPT OS DOIS PRIMEIROS ARGUMENTOS SÃO OS ARQUIVOS DE ENTRADA APÓS ISTO OS DOIS ARGUMENTOS SEGUINTES SERÃO O IDENTIFICADOR DA COLUNA QUE TERÁ INFORMAÇÃO COMUM NOS DOIS ARQUIVOS DE ENTRADA E POR FIM (OPCIONAL) PODE SELECIONAR A OPÇÃO --print.out.label ESTA OPÇÃO CONCERVARÁ O CABEÇALHO IDENTIFICADOR DAS COLUNAS E POR ULTIMO SELECIONAR O ARQUIVO DE SAÍDA**

Deste modo as informações: significância, gene, valores brutos (de expresão para cada amostra nas duas condições biológicas), valor normalizado por condição biológica, log2FoldChange, valores de "P" e "Q" e anotação relacionada ao gene foram organizadas conforme tabela abaixo:


A tebela acima também traz informações referentes ao arquivo abundance_A e abundance_B que carregam a informação inicial das proporções utilizadas pela simulação, deste modo é possível comparar, se os resultados são fiéis ao presuposto, foi observado que as diferenças iniciais para os genes HUE53_RS00105 e *125 eram iguais e após a simulação os resultados indicaram redução e aumento de expressão, este fato é esperado por se tratar de uma simulação e com numero reduzido de genes, mas ao observar os valores de log2(fold-change) pode se notar que é um valor bem proximo de "0"

**Caso a leitura da tabela esteja dificil pela falta de proporção entre as colunas separadas por tabulação, pode se utilizar a linha de comando sed que ira "normalizar" os espaços conforme abaixo**

column -t -s$'\t' ./output/edger/CTR-TST.txt | less -S

#### **Normalização utilizando os scripts run-DESeq2.R run-edgeR.R**

Após o script de-normaliseR.R gerar os valores brutos de expressão para cada amostra e condição biológica, estes dados podem ser utilizados nos dois scripts acima em 'R' que realizarão o cálculo de normalização levando em conta os valores de TPM (transcripts per Megabase) e tamanho das sequências e realizará também o claculo do log2(fold-change) e testes de significância com o ajuste de FDR (*false discovery rate*).

Primeiramente para o rodar os scripts criar arquivo na pasta do projeto chamado groups.txt.

cd /state/partition1/<usuário>/<projeto>/ \
nano groups.txt

Adicionar seguinte conteúdo:
**LEMBRANDO QUE AS INFORMAÇÕES NAS COLUNAS DEVEM SER SEPARADAS POR TABULAÇÃO E OS CABEÇALHOS DEVEM TER EXATAMENTE OS NOMES ACIMA**

O arquivo de entrada utilizado neste script, são os valores de expressão não normalizados (ou seja o número de reads que mapearam para o gene) estes valores foram gerados pelo script de-normalise-cuffnormR.R no arquivo isoforms_raw.count.table.txt localizado no diretorio <linhador>_<montador>_cuffnorm/

-Criar diretório de saída na pasta output com o nome de sua preferencia selecionei deseq2 e edgeR.

Para que não seja necessário entrar no ambiente R para rodar o script foi utilizado o script run-deseq2.R e run-edgeR.R.

**Caso tenha dúvidas sobre a construção da linha de comando, o manual tanto do script run-DESeq2.R como o run-edgeR.R pode ser acessado apenas digitando o nome da função e a tecla "ENTER"**

`--in=<inputdir> --groups=<diretóriodeGROUPS --out=outdir/<Deseq2 ou edgeR.R --ncores=(numero de cores para a análise para o servidor que utilizamos sempre recomenda-se 6 pois cada computador é compartilhado com dois alunos)`

Exemplo de tabela após fundir todas as informações incluindo a normalização pelos scripts deseq2 e edgeR.R.


Ao observar a tabela nota-se que o resultado do script DESeq2.R apresentou resultados que melhor se assemelham as abundancias inicias, isto se deve pela normalização do valor P com o criterio de falsa descoberta, portanto os genes que tinham niveis inicias de abundancia iguais por mais que apresentem log2FC diferente de zero, este resultado não foi significante:


Para observar o resultado dos scripts em 'R' (deseq e edgeR) utilizar o comando `column -t -s$'\t' ./output/<deseq ou edgeR/CTR-TST.txt | less -S`.


Os valores dentro dos cículos apresentam resultados de reads que foram mapeadas erroneamente, mas pode-se observar que os valores normalizados não ultrapassam 30 e o valor de 'P' é acima de 0,005 portanto estas falhas serão discartadas na expressão diferencial, este mapeamento erroneo de algumas reads se deve ao fato de durante a fragmentação do mRNA alguns pedações apresentam Kmers bem frequentes possibilitando o alinhamento da read em outra região do genoma.

# Montagem "*de-novo*" utilizando o pipeline "Trinity"

O pipeline Trinity conta com 3 programas Ichtworm, Crysalida e Buterfly, estes programas vão realizar o alinhamento e cálculos de expressão diferencial de transcritomas sem referência. Mas antes vamos conferir os parametros que estão embutidos no scrip para isto abrir o scrip com editor de texto "vim" ou "nano" e conferir o bloco onde os parametros são adicionados conforme foto abaixo:

**Lembrete: Conferir o [manual do trinity](https://github.com/trinityrnaseq/trinityrnaseq/wiki) para conferir as opções que se deseja inserir na análise**

Para a realização desta etapa utilizaremos o script rnaseq-denovo.sh copiado do diretório /data/tmp/biaat/
Primeiramente foi criado pasta trinity dentro do diretório output `mkdir ./output/trinity` está pasta sera utilizada para armazenar os resultados do pipeline

Para executar o script foi montada linha de comando conforme abaixo:
TRINITY 

Mas antes vamos conferir os parametros que estão embutidos no scrip para isto abrir o scrip com editor de texto "vim" ou "nano" e conferir o bloco onde os parametros são adicionados conforme foto abaixo:


O argumentos da linha de comando para execução do script são os seguintes:
**Primeiro argumento** é o diretório input que pode ser a saida dos programas prinseq e atropos, no caso indiquei o conteúdo do diretorio processed/cleaned , o **segundo argumento** é o diretório de saída que foi criado, o **terceiro argumento** é o numero de Threads e o **ultimo argumento** é a memoria utilizada no processo remomenda-se 12G.

### Parte 1- Visualisação dos resultados

Para conferir o resultado do pipeline visualisar o conteúdo do diretório output/trinity/trinity_assembled/

Esta pasta contém o arquivo Trinity.fasta que tem como conteudo os *contigs* gerados pelo processo de alinhamento das reads sem referência segue foto do arquivo.


Pode se utilizar o comando `grep '^>' Trinity.fasta` para contar quantos *contigs* foram gerados, lembrando que como o organismo escolhido foi procarioto cada gene deve gerar um contig sem referência.


No meu caso como foram selecionados 6 genes, 6 *contigs* foram observados.

### Parte 2- Anotação dos contigs gerados.
Agora é necessário atribuir identidade as sequencias geradas, como não há referencia neste tipo de análise, será construido um indice que será utilizado como base de dados, como foi feito uma simulação o arquivo (transcriptoma.fa) bastará, portanto os genes presentes estão contidos neste arquivo `grep '^>' transcriptoma.fa`
Então para indexar este arquivo como base de dados para realizar nBLAST na etapa seguinte, a indexação é feita da seguinte maneira:
`makeblastdb -title "ReferenceTranscriptoma" -dbtype nucl -out RefTrans -in transcriptoma.fa`

#### nBLAST
O comando a seguir realizará o alinhamento do arquivo contendo as sequencias em formato fasta com a base de dados indexada, e identificar maior correspondência.
`blastn -max_hsps 1 -max_target_seqs 1 -query ./output/trinity/trinity_assembled/Trinity.fasta -db ./RefTrans -out ./Trinity_X_RefTrans.txt -evalue 1e -5 -outfmt "6 qseqid pident lenght mismatch gapopen qstart qend sstar send evalue bitsocore"`

-max_hsps: numero maximo de "hits" ou seja mostrar resultados com alto score de alinhamento
-max_target_seqs: numero de correspondentes que deseja mostrar
-query: sequencia(s) que se deseja realizar a consulta formato fasta
-db: arquivo com a base de dados indexada
-out: arquivo de saída para análise
-outfmt: consultar por -outfmt em `blastn -help`

Deste modo o arquivo final gerado pela consulta apresentará o seguinte formato:

**NOTA: são os mesmos genes selecionados pela simulação**

##### Verificando se as sequecias são codificadoras

Para verificar se as sequencias montadas são codificadoras utiliza-se a uma ferramenta do pacote trinity chama TansDecoder, a execução é feita da seguinte maneira:

TransDecoder.LongOrfs -t ./output/trinity/trinity_assembled/Trinity.fasta -m 90

O parâmetro m referese ao numero minimo de aminoácidos das sequencias ao final da execução da linha de comamando indicará o diretorio de saída.
Este comando irá transcrever os contigs em sequencias de proteínas gerará o arquivo trinity.fasta.transdecoder.pep (visulização do arquivo abaixo).

##### Predição
Agora será necessário fazer a predição das proteínas (modo alternativo realizar um blastp) mas do modo apresentado neste tópico a predição será feita pelo comando TransDecoder.Predict

TransDecoder.Predict -t <trinity_assembled>/Trinity.fasta -m 90 --single_best_only

Para obter ajuda sobre os parâmetros que pode se passar para esta execução digitar TransDecoder.Predict e aparecerá os opções, caso tenha feito blastp anteriormente pode passar esta opção, passar o argumento --single_best only (apenas a melhor correspondência).

O arquivo de saída irá para o diretorio que será criado pelo comando (caso queira refazer a análise utilizar a opção rm -fr TransDecoder.dir*).
##### Estimatia de abundancias (expressão)
Para observar o resultado do comando checar o arquivo Trinity.fasta.transdecoder.pep. Agora criar o diretóório abundance dentro de output/trinity, copiar o script rnaseq-trinity-abundance.sh da pasta bioaat.
Este script após gerar as estimativas de abundancia fora a nomalisação e os testes de expressão diferencial baseados no log2fold-change com o script DESeq2.R que está embutido no script de análise de abundancia.
O arquivo de saída contendo as estimativas de abundância para cada contig, amostra e condição biológica estará no diretório 'output/trinity/abundance/CTR-TXT.txt'

Abaixo segue foto do resultado com a visualização do arquivo CTR-TST.txt


### Anotação das sequencias de aminoácidos geradas com EggNogMapper e InterProScan
O programa EggNogMapper conta com uma base de dados curada de proteínas esta base será utilizada para inferir anotação as proteínas preditas na etapa anterior através da busca por similaridade em genes ortólogos.
O InterProScan realizará anotação funcional dos peptideos através da identificação do domínio e relação com informações do Gene Ontology.
Primeiro rodaremos o programa EggNoggMapper e para isto ele necessita de um ambiente python ativado, executar o comando <pyenv activate eggnog-mapper>.
Aparecerá no identificador fixo do usuario a descrição eggnog-mapper isto significa que ja esta em um ambiente python especifico para rodar o programa agora para realizar a anotação o software executar a linha de comando:

emapper.py -m diamond --no_annot --no_file_comments --cpu 6 -i Trinity.fasta.transdecoder.pep --dmnd_iterate no --override

No caso o programa faz interação com um outro chamado diamond entretanto o servidor que estamos trabalhando não comporta o uso deste software então o argumento 6 entra para inibir a interação e o último comando é para sobreescrever caso os arquivos ja tenham sido formados anteriormente e deseja-se refazer a análise as informações do resultado serão compiladas no arquivo "emapper.out.emapper.annotations" este arquivo e de dificil visualização por incompatibilidade das distancias entre tabulaturas então novamente usa-se o comando `column -t -s$'\t' <arquivo> | less -S` para ser extrair as informações mais importantes do arquivo foi criado o script twocolsdescmerge.pl ele funciona concaternado ao comando:

cut -f 1,8 emapper.out.emapper.annotations | sed 's/_i[0-9]\+\.p1//' | ./twocolsdescmerger.pl > gene_annotation.txt

O arquivo gerado será uma tabela separada por tabulação com a coluna um os identifcadores dos contigs de acordo com o arquivo Trinity.fasta e a coluna ao lado a descrição do gene de acordo com a base do EggNogMapper abaixo segue foto do arquivo.


Agora podemos juntar o resultado da anotação pelo programa com o arquivo CTR-TST.txt no diretorio output/trinity/abundance/DEG/ para isto sera utilizado novamente o script facilitador mergeR.R lembrando que caso o nome da coluna do arquivo "gene_annotation.txt" tiver "#" editar o arquivo retirando o simbolo pois a linguagem em R utilizado como argumento do script mergeR.R não aceita este simbolo.
Executar então a seguinte linha de comando:
`mergeR.R --x=<*abundance/DEG/CTR-TST.txt --by.x"x" --y=./gene_annotation.txt --by.y="query" --out=*abundance/DEG/CTR-TST.annot.txt --print.out.label`
Agora visualizar o arquivo pelo comando `column -t -s$'\t' <arquivo> | less -S`


# Gerando graficos do tipo Heatmap pelo RStudio com os resultados das análises
Neste tópico trataremos sobre as intruções iniciais para elaboração de gráfico do tipo heatmap para representar os dados de expressão obtidos pela simulação.
A linguagem R é muito utilizada para trabalhar com tabelas de dados, e tem diversos pacotes para instalação que tem funções específicas para análise de sequencias, (pacotes de R voltados para análise de biopolimetos, mais informações em: https://www.bioconductor.org/).

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:

`NOMEDESEJADO.df <- read.delim(file="caminho/do/arquivo/.txt", sep="\t", header=TRUE, stringsAsFactors = FALSE)`

Exibir somente as 2 primeiras linhas para checagem
head(NOME.df,2)

Dimensão do data.frame, onde o primeiro valor refere-se à quantidade de linhas e o segundo valor à quantidade de colunas:

dim(NOME.df)

selecionar um subconjunto de linhas do data.frame deg.df em um outro data.frame com log maior ou igual a 2 e menor ou igual a -2:
`("subset.NOME.df") logFC >= 2 ou logFC <= -2`

Além de possuírem um p-valor corrigido ("FDR") <= 0.05 concatenando "regras":

subset.NOME.df <- subset(NOME.df, (((logFC >= 2) | (logFC <= -2)  ) & (FDR <=0.05))  )

dim(subset.NOME.df)

head(subset.NOME.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.NOME.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.NOME.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), seq(-0.49,0.5,length=100), seq(0.51,1,length=100))

Criação de uma imagem de tamanho 5 x 5 polegadas

`png("./output/abundance/DEG/heatmap_DEGS_SAMPLEA-SAMPLEB.png",width = 5*300, height = 5*300, res = 300, pointsize = 8)
expression_data_to_plot, main = "Correlation distance (z-score)", density.info="none", trace="none", margins =c(12,12), col=my_palette, breaks=col_breaks, symbreaks=TRUE, dendrogram="both", distfun = function(x) as.dist(1-cor(t(x))), hclustfun = function(x) hclust(x, method="centroid"))

dev.off() # fecha o arquivo da imagem PNG

# Fluxograma de análises de Transcritomas