1° Passo:

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

• anaconda3: pasta contendo todo o ambiente dos programas que serão utilizados nesta aula;
• data: pasta com todos os arquivos de entradas necessários para realizarmos a montagem dos transcritos.

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:

• rnaseq_data: pasta com todas as reads que serão utilizadas;
• uniprot_sprot.pep: arquivo contendo os fastas de proteínas do banco de dados SwissProt - Uniprot;
• uniprot_sprot.dmnd: arquivo de banco de dados criado pelo Diamond a partir dos fastas do Uniprot;
• diamond: o alinhador que iremos utilizar para anotar os transcritos montados;
• samples.txt: arquivo que descreve as amostras e as condições utilizadas pelo experimento;
• filtrar.sh: um scriptzinho para poder pegar as sequencias dos genes lá no final.

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:

• --left (rnaseq_data/Sp_log.left.fq.gz,rnaseq_data/Sp_plat.left.fq.gz): indicamos o caminho de todas as reads no sentido forward;
• --right (rnaseq_data/Sp_log.right.fq.gz,rnaseq_data/Sp_plat.right.fq.gz): indicamos o caminho de todas as reads no sentido reverse;
• Não se esqueça de que estes arquivos devem ser separados por vírgulas, e não espaçamentos.
• --SS_lib_type (RF): seleciona o tipo da sua biblioteca, se é paired-end ou single-end;
• --CPU (1): quantidade de threads que iremos utilizar;
• --max_memory (1G): quantidade máxima de memória que o Trinity irá utilizar;
• --output (trinity_saida/): local no qual o Trinity irá salvar todos os dados relacionados à montagem dos transcritos.
• Importante que esta pasta esteja dentro do seu diretório de trabalho e nomeada como trinity_saida.

        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:

• O identificador “>” é utilizado como todo arquivo fasta, ser um cabeçalho para a sequência;
• TRINITY_DN94, como no exemplo, é o nome do transcrito montado;
• _g1: indica que é o gene 1;
• _i1: indica que é a isoforma 1.

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.

• Nx: Todos os contigs são ordenados, a partir do seu tamanho, do maior para o menor; o menor contig montado dentre os maiores compreendendo x% dos nucleotídeos totais, representa o seu Nx.

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:

• align_and_estimate_abundance.pl: script do pacote Trinity que ira realizar todo o pipeline;
• --seqType (fq): indicamos qual o formato dos arquivos das reads;
• --left (rnaseq_data/Sp_log.left.fq.gz): indicamos o caminho de todas as reads no sentido forward;
• --right (rnaseq_data/Sp_log.right.fq.gz): indicamos o caminho de todas as reads no sentido reverse;
• --transcripts (trinity_saida/Trinity.fasta): o caminho do arquivo contendo todos os transcritos montados pelo Trinity;
• --est_method (salmon): indicamos qual o método/software de estimação iremos utilizar;
• --trinity_mode: irá gerar um arquivo que irá mapear todas as isoformas em seus respectivos genes;
• --prep_reference: irá criar índices no arquivo dos transcritos para agilizar os processos;
• --output_dir (salmon_saida/sp_log/): indicamos qual o local dos arquivos de saída.

        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:

• cd salmon_saida/sp_log

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

• quant.sf: quantificação da abundância dos transcritos;
• quant.sf.genes: quantificação de abundância a nível de 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:

• --est_method (salmon): indicamos qual o método/software de estimação foi utilizado para o cálculo de abundância;
• --out_prefix (salmon_saida/salmon): indicamos a parta e o nome dos arquivos de saída;
• --gene_trans_map (trinity_saida/Trinity.fasta.gene_trans_map): aquele índice contento cada transcrito pareado com a sua isoforma gerado no passo 5;
• --name_sample_by_basedir: utilizado para poder nomear cada condição de acordo com o nome da pasta que as contém.

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

        Arquivos importantes:

• Aqueles que possuem isoform são os resultados a nível dos transcritos e aqueles que possuem gene são a nível de genes;
• Arquivos terminados em .counts.matrix possuem a estimativa da contagem dos fragmentos de RNA-Seq;
• .TPM.not_cross_norm: matriz contendo os valores de expressão TPM que não estão normalizadas por amostra;
• .TMM.EXPR.matrix: matriz contendo os valores de expressão TMM normalizados.

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:

