1° PARTE: MONTAGEM E ANÁLISE EXPRESSÃO GÊNICA

1° Passo:

Ao entrar na bioinfo, é possível verificar que o seu home é composto de duas pastas:

Vamo criar então o nosso local de trabalho com o comando mkdir:

        mkdir <seunome>

2° Passo:

        Vamos copiar todos os arquivos necessários para podermos realizar esta prática. Para isso, vamos para o seu local de trabalho:

        cd <seunome>

        Agora vamos copiar os arquivos da pasta data:

        cp -r ../data/* .

        

        Muito importante o ponto final depois do espaço do asterisco!!

        

Dentro de sua pasta é possível encontrar os seguintes arquivos:

3° Passo:

        Nesta aula utilizaremos dados de RNA-Seq correspondentes ao fungo Schizosaccharomyces pombe, contendo reads de 76 pares de base paired-end correspondentes a duas amostras: Sp_log (crescimento logarítmico) e Sp_plat (fase de platô). Dentro da pasta contendo essas reads, é possível encontrar os arquivos com final “left.fq” e “right.fq” no formato FASTQ formatado no estilo do sequenciador Illumina:

        Para poder realizar a montagem dos transcritos, vamos utilizar todas as reads para montar um único arquivo fasta contendo todos os transcritos montados. Isto é feito para futuramente utilizarmos este arquivo para realizar a análise de genes diferenciais. Iremos rodar o Trinity com o seguinte comando:

---------------------------------------------------------------------------------------------------------------------------

Trinity --seqType fq --SS_lib_type RF --left rnaseq_data/Sp_log.left.fq.gz,rnaseq_data/Sp_plat.left.fq.gz --right rnaseq_data/Sp_log.right.fq.gz,rnaseq_data/Sp_plat.right.fq.gz --CPU 1 --max_memory 1G --output trinity_saida/

---------------------------------------------------------------------------------------------------------------------------

Detalhes:

        Após completar a montagem, será criado uma pasta chamada trinity_saida, tal como o seu output foi descrito. Dentro desta pasta é possível encontrar todos os dados relacionados à montagem dos transcritos.

4° Passo:

        Vamos agora verificar o arquivo de montagem dos transcritos com o comando:

cd trinity_saida

less Trinity.fasta

Informações importantes:

CURIOSIDADE:

Como os transcritos são processados em forma de gráficos de De Bruijn, nós podemos visualizar utilizando o programa Bandage (https://rrwick.github.io/Bandage/).

Vídeo explicativo caso queira tentar em casa: https://www.youtube.com/watch?v=VuRN28XyFcI

5° Passo:

        Vamos analisar agora os dados relacionados a esta montagem utilizando um script do Trinity chamado TrinityStats.pl. Para isso, vamo para a pasta do trinity_saida:

        cd trinity_saida

Agora digite o comando:

TrinityStats.pl Trinity.fasta

Resultado:

################################

## Counts of transcripts, etc.

################################

Total trinity 'genes':  396

Total trinity transcripts:      400

Percent GC: 39.16

########################################

Stats based on ALL transcript contigs:

########################################

        Contig N10: 2589

        Contig N20: 2148

        Contig N30: 1818

        Contig N40: 1594

        Contig N50: 1329

        Median contig length: 548.5

        Average contig: 839.30

        Total assembled bases: 335720

#####################################################

## Stats based on ONLY LONGEST ISOFORM per 'GENE':

#####################################################

        Contig N10: 2589

        Contig N20: 2118

        Contig N30: 1818

        Contig N40: 1584

        Contig N50: 1322

        Median contig length: 545.5

        Average contig: 830.32

        Total assembled bases: 328805

Neste reporte é possível ter acesso à informações como: Número de genes encontrados, número de transcritos, a porcentagem de GC global e a média do tamanho dos contigs.

6° Passo:

        Vamos agora quantificar e contar a abundância de cada transcrito e gene que foram montados utilizando o Salmon para cada condição com os seguintes comandos:

        Volte agora para o seu diretório de trabalho:

        cd ..

E digite os comandos:

---------------------------------------------------------------------------------------------------------------------------

/home/bioufmg2/anaconda3/pkgs/trinityrnaseq-v2.12.0/util/align_and_estimate_abundance.pl --seqType fq --left rnaseq_data/Sp_log.left.fq.gz --right rnaseq_data/Sp_log.right.fq.gz --transcripts trinity_saida/Trinity.fasta --est_method salmon --trinity_mode --prep_reference --output_dir salmon_saida/sp_log/

---------------------------------------------------------------------------------------------------------------------------

/home/bioufmg2/anaconda3/pkgs/trinityrnaseq-v2.12.0/util/align_and_estimate_abundance.pl --seqType fq --left rnaseq_data/Sp_plat.left.fq.gz --right rnaseq_data/Sp_plat.right.fq.gz --transcripts trinity_saida/Trinity.fasta --est_method salmon --trinity_mode --prep_reference --output_dir salmon_saida/sp_plat/

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

        Ao final do processo, irá ser criado uma pasta chamada salmon_saida com duas pastas: Sp_log e Sp_plat, contendo todos os resultados obtidos pelo Salmon:

7° Passo:

        

Vamos para a pasta contendo os resultados de uma das condições para podermos analisar os arquivos de saída gerados pelo Salmon:

Dentro da pasta podemos encontrar 2 arquivos importantes, o quant.sf e o quant.sf.genes

        Estes arquivos são separados por tabulações e possuem colunas das quais podemos identificar o nome do transcrito/gene, o tamanho, e o valores de TPM e de suas reads:

        Utilizando o comando:

less quant.sf.genes

Obtemos o seguinte resultado (não se preocupe se o seu resultado estiver diferente da imagem, o importante são as colunas e seus valores):

        Fique a vontade agora para explorar os outros arquivos.

8° Passo:

        Agora vamos comparar os genes e transcritos entre as amostras utilizando um script do pacote Trinity chamado: abundance_estimates_to_matrix.pl. Primeiramente vamos para a nossa pasta  de trabalho:

                

        cd ../../

Este script irá calcular tanto a nível de genes quanto de transcritos com o seguinte comando:

---------------------------------------------------------------------------------------------------------------------------  /home/bioufmg2/anaconda3/pkgs/trinityrnaseq-v2.12.0/util/abundance_estimates_to_matrix.pl --est_method salmon --out_prefix salmon_saida/salmon salmon_saida/sp_log/quant.sf salmon_saida/sp_plat/quant.sf --gene_trans_map trinity_saida/Trinity.fasta.gene_trans_map --name_sample_by_basedir

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

Agora com ls na pasta do salmon_saida podemos ver que apareceram novos arquivos:

        Arquivos importantes:

Vamos para a salmon_saida:

        cd salmon_saida

E digitar:

 less salmon.gene.counts.matrix

Podemos observar que as colunas dos valores foram nomeadas de acordo com a pasta que elas estavam, facilitando o trabalho, nos próximos passos, para identificar qual dado é de qual amostra é condição. Chegaremos lá!

Também fique a vontade para ver os outros arquivos e perceber que todos eles seguem o mesmo padrão.

9° Passo:

        Partiu agora descobrir os genes que estão diferencialmente expressos e montar uns gráficos para melhor visualização destes dados que criamos. Vamos utilizar o programa edgeR para este cálculo e outro script do Trinity chamado run_DE_analysis.pl. 

        Vamos agora sair da pasta do salmon_saida e vamos para a nossa pasta de trabalho:

        cd ..

Certifique-se que está na sua pasta de trabalho para evitar problemas. Caso esteja perdido, utilize o comando pwd para ver o local do seu terminal.

Bora rodar:

---------------------------------------------------------------------------------------------------------------------------

/home/bioufmg2/anaconda3/pkgs/trinityrnaseq-v2.12.0/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix salmon_saida/salmon.gene.counts.matrix --method edgeR --output edgeR_saida/ --samples_file samples.txt --dispersion 0.1

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

O arquivo samples.txt é muito importante nessa etapa uma vez que ele indica qual a condição e suas replicatas para o programa. Ele possui somente duas colunas, com suas linhas separadas por tabulações, e os nomes destas colunas devem estar de acordo com o nomes do passo 7, dentro do arquivo que possui o final .counts.matrix. Exemplo:

#condition        #samples

condA                    repA1

condA                repA2

condB                repB1

condB                repB2

10° Passo:

        Com os cálculos do edgeR realizados, vamos agora ver os arquivos de saída. Vamos para a pasta do edgeR_saida:

        cd edgeR_saida

        Arquivos importantes:

Para ficar facilitar a visualização no navegador, digite o comando:

---------------------------------------------------------------------------------------------------------------------------

cp salmon.gene.counts.matrix.Log_Phase_vs_Plat_Phase.edgeR.DE_results.MA_n_Volcano.pdf ../MA_Volcano.pdf

---------------------------------------------------------------------------------------------------------------------------

        Abra agora uma nova página no seu navegador no link:

http://bioinfo.icb.ufmg.br/bioufmg/look/

        Será necessário um login (bioufmg) e senha (carrarataxis). Ao acessar, entre na pasta look e depois a pasta com o seu nome.

Ao abrir o arquivo PDF (MA_Volcano.pdf) podemos verificar que as características com um valor de FDR < 0.5 estão coloridas de vermelho.

11° Passo:

        Quantos genes estão diferencialmente expressos? Vamos descobrir com o seguinte comando (dentro da pasta edgeR_saida! pwd se estiver perdido):

---------------------------------------------------------------------------------------------------------------------------

sed '1,1d' salmon.gene.counts.matrix.Log_Phase_vs_Plat_Phase.edgeR.DE_results | awk '{ if ($7 <= 0.05) print;}' | wc -l

---------------------------------------------------------------------------------------------------------------------------

Este comando irá quantificar quantos genes possuem um valor de FDR < 0.05 (coluna 7). Sinta-se livre para brincar com este valor.

12° Passo:

        

        Vamos agora pegar todas as sequências de todos os genes que possuem o FDR <= 0.05 utilizando o script filtrar.sh. Para isso, vamos voltar para a nossa pasta de trabalho:

        cd ../

        Agora só rodar:

        sh filtrar.sh

        Após dar um ls, é possível ver que foram criados 2 arquivos:

        DE_genes.txt: contém o nome de cada gene;

        DE_genes.fasta: contém o header e a sequência do gene;

---------------------------------------------------------------------------------------------------------------------------

2° PARTE: ANOTAÇÃO  E EXPLORAÇÃO DOS GENES DIFERENCIALMENTE EXPRESSOS

1° Passo:

        Iremos agora extrair esses genes e gerar os famosos gráficos de heatmap. Para isso utilizaremos o analyze_diff_expr.pl do Trinity, mas precisamos estar dentro da pasta edgeR_saida:

cd edgeR_saida/

---------------------------------------------------------------------------------------------------------------------------

/home/bioufmg2/anaconda3/opt/trinity-2.1.1/Analysis/DifferentialExpression/analyze_diff_expr.pl  -P 1e-3 -C 2 --matrix ../salmon_saida/salmon.gene.TMM.EXPR.matrix --samples ../samples.txt

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

Dando um ls nós vemos vários arquivos novos que apareceram:

        Arquivos importantes:

Vamos agora brincar com estes arquivos.

2° Passo:

        Vamos conhecer os genes super regulados em uma amostra com um:

---------------------------------------------------------------------------------------------------------------------------

less salmon.gene.counts.matrix.Log_Phase_vs_Plat_Phase.edgeR.DE_results.P1e-3_C2.Log_Phase-UP.subset

---------------------------------------------------------------------------------------------------------------------------

Lembre-se de estar dentro da pasta edgeR_saida. Podemos perceber que este arquivo possui a separação em suas linhas por tabulações, facilitando a manipulação destes arquivos. Fique a vontade para conhecer os genes da outra amostra!

3° Passo:

        

        Com um  less diffExpr.P1e-3_C2.matrix.log2.centered.dat nós podemos perceber que este arquivo também é tabulado:

4° Passo:

        Os resultados dos seus heatmaps gerados devem se parecer com os mostrados abaixo.

        Para facilitar visualização no navegador, digite os comandos dentro da pasta edgeR_saida:

---------------------------------------------------------------------------------------------------------------------------

cp diffExpr.P1e-3_C2.matrix.log2.centered.genes_vs_samples_heatmap.pdf ../heatmap_cluster.pdf

---------------------------------------------------------------------------------------------------------------------------

cp diffExpr.P1e-3_C2.matrix.log2.sample_cor_matrix.pdf ../heatmap_cor_matrix.pdf

---------------------------------------------------------------------------------------------------------------------------

5° Passo:

        Agora vamos terminar as análises de genes diferenciados e fazer uma clusterização dos mesmos para melhor visualização de grupos que com um padrão similar nas amostras!

        Vamos dar o seguinte comando, ainda dentro da pasta edgeR_saida:

---------------------------------------------------------------------------------------------------------------------------

../../anaconda3/opt/trinity-2.1.1/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl -R diffExpr.P1e-3_C2.matrix.RData --Ktree 3

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

        Foi criado agora uma nova pasta, vamos para ela.

6° Passo:

        Dentro da pasta edgeR_saida vamos dar um:

        cd diffExpr.P1e-3_C2.matrix.RData.clusters_fixed_Ktree_3/

Lembrando que o nome da pasta depende dos valores de p-valor, FC e Ktree utilizados nos passos do script analyze_diff_expr.pl e define_clusters_by_cutting_tree.pl!

7° Passo:

        Observando os arquivos desta pasta podemos verificar que foi adicionado um novo arquivo pdf e outros de texto.

        Para facilitar, vamos dar o seguinte comando, dentro da pasta diffExpr.P1e-3_C2.matrix.RData.clusters_fixed_Ktree_3/:

---------------------------------------------------------------------------------------------------------------------------

cp my_cluster_plots.pdf ../../my_cluster_plots.pdf

---------------------------------------------------------------------------------------------------------------------------

O seu PDF deve se parecer com algo assim:

8° Passo:

        Para poder saber exatamente quem são esses grupos, basta somente verificar cada arquivo que foi gerado a partir do algoritmo de clusterização com o comando que preferir!

Para cada cluster é gerado este tipo de arquivo tabulado. Caso tenha curiosidade, você pode explorar os outros sub-clusters criados.

Por fim, volte a sua pasta de trabalho com:

        cd ../../

9° Passo:

        Agora vamos utilizar o Diamond para poder identificar a nomenclatura de cada gene para podermos conhecê-los. Para isso, certifique-se que esteja na sua pasta de trabalho. Caso esteja perdido: pwd

        Para rodar o Diamond, rode o seguinte comando:

---------------------------------------------------------------------------------------------------------------------------

Diamond blastx -d ../data/uniprot_sprot.dmnd -q DE_genes.fasta -v --outfmt 6 -e 1e-10 --max-target-seqs 1 --outfmt 6 -o blastx.outfmt6 -b4

---------------------------------------------------------------------------------------------------------------------------

        Detalhes:

10° Passo:

        Vamos agora conhecer os nossos transcritos a partir do arquivo blastx.outfmt6. Se dermos um less neste arquivo podemos perceber que é um arquivo texto tabulado:

Nesse caso temos as seguintes colunas:

11° Passo:

        Vamos agora identificar quais genes essas proteínas pertencem utilizando a plataforma online do Uniprot na aba de Get Gene Name. Primeiramente vamos printar somente a segunda coluna do nosso blast, referente aos nomes das proteínas:

awk '{print $2}' blastx.outfmt6

        Selecione todas as linhas contendo estes nomes:

Lembrando que o PuTTy copia automaticamente ao selecionar alguma coisa no terminal.

12° Passo:

        Acesse agora, numa nova aba do seu navegador, o site uniprot.org.

        

        Logo no cabeçalho do site vamos clicar em Retrieve/ID mapping:

        Cole o texto que você copiou do PuTTY (ctrl + v) na região do Provide your identifiers; selecione, na parte de Select Options, To Gene Name. Por fim, clique em Submit.

Irá abrir então uma página contendo todos os gene names dos seus identificadores:

13° Passo:

        Retorne à página anterior do seu browser para fazermos uma nova busca utilizando os mesmos identificadores. Desta vez, vamos selecionar To KEGG dentro do grupo de Genome:

        

        Agora irá abrir uma página contendo os identificadores KEGG que dão acesso diretamente a esse banco de dados:

14° Passo:

        Vamos explorar alguma dessas vias. Clique na primeira opção de KEGG com o botão direito e vá em Abrir em uma nova guia:

        Na aba que irá abrir, será do banco de dados do KEGG:

        A partir daqui nós podemos visualizar todas as vias metabólicas que este gene está envolvido, como por exemplo, a spo00010 Glycolysis/Gluconeogenesis. Vamos fazer o mesmo procedimento poder abrir este link em uma nova aba do navegador:

        Nesta nova aba aberta agora podemos explorar em qual região da via este gene está com sua expressão diferenciada. Neste caso, está em vermelho o número 3.1.3.11 dentro da parte de Pentose phophate pathway.

        Agora você pode explorar todas as vias e genes que você identificou como diferencialmente expressos. Boa sorte.