HackMDBioinformática II: Análise de Transcriptomas
Atualizado em Nov/2021
URL curta: http://bit.ly/bioinfotranscriptomics
Sumário
Simulação
A simulação a seguir, necessita do pacote E-Direct com as ferramentas E-utilities. Para maiores informações, consultar a documentação Entrez Direct E-utilities on the UNIX Command Line.
Requisitos
Seguir as instruções para instalação do e-direct.
- Caso precise atualizar o conjunto de programas do e-direct, remova primeiro a pasta correspontende
- Copiar e colar o código a seguir
- Adicionar o caminho no $PATH caso não tenha solicitado que seja feito automaticamente no final do passo anterior
- Para testar:
Construção do transcriptoma para simulação
Teremos que utilizar um transcriptoma com um genoma de referência. O exemplo a seguir utilizará genes de Trichoderma harzianum (Agente modelo para uso como biocontrole de fitopatógenos).
Obtenha o Taxonomy ID para esta espécie: 5544
Dessa forma, para as consultads nos bancos de dados do NCBI, utilizar: txid5544[Organism:exp]
A seguir vamos utilizar as ferramentas esearch/efetch para fazer o download e construir o arquivo fasta contendo os transcritos referência para a simulação dos experimentos de investigação do transcriptoma.
Para conhecer os IDs de interesse, consultar a página do banco de dados Gene do NCBI.
e depois de selecionar o gene de interesse:
Identificar o número de acesso do transcrito. No caso de transcritos no banco de dados de referência RefSeq, os IDs iniciam com o prefixo "NM_" ou "NR_", respectivamente para transcritos codificadores de proteínas (mRNAs) e transcritos não codificadores de proteínas (ncRNAS). Os IDs também podem ser respectivamente "XM_" ou "XR_", caso essas entradas tenham somente uma predição computacional, sem ainda terem passado pela etapa de curadoria.
Procurar pelos números de acesso dos transcritos. No caso do gene acima, há 2 transcritos (2 isoformas).
Números de acesso do gene M431DRAFT_343923:
XM_024914431.1 → XP_024777235.1 hypothetical protein M431DRAFT_343923 [Trichoderma harzianum CBS 226.95]
XM_024914432.1 → XP_024777236.1 hypothetical protein M431DRAFT_343923 [Trichoderma harzianum CBS 226.95]
Selecionar mais transcritos de outros genes:
Com organismos procariotos é necessário copiar e colar manualmente o resultado da busca do NCBI no banco de dados gene com as coordenadas do gene referentes ao transcrito.
Veja exemplo do gene COR42_RS23510 que codifica proteína DUF3018 family protein de Xanthomonas citri pv. citri strain Xcc29-1 plasmid pXAC64.
O número de início (start) e fim (stop) têm que ser ajustados. No caso acima, COR42_RS23510 possui as seguintes coordenadas 31468..31686 complement, dessa forma, devemos considerar 1 base a menos no início e uma base a menos no fim, ou seja, neste caso 31467 e 31685, respectivamente para início e fim, assim como o fato de estar na fita complementar (minus).
E no caso do gene COR42_RS23335 que codifica proteína recombinase family protein de Xanthomonas citri pv. citri strain Xcc29-1 plasmid pXAC33 (NCBI Genome).
No caso de COR42_RS23335, o mesmo ajuste deve ser feito.
A fim de simplificar, podemos fazer o mesmo processo utilizando uma estrutura de repetição.
Para facilitar ainda mais, podemos criar um arquivo com os números de acesso na primeira coluna e utilizá-lo com o código a seguir.
A seguir o arquivo construído com os IDs dos transcritos e o nome dos genes na segunda coluna (As colunas devem ser separadas por <TAB> e não espaços simples). Isso auxiliará no futuro a identificar quais transcritos relacionam-se com quais genes.
Arquivo ACCS.txt:
XM_024918093.1 M431DRAFT_502101 XM_024914151.1 M431DRAFT_300419 XM_024918802.1 M431DRAFT_506422 XM_024922757.1 M431DRAFT_78900 XM_024914431.1 M431DRAFT_343923 XM_024914432.1 M431DRAFT_343923
Para procariotos:
Arquivo ACCS.txt:
COR42_RS23510 NZ_CP024033.1 31467 31685 minus COR42_RS23335 NZ_CP024032.1 30895 31497 plus
Se executar o comando a seguir e encontrar uma linha para cada transcrito, o processo foi bem sucedido.
Criando os arquivos de abundância
Os arquivos abundance_A.txt e abundance_B.txt devem ser criados manualmente, respectivamente para a condição biológica A e condição biológica B. Os mesmos números de acesso usados na construção do arquivo transcriptoma.fa devem ser utilizados (coluna 1) juntamente com a proporção de cada transcrito (coluna 2) considerando que pode haver diferenças de abundância entre essas duas condições biológicas.
ATENÇÃO: As duas colunas devem ser separadas por TAB, sem linhas adicionais e sem espaços.
- Para eucariotos
abundance_A.txt
XM_024918093.1 0.50 XM_024914151.1 0.05 XM_024918802.1 0.05 XM_024922757.1 0.25 XM_024914431.1 0.15 XM_024914432.1 0
abundance_B.txt
XM_024918093.1 0.05 XM_024914151.1 0.05 XM_024918802.1 0.05 XM_024922757.1 0.25 XM_024914431.1 0.10 XM_024914432.1 0.50
- Para procariotos
abundance_A.txt
COR42_RS23510 0.8 COR42_RS23335 0.2
abundance_B.txt
COR42_RS23510 0.2 COR42_RS23335 0.8
Execute os comandos abaixo, o resultado da soma dos valores da coluna 2 tem que ser 1 para ambos os arquivos:
Genoma reduzido
Caso o genoma seja muito grande, é possível criar uma versão reduzida do genoma somente com os cromossomos dos transcritos selecionados. Para as sequências em fasta:
E para as coordenadas no GFF:
Gerando os arquivos .fastq
Para exemplificar, abaixo estão apenas os comandos para gerar a primeira réplica da simulação da corrida para a amostra A utilizando o (abundance_A.txt).
- Gerar "25000" sequências de fragmentos aleatórios em média de "300 pb" com desvio padrão de "30 pb" a partir das sequências do transcriptoma montado anteriormente "transcriptoma.fa", considerando a abundância contida em "abundance_A.txt"
- Renomear as sequências para "SAMPLEA" juntamente com um número hexadecimal sequencial começando por "0000000001" e posicionar todas as bases em uma única linha (até 1000 pb). A linha de SED irá substituir a descrição adicional.
- Fazer a simulação considerando leituras "paired" 2x "151" bp , com parâmetros configurados ([https://www.ebi.ac.uk/goldman-srv/simNGS/runfiles5/151cycPE/s_4_0099.runfile]) usando simNGS:
- Desintercalar as leituras em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq no diretório ../raw
- Remover arquivos desnecessários:
- Verificar o número de sequências geradas
- Repetir o mesmo processo para a réplica 2 da condição biológica A (SAMPLEA2) utilizando o arquivo abundance_A.txt
…
- Repetir o mesmo processo para a réplica 1 e réplica 2 para a condição biológica B (SAMPLEB).
…
Impacientes
Os impacientes que já tiverem os arquivos ACCS.txt, abundance_A.txt e abundance_B.txt e podem acompanhar a partir deste ponto. Os impacientes que não tiverem os arquivos, manter a calma e ler com atenção as seções anteriores.
As etapas de simulação realizadas anteriormente podem ser feitas utilizando os scripts sim.sh para eucariotos ou sim-prok.sh para procariotos.
Obtenção das referências para análises
Sequências e Anotações
Sequências no formato FASTA e anotações gênicas no formato GFF.
Página do genoma para Trichoderma harzianum. Utilizar a consulta pelo Taxonomy ID ().
- Utilize gunzip para descompactar
Como fazer o download de corridas do SRA:
ATENÇÃO: NÃO É NECESSÁRIO FAZER O DOWNLOAD DOS DADOS DO SRA CASO ESTEJA PARTINDO DOS DADOS PROVENIENTES DA ETAPA ANTERIOR DE SIMULAÇÃO.
ESTA SEÇÃO É PARA OS QUE DESEJAM APRENDER UMA FORMA DE OBTER OS DADOS DE EXPERIMENTOS REAIS E SEGUIR O MESMO PROTOCOLO. OBVIAMENTE SERÃO NECESSÁRIAS ADAPTAÇÕES NOS NOMES DAS CONDIÇÕES BIOLÓGICAS E DOS ARQUIVOS.
APÓS ESTA SEÇÃO, É POSSÍVEL SEGUIR O TUTORIAL COM OS DADOS DA SIMULAÇÃO.
ATENÇÃO: NÃO FAZER O DOWNLOAD DOS ARQUIVOS DO NCBI SRA NO SERVIDOR HAMMER!
Os comandos a seguir podem ser executados em um outro computador com mais espaço.
- Criar diretório apropriado
- Consultar por organismo Trichoderma harzianum:
txid5544[Organism:noexp]
- Identificar o ID da corrida (prefixo SRR) a partir do experimento (prefixo SRX):
- Utilizar o comando prefetch do NCBI Toolkit:
- O arquivo .sra será baixado para o diretório ~/ncbi/public/sra/
- Obter os arquivos .fastq referentes a esta corrida. O parâmetro –split-3 indica que se deseja separar as sequências em arquivos _1.fastq e _2.fastq. O parâmetro –origfmt indica que se deseja utilizar o nome original das reads.
Read 22339605 spots for SRR8707026
Written 22339605 spots for SRR8707026
Limpeza dos arquivos referência
Voltar para diretório ./refs:
Limpeza dos cabeçalhos FASTA
Criar arquivo com o script (cleanfasta.sh) abaixo:
Ctrl+c & Ctrl+v
Dar permissão:
Executar:
Ajustes do arquivo GFF:
Executar o script fixNCBIgff.sh do repositório Github:
O script fixNCBIgff.sh possui dependências:
- fixGFFgenestruct.pl
- fixGFFdupID.pl
- sortGFF.pl
Veja descrição das colunas do arquivo GFF e GFF2/GTF.
Testes com GFF
- Para a conversão de GFF versão 3 para um arquivo GTF:
- Para obtenção de um arquivo com a sequência FASTA do trasncriptoma:
- Para obtenção de um arquivo com a sequência FASTA do proteoma:
Testes de alinhamento com Bowtie2
NÃO É NECESSÁRIO REALIZAR ALINHAMENTOS COM BOWTIE2 PARA AS ANÁLISES USANDO O GENOMA REFERÊNCIA. É APENAS TESTE.
Baixar o Integrative Genomics Viewer (IGV) para visualização dos alinhamentos.
- Indexação do genoma
- Alinhamento
- Conversão de SAM para BAM
Devem ser copiados os seguintes arquivos para serem abertos no IGV:
bowtie2-test-sorted.bam
bowtie2-test-sorted.bam.bai
Mapeamento de RNA-Seq reads x Genoma:
Observar que há falta de cobertura nas regiões correspondentes a introns.
Testes de alinhamento com TopHat
ATENÇÃO: O TopHat tem problemas de compatibilidade com versões mais recentes do bowtie2. Verifique a compatibilidade antes de utilizá-los.
No caso do TopHat basta indexar o arquivo .bam resultante para incorporá-lo ao IGV.
Notar o alinhamento das leituras que atravessam os introns e alinham concordantemente com a estrutura gênica contida no arquivo GTF.
Testes de alinhamento com STAR
- Criação do índice
- Entrar no diretório
- Criar diretório
- Criação do índice
Alinhamento:
O resultado no formato .sam deve ser convertido para .bam para que possa ser utilizado no IGV.
Procura por: XM_024918093.1
Notar o alinhamento das leituras que atravessam os introns e alinham concordantemente com a estrutura gênica contida no arquivo GTF.
Testes com fastqc (pre)
Testes com prinseq-lite.pl
Poda de qualidade baseada em janela deslizante.
AATAGATACGATGAGAGACCAGTAGACGATCA
||||||
|||||10
||||20
|||25
||33
|30
30
mean(c(25, 20, 10)) = 18.33
mean(c(33, 25, 20)) = 26
mean(c(30, 33, 25)) = 29.33
mean(c(30, 30, 33)) = 31
-trim_qual_window 3
-trim_qual_step 1
-trim_qual_right 30
-trim_qual_rule lt
-trim_qual_type mean
AATAGATACGATGAGAGACCAGTAGACGA
Testes com fastqc (pos)
Criando script com loop
Utilizar nano ou vim para criar o arquivo ./scriptfor.sh com conteúdo abaixo:
Salvar e sair do editor. Depois, dar permissão de execução e, por fim, executar:
Script FastQC(pre)+PrinSeq+FastQC(pos)
O script preprocess1.sh foi desenvolvido para uma análise simples considerando que há arquivos *_R1.fastq e *_R2.fastq no diretório ./raw. O script deve ser gravado no diretório que contém o subdiretório ./raw :
… e executado da seguinte forma:
O script preprocess2.sh tem a mesma funcionalidade que o anterior, porém é muito mais generalista, pois permite a execução a partir de qualquer local, desde que sejam informados o diretório de entrada contendo os arquivos *_R1.fastq e *_R2.fastq (1º argumento) e o diretório de saída (2º argumento), onde serão gravados os diretórios referentes às análises com FastQC e PrinSeq.
]
Testes com atropos (adapter trimming)
Para chamar o comando atropos no servidor hammer é necessário ativar o ambiente Python compatível com o programa.
Com o ambiente carregado é possível invocar o comando e visualizar as opções disponíveis para a poda (trimagem).
Para desativar o ambiente e retornar a o ambiente padrão:
ATENÇÃO: Não desativar o comando até terminar a execução do atropos
A primeira execução será no modo insert:
Consulte o manual para descrever o uso de cada um desses parâmetros usados na linha de comando. Também é possível consultar por meio da chamada ao help:
A saída do atropos no modo insert das reads que não foram podadas serão as entradas do atropos no modo adapter:
Concatenando os resultados de ambos os modos de execução:
Para depois remover os arquivos já concatenados no arquivo com sufixo .atropos_final.fastq.
Os adaptadores utilizados nesses comandos referem-se à seguinte construção:
E às sequências:
-g AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
-G GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
-a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
Montando o script de pré-processamento final
A ordem de execução será a seguinte:
FastQC (pré) -> ATROPOS (insert) -> ATROPOS (adater) -> resultados concatanados ATROPOS (insert) & ATROPOS (adapter) -> PrinSeq -> FastQC (pós)
Dessa forma, o script preprocess2.sh foi atualizado para incorporar as etapas de pré-processamento com ATROPOS no script preprocess3.sh.
Montando o script de análise RNA-Seq com referência
O script rnaseq-ref.sh foi criado utilizando o script preprocess3.sh para fazer o pré-processamento e posteriormente executar o STAR com as leituras Paired-End e as leituras Single-End (singletons), além de combinar os resultados, ordenar e executar o script SAM_nameSorted_to_uniq_count_stats.pl para coletar as estatísticas de alinhamento.
Para verificar que o script foi executado com sucesso, confira com a listagem dos arquivos resultantes dos alinhamentos e que foram posteriormente ordenados com samtools.
Verifique a porcentagem geral de alinhamento em pares de modo adequado (proper pairs):
Execução do pipeline
O script rnaseq-ref.sh deve ser executado da seguinte forma:
Caso tenha problemas de execução, fazer as seguintes conferências:
CHECKLIST do sucesso:
:smiley:
- Os seguintes diretórios e arquivos de entrada existem: ./raw, ./output, ./refs/genome.gtf e ./refs/genome.fa. No caso dos arquivos (genome.gtf e genome.fa), eles devem possuir conteúdo, representado na linha abaixo em número de Bytes, KBytes ou MegaBytes:
- Os arquivos contidos em ./raw possuem nomenclatura adequada e devem estar em pares R1 e R2:
NOMEDAAMOSTRA_R[12].fastq
Exemplo:
SAMPLEB1_R1.fastq e SAMPLEB2_R2.fastq
- Todos os arquivos ./raw possuem conteúdo, seja por número de Bytes, número de linhas ou quantidade de sequências, respectivamente:
… usando um outro código :neutral_face: shell para contar sequências em arquivo .fastq:
Avaliando a presença de contaminantes
Supondo que possa haver contaminantes na amostra, fazemos um mapeamento das leituras com sequências do contaminante.
No diretório ./refs obter e formatar a(s) sequência(s) do contaminante. Neste caso Bacillus subtilis:
Descompactar:
Construir o índice:
Alinhar:
Redução do genoma
Caso estejam limitados com relação à memória e processamento do servidor e tenham que reduzir o tamanho do genoma, devem fazer isso com o arquivo FASTA das sequências genômicas e também o arquivo GTF/GFF.
O primeiro passo é criar um arquivo contendo uma lista com os nomes das sequências desejadas.
Vamos selecionar somente os cromossomos referentes aos transcritos selecionados na simulação. Para começar a filtrar as sequências dos cromossomos, vamos montar uma lista não redundante com os números de acesso das sequências de RNA.
Depois utilizá-la para filtrar as sequências dos cromossomos:
Para o arquivo com as sequências no formato FASTA:
E para o arquivo GFF/GTF:
Montagem dos transcritomas
Para seguir, é necessário ter as sequências das leituras alinhadas com o genoma referência. Neste caso, com o alinhador STAR.
./output/star_out_final/
├── SAMPLEA1
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
├── SAMPLEA2
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
├── SAMPLEB1
│ ├── Aligned.out.bam
│ ├── Aligned.out.sorted.bam
│ ├── Aligned.stats.txt
│ ├── samtools.merge.log.err.txt
│ ├── samtools.merge.log.out.txt
│ ├── samtools.sort.log.err.txt
│ └── samtools.sort.log.out.txt
└── SAMPLEB2
├── Aligned.out.bam
├── Aligned.out.sorted.bam
├── Aligned.stats.txt
├── samtools.merge.log.err.txt
├── samtools.merge.log.out.txt
├── samtools.sort.log.err.txt
└── samtools.sort.log.out.txt
4 directories, 28 files
Execução do cufflinks para montagem dos transcritos baseada em alinhamento de leituras (reads) com referência genômica.
O comando a seguir é um exemplo com a execução para SAMPLEA1.
Para a execução considerando todas as amostras, obter o script
rnaseq-ref.sh do GitHub.
O script auxiliar checkMinMaxIntronSize.sh baseado no script introntab.pl o qual foi modificado do script original (introntab.pl) obtido de wFleaBase.
Este número 0.9 que está na linha de comando representa o 90% percentil da distribuição de tamanhos de introns. O resultado deste script são dois valores separados por vírgula, onde o primeiro valor refere-se ao tamanho mínimo e o segundo ao tamanho máximo do intron no arquivo GFF referência. O resultado da linha de comando acima é 127,9697. Os valores obtidos podem ser aproximados, como foi feito no exemplo acima, 127 para 100 e 9697 para 10000.
Veja este link sobre a orientação das sequências em protocolo RNA-Seq fita-específica.
Com a execução do cufflinks, a expectativa é que os arquivos transcripts.gtf tenham sido gerados com sucesso.
./output/cufflinks/
├── SAMPLEA1
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
├── SAMPLEA2
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
├── SAMPLEB1
│ ├── genes.fpkm_tracking
│ ├── isoforms.fpkm_tracking
│ ├── skipped.gtf
│ └── transcripts.gtf
└── SAMPLEB2
├── genes.fpkm_tracking
├── isoforms.fpkm_tracking
├── skipped.gtf
└── transcripts.gtf
4 directories, 16 files
É importante verificar se os arquivos têm conteúdo.
Para visualizar o resultado da quantificação de expressão em FPKM para todos os genes:
… ou somente para os genes da simulação:
Com a linha de comando abaixo conseguimos acesso aos IDs dos genes referentes aos RNAs selecionados na simulação.
Vamos utilizar essa lista de genes para filtrar o transcripts.gtf gerado com cufflinks.
Podemos montar uma outra lista, agora com números de acesso RefSeq dos RNAs selecionados anteriormente. Este ID é próprio do GFF e não corresponde neste caso aos números de acesso RefSeq.
… usando essa lista de IDs de RNAs (do GFF/GTF) para obter as coordenadas desejadas.
Montagem do transcritoma referência
O transcritoma referência será constituído por uma fusão dos transcritomas individuais gerados com cufflinks. O programa cuffmerge é responsável por esta fusão.
Para obter uma contagem de ocorrências de montagem de acordo com o código de classe do cufflinks, encontrado na documentação do programa cuffcompare:
Quantificação baseada na nova referência (merged.gtf)
Para a quantificação de cada uma das réplicas executamos
o programa cuffquant. Os arquivos são obtidos a partir do
alinhamento realizado anteriormente com STAR (* verificar utilização de ) ou TopHat e :
Ele resultará em arquivos .cxb que poderão ser repassados
ao programa cuffdiff, para o cálculo da expressão gênica diferencial,
ou ao programa cuffnorm, para a obtenção de tabelas de contagens em formato texto.
Obtenção das matrizes de contagem
O programa cuffnorm será executado para a obtenção das matrizes de contagens de reads por gene/isoforma/cds/tss e normalizadas.
Para obter os dados de contagem sem a normalização realizada pelo cuffnorm, é necessário utilizar o script de-normalize-cuffnorm.R.
Análises de Expressão Gênica Diferencial
A expressão diferencial será realizada por gene/isoforma/cds/tss e normalizadas a partir do resultado do cuffquant com o programa cuffdiff.
Acrescentando informações aos resultados de expressão gênica diferencial
A intenção aqui é mostrar um método para acrescentar informações, ou fundir dados, baseado em coluna comum (usando um IDentificador comum) entre os dois arquivos formatados como tabelas, ou matrizes, ou, utilizando uma terminologia de R, data.frames.
Neste exemplo, a seguir, vamos fundir os dados resultantes do cuffdiff (gene_exp.diff) com descrições do gene obtidas de uma consulta do NCBI.
Para a consulta ao NCBI, vamos utilizar diretamente via linha de comando,
as ferramentas E-utilities do NCBI.
A seguir, será realizada a busca pelos genes de Trichoderma harzianum CBS 226.95 utilizando
o Taxonomy ID específico desta linhagem (983964):
Nem todas as colunas são necessárias, podemos selecionar somente as de interesse. Neste caso, selecionaremos somente as colunas 3, 6 e 8, respectivamente, GeneID, Symbol e description, para um novo arquivo ./refs/gene_info.txt.
O cufflinks/cuffmerge, durante o seu processamento, pode fundir genes considerados inicialmente distintos, ou seja, que possuem IDs distintos (exemplo: M431DRAFT_12071 e M431DRAFT_502167). Isso pode ser observado no arquivo gene_exp.diff. Neste caso, o script mergeR.R, criado para realizar a tarefa final de fundir os dados, não irá funcionar, afinal, na linha correspondente, os identificadores estarão unidos por vírgula (exemplo: M431DRAFT_12071,M431DRAFT_502167) impossibilitando o relacionamento. Essa impossibilidade se dá pelo fato de que no arquivo das descrições os dois identificadores estarão cada um em uma linha distinta no arquivo ./refs/gene_info.txt.
Após gerar o novo arquivo, executar a tarefa de fusão dos dados tabulares. Note que utilizaremos o parâmetro –all.x para que todas as linhas do arquivo ./output/cuffdiff/gene_exp.diff.spplitted.txt estejam no resultado final, no arquivo ./output/cuffdiff/gene_exp.diff.spplitted.merged.txt, mesmo que não tenha uma linha correspondente em ./refs/gene_info.txt. Essa última situação pode acontecer quando o cufflinks cria um novo mapeamento gênico, que não corresponde a nenhum outro contido no GFF/GTF original, usado como guia para alinhamento e montagem do transcritoma.
Montagem de novo de transcriptomas
Nesta prática de montagens utilizaremos o programa Trinity.
Os dados de entrada para este programa serão os mesmos utilizados com o protocolo baseado em referência usando Cufflinks (Tuxedo). Esses arquivos processados estão no diretório ./output/processed/prinseq.
Para preparar os dados para montagem com Trinity além de construir a linha de comando, foi desenvolvido o script rnaseq-denovo.sh para realizar a montagem de novo.
O script rnaseq-denovo.sh vai realizar um processo inicial de renomear as leituras para que os cabeçalhos contenham o sufixo /1 e /2, respectivamente para as reads R1 e R2.
Tanto as sequências pareadas quanto as que perderam um dos membros do par (singletons) serão renomeadas no script.
O script rnaseq-denovo.sh assume que as sequências foram previamente processadas com prinseq, e contêm os sufixos _1.fastq/_2.fastq e _1_singletons.fastq/_2_singletons.fastq
O pipeline de montagem com Trinity funciona nas versões atuais com python 3.6, portanto, é necessário carregar o ambiente com esta versão de python, antes de chamar o pipeline rnaseq-denovo.sh.
O Trinity também tem uma opção para executar a montagem utilizando como guia os alinhamentos das leituras X genoma referência. Neste caso, as entradas são os alinhamentos obtidos com STAR, nas etapas anteriores quando executadas os procedimentos seguindo o protocolo baseado em referência (Tuxedo). Para executar este modo de montagem, foi criado o script rnaseq-ref-trinity.sh. Neste script, o diretório que serve como entrada é o ./output/star_out_final, onde estão os arquivos com nome Aligned.out.sorted.bam. A primeira etapa deste script antes de realizar a montagem é fundir todos esses alinhamentos em um único arquivo (./output/All.sorted.bam).
Para avaliação das montagens
Indexação do transcriptoma utilizado para a simulação:
Para conferir se esta indexação foi bem sucedida, neste caso, aparecerão os seguintes arquivos:
- RefTrans.nhr
- RefTrans.nin
- RefTrans.nsq
Para executar a comparação vamos utilizar o blastn:
Primeiro, para a montagem com Trinity de novo:
Visualizar os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos:
Contar quantas ocorrências de sequências da base de dados (subject) foram encontradas:
Depois, para a montagem com Trinity Genome Guided:
Visualizar os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos:
Contar quantas ocorrências de sequências da base de dados (subject) foram encontradas:
ATENÇÃO: Se o número de transcritos montados for menor do que o esperado pela simulação, reveja os parâmetros, por exemplo, o tamanho mínimo de contigs gerados por padrão no script é de 600 (–min_contig_length 600), e se o tamanho do transcrito esperado for menor que este tamanho, ele será omitido do resultado final.
Para estimar a abundância e encontrar os genes/isoformas diferencialmente expressos
A seguir, vamos obter os valores de Log Fold-Change e os p-valores ajustados referentes às comparações pareadas entre SAMPLEA e SAMPLEB.
Para isso utilizaremos scripts auxiliares do programa Trinity.
Para começar a análise, a variável ambiente do linux ${TRINITY_HOME} deve apontar para o caminho do diretório base do programa Trinity.
Para checar:
…e para testar se os scripts auxiliares align_and_estimate_abundance.pl e abundance_estimates_to_matrix.pl funcionam e verificarmos as opções de parâmetros:
Primeiro este:
… depois este:
Além disso utilizaremos o script run-DESeq2.R para executar uma análise comparativa pareada de expressão gênica diferencial.
Verificar se ele está disponível para execução:
Para executar esse protocolo de análise de expressão gênica diferencial de uma forma automatizada, foi criado o script rnaseq-trinity-abundance.sh.
Para visualizar o resultado da anális com DESeq2 de modo alinhado:
Acessando o RStudio para algumas análises dos resultados
Para acessar a interface do RStudio no servidor utilizar este link.
Na interface do RStudio inserir o script para selecionar os genes diferencialmente expressos:
Troque o caminho da pasta ~/ pelo nome da sua pasta contendo as análises.
Há duas versões utilizando os dados padronizados por linha (Z-score). Na primeira versão, acontece primeiro a padronização manual dos valores da estimativa de abundância em log2. Em seguida, há a definição dos valores onde haverá transições de cores.
A segunda versão é mais automatizada, e, portanto, não é preciso especificar as transições. Dessa forma, ela é mais rápida de escrever, porém não permite uma especificação mais precisa das quebras para a transição de cores.