• --matrix (salmon_saida/salmon.gene.counts.matrix): indicamos de qual matriz contém os dados de abundâncias dos transcritos de genes. Neste caso, iremos calcular somente os genes diferenciados e não suas isoformas.
• --method (edgeR): indicamos qual programa usar;
• --samples_file (samples.txt): um arquivo que descreve as replicatas e condições dos seus experimentos. Já já iremos vê-lo;
• --dispersion (0.1): termo utilizado na função binomial negativa do edgeR para estimar a contagem dos genes.

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:

• O nome do arquivo segue uma lógica: {prefix}.sampleA_vs_sampleB.{method};
• .DE.results: arquivo contendo os resultados das análises, incluindo o fold change e a significância estatística (FDR);
• .MA_n_Volcano.pdf: nossos primeiros gráficos, sendo um de MA e um Volcano!

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/

        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;

13° 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:

• Estamos rodando o blastx do Diamond;
• -d (../data/uniprot_sprot.dmnd): indicamos o caminho do banco de dados que será utilizado para blastar as sequências. Neste caso, estamos utilizando o SwissProt-Uniprot;
• -q (DE_genes.fasta): indicamos a query, que neste caso é o arquivo que contém a sequência dos genes que possuem um valor de FDR <= 0.05;
• -v: modo verbose, para poder mostrar letrinhas no terminal;
• -e (1e-10): valor de e-value;
• --max-target-seqs (1): para poder anotar somente o best-hit de cada transcrito;
• --outfmt (6): define o modelo do qual será o arquivo de saída, que neste caso, é igual ao formato de saída do programa Blast;
• -o (blastx.outfmt6): arquivo de output.
• -b4: argumento para poder não limitar a memória no uso do Diamond.

14° 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:

• Coluna 1: Nome da query;
• Coluna 2: Nome da referência utilizada;
• Coluna 3: porcentagem de identidade;
• Coluna 4: tamanho do alinhamento;
• Coluna 5: número de mismatches;
• Coluna 6: número de gaps;
• Coluna 7; posição de início na sequência da query;
• Coluna 8: posição final na sequência da query;
• Coluna 9: posição de início na sequência da referência;
• Coluna 10: posição final na sequência da referência;
• Coluna 11: e-value;
• Coluna 12: bit score.

15° 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.

16° 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:

17° 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:

18° 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.

PASSOS ADICIONAIS OPTATIVOS:

19° 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:

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

/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:

• -P (1e-3): define o p-valor a ser utilizado para ser considerado um gene diferencialmente expresso;
• -C (2): define o valor do fold change;
• --matrix (../salmon_saida/salmon.gene.TMM.EXPR.matrix): indicamos de qual matriz contém os dados de abundâncias dos transcritos de genes gerados no passo 7 contendo os valores de TMM entre as amostras;
• --samples (../samples.txt): aquele mesmo arquivo que descreve as replicatas e condições do seus experimentos.

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

        Arquivos importantes:

• .{sampleA}-UP.subset: features que estão reguladas positiviamentes na amostra A;
• .{sampleB}-UP.subset: features que estão reguladas positiviamentes na amostra B;
• diffExpr.{pvalor}_{valorFC}.matrix.log2.dat: todas a features encontradas que estão diferencialmente expressas em todas as comparações entre as amostras;
• diffExpr.{pvalor}_{valorFC}.matrix.log2.sample_cor.dat: uma matriz de correlação de Pearson para comparações entre as amostras baseada nas features diferencialmente expressas encontradas;
• diffExpr.{pvalor}_{valorFC}.matrix.log2.centered.genes_vs_samples_heatmap.pdf: nome muito grande para poder falar que é um gráfico de heatmap clusterizado dos genes diferencialmente expressos em cada replicata de suas amostras;
• diffExpr.{pvalor}_{valorFC}.matrix.log2.sample_cor_matrix.pdf: outro heatmap clusterizado mostrando as matrizes de correlações entre as amostras;

Vamos agora brincar com estes arquivos.

20° 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!

21° Passo:

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

22° 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

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

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

23° 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:

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

define_clusters_by_cutting_tree.pl -R diffExpr.P1e-3_C2.matrix.RData --Ktree 3

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

        Detalhes:

• -Ktree: Indica que vamos utilizar valor de k-means para poder fazer esta clusterização. É possível brincar com este valor até achar um que melhor represente a sua amostra!

        Foi criado agora uma nova pasta, vamos para ela.

24° 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!

25° 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:

26° 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.