Análise de transcriptome de alto rendimento para investigar interações host-pathogen

Summary

O protocolo aqui apresentado descreve um pipeline completo para analisar dados de transcriptome de sequenciamento de RNA desde leituras brutas até análises funcionais, incluindo etapas de controle de qualidade e pré-processamento para abordagens analíticas estatísticas avançadas.

Abstract

Patógenos podem causar uma grande variedade de doenças infecciosas. Os processos biológicos induzidos pelo hospedeiro em resposta à infecção determinam a gravidade da doença. Para estudar tais processos, os pesquisadores podem usar técnicas de sequenciamento de alto rendimento (RNA-seq) que medem as mudanças dinâmicas do transcriptome hospedeiro em diferentes estágios de infecção, desfechos clínicos ou gravidade da doença. Essa investigação pode levar a uma melhor compreensão das doenças, bem como a descoberta de potenciais alvos e tratamentos medicamentosos. O protocolo aqui apresentado descreve um pipeline completo para analisar dados de sequenciamento de RNA desde leituras brutas até análises funcionais. O gasoduto é dividido em cinco etapas: (1) controle de qualidade dos dados; (2) mapeamento e anotação de genes; (3) análise estatística para identificar genes expressos diferencialmente e genes co-expressos; (4) determinação do grau molecular da perturbação das amostras; e (5) análise funcional. A etapa 1 remove artefatos técnicos que podem afetar a qualidade das análises a jusante. Na etapa 2, os genes são mapeados e anotados de acordo com os protocolos padrão da biblioteca. A análise estatística na etapa 3 identifica genes que são expressos diferencialmente ou co-expressos em amostras infectadas, em comparação com as não infectadas. A variabilidade da amostra e a presença de potenciais outliers biológicos são verificadas usando o grau molecular de abordagem de perturbação na etapa 4. Por fim, a análise funcional na etapa 5 revela os caminhos associados ao fenótipo da doença. O pipeline apresentado tem como objetivo apoiar pesquisadores através da análise de dados RNA-seq a partir de estudos de interação hospedeiro-patógeno e impulsionar futuros experimentos in vitro ou in vivo , que são essenciais para entender o mecanismo molecular das infecções.

Introduction

Arboviroses, como dengue, febre amarela, chikungunya e zika, têm sido amplamente associadas a diversos surtos endêmicos e surgiram como um dos principais patógenos responsáveis por infectar humanos nas últimas ^décadas1,2. Indivíduos infectados com o vírus chikungunya (CHIKV) frequentemente têm febre, dor de cabeça, erupção cutânea, poliartalgia e ^artrite3,4,5. Vírus podem subverter a expressão genética da célula e influenciar várias vias de sinalização de hospedeiro. Recentemente, estudos de transcrição sanguínea utilizaram RNA-seq para identificar os genes expressos diferencialmente (DEGs) associados à infecção aguda de CHIKV em comparação com convalescência6 ou controles ^saudáveis7. Crianças infectadas pelo CHIKV tinham genes regulamentados que estão envolvidos na imunidade inata, como os relacionados a sensores celulares para RNA viral, sinalização JAK/STAT e vias de sinalização de receptores semelhantes a ^pedágio6. Adultos gravemente infectados com CHIKV também mostraram indução de genes relacionados à imunidade inata, como os relacionados a monócitos e ativação celular dendrítica, e às respostas ^antivirais7. As vias de sinalização enriquecidas com genes regulados incluíam as relacionadas à imunidade adaptativa, como ativação celular T e diferenciação e enriquecimento em células T e ^B7.

Vários métodos podem ser usados para analisar dados de transcriptome de genes hospedeiros e patógenos. Muitas vezes, a preparação da biblioteca RNA-seq começa com o enriquecimento de transcrições maduras poli-A. Esta etapa remove a maior parte do RNA ribossômico (rRNA) e, em alguns dos casos, RNAs virais/bacterianas. No entanto, quando a questão biológica envolve a detecção da transcrição do patógeno e o RNA são sequenciados independente da seleção anterior, muitas outras transcrições diferentes poderiam ser detectadas por sequenciamento. Por exemplo, os mRNAs subgênero têm se mostrado um fator importante para verificar a gravidade das ^doenças8. Além disso, para certos vírus, como CHIKV e SARS-CoV-2, até mesmo bibliotecas poli-A enriquecidas geram leituras virais que podem ser utilizadas em análises a jusante9,10. Quando focados na análise do transcriptome hospedeiro, os pesquisadores podem investigar a perturbação biológica entre amostras, identificar genes expressos diferencialmente e caminhos enriquecidos e gerar módulos de co-expressão7,11,12. Este protocolo destaca análises de transcriptome de pacientes infectados pelo CHIKV e indivíduos saudáveis utilizando diferentes abordagens bioinformáticas (Figura 1A). Foram utilizados dados de um estudo publicado anteriormente^{, composto} por 20 indivíduos saudáveis e 39 CHIKV com infecções agudas para gerar os resultados representativos.

Protocol

As amostras utilizadas neste protocolo foram aprovadas pelos comitês de ética tanto do Departamento de Microbiologia do Instituto de Ciências Biomédicas da Universidade de São Paulo quanto da Universidade Federal de Sergipe (Protocolos: 54937216.5.0000.5467 e 54835916.2.0000.5546, respectivamente).

1. Instalação de desktop docker

NOTA: As etapas para preparar o ambiente Docker são diferentes entre os sistemas operacionais (OSs). Portanto, os usuários de Mac devem seguir as etapas listadas como 1.1, os usuários de Linux devem seguir as etapas listadas como 1.2 e os usuários do Windows devem seguir as etapas listadas como 1.3.

1. Instale no MacOS.
   1. Acesse o site Get Docker (Tabela de Materiais), clique no Docker Desktop para Mac e clique no link Download do Docker Hub .
   2. Baixe o arquivo de instalação clicando no botão Obter Docker .
   3. Execute o arquivo .dmg Docker para abrir o instalador e, em seguida, arraste o ícone para a pasta Aplicativos . Localize e execute o Docker.app na pasta Aplicativos para iniciar o programa.
      NOTA: O menu específico do software na barra de status superior indica que o software está em execução e que ele está acessível a partir de um terminal.
2. Instale o programa de contêineres no Sistema Operacional Linux.
   1. Acesse o site Get Docker Linux (Tabela de Materiais) e siga as instruções para instalar usando a seção de repositório disponível no link do Repositório Linux Docker .
   2. Atualize todos os pacotes Linux usando a linha de comando:
      sudo apt-get atualização
   3. Instale os pacotes necessários no Docker:
      sudo apt-get instalar apt-transport-https ca-certificates curl gnupg lsb-release
   4. Crie um arquivo de chaveamento de arquivo de software:
      curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo gpg --dearmor -o /usr/share/keyrings/docker-archive-keyring.gpg
   5. Adicione informações do Docker deb no arquivo source.list:
      echo "deb [arch=amd64 signed-by=/usr/share/keyrings/docker-archive-keyring.gpg] https://download.docker.com/linux/ubuntu $(lsb_release -cs) estável" | sudo tee /etc/apt/sources.list.d/docker.list > /dev/null
   6. Atualize todos os pacotes novamente, incluindo os recentemente adicionados:
      sudo apt-get atualização
   7. Instale a versão desktop:
      sudo apt-get instalar docker-ce docker-ce-cli containerd.io
   8. Selecione a área geográfica e o fuso horário para terminar o processo de instalação.
3. Instale o programa de contêiner no sistema operacional Windows.
   1. Acesse o site get docker (Tabela de Materiais) e clique em Get Started. Encontre o instalador do Docker Desktop para Windows. Baixe os arquivos e instale-os localmente no computador.
   2. Após o download, inicie o arquivo de instalação (.exe) e mantenha os parâmetros padrão. Certifique-se de que as duas opções Instale componentes necessários do Windows para WSL 2 e Adicione atalho à área de trabalho estão marcadas.
      NOTA: Em alguns casos, quando este software tenta iniciar o serviço, ele mostra um erro: a instalação do WSL está incompleta. Para descobrir esse erro, acesse o site WSL2-Kernel (Tabela de Materiais).
   3. Baixe e instale o mais recente kernel WSL2 Linux.
   4. Acesse o terminal PowerShell como administrador e execute o comando:
      dism.exe /online /enable-feature /featurename:Microsoft-Windows-Subsystem-Linux /all /norestart
   5. Certifique-se de que o software Docker Desktop seja instalado com sucesso.
4. Baixe a imagem do repositório da CSBL no hub Docker (Tabela de Materiais).
   1. Abra a área de trabalho do Docker e verifique se o status está "funcionando" no canto inferior esquerdo da barra de ferramentas.
   2. Vá para a linha de comando do terminal Windows PowerShell. Baixe a imagem do Contêiner Linux para este protocolo no repositório da CSBL no hub Docker. Execute o seguinte comando para baixar a imagem:
      docker puxar csblusp/transcriptome
      NOTA: Depois de baixar a imagem, o arquivo pode ser visto no Desktop Docker. Para criar o contêiner, os usuários do Windows devem seguir o passo 1.5, enquanto os usuários de Linux devem seguir o passo 1.6.
5. Inicialize o contêiner do servidor no sistema operacional Windows.
   1. Visualize o arquivo de imagem Do Docker no gerenciador de aplicativos de desktop da barra de ferramentas e acesse a página Imagens.
      NOTA: Se a imagem do pipeline foi baixada com sucesso, haverá uma imagem csblusp/transcriptome disponível.
   2. Inicie o recipiente a partir da imagem csblusp/transcriptome clicando no botão Executar . Expanda as Configurações Opcionais para configurar o recipiente.
   3. Defina o Nome do Contêiner (por exemplo, servidor).
   4. Associe uma pasta no computador local com a pasta dentro do docker. Para fazer isso, determine o Caminho do Host. Defina uma pasta na máquina local para armazenar os dados processados que serão baixados no final. Defina o caminho do contêiner. Defina e vincule a pasta do contêiner csblusp/transcriptome ao caminho da máquina local (use o nome "/opt/transferdata" para o Caminho do Contêiner).
   5. Depois disso, clique em Executar para criar o recipiente csblusp/transcriptome.
   6. Para acessar o terminal Linux a partir do recipiente csblusp/transcriptome, clique no botão CLI.
   7. Digite o terminal bash para ter uma experiência melhor. Para isso, execute o comando:
      Bash
   8. Após executar o comando bash, certifique-se de que o terminal mostra (root@:/#):
      root@ac12c583b731:/ #
6. Inicialize o contêiner do servidor para o Sistema Operacional Linux.
   1. Execute este comando para criar o contêiner Docker com base na imagem:
      docker run -d -it --rm --name server -v :/opt/transferdata csblusp/transcriptome
      NOTA: : defina um caminho da máquina de pasta local.
   2. Execute este comando para acessar o terminal de comando do contêiner Docker:
      docker exec -it servidor bash
   3. Garanta a disponibilidade de um terminal Linux para executar quaisquer programas/scripts usando a linha de comando.
   4. Após executar o comando bash, certifique-se de que o terminal mostra (root@:/#):
      root@ac12c583b731:/ #
      NOTA: A senha raiz é "transcriptome" por padrão. Se desejar, a senha raiz pode ser alterada executando o comando:
      passwd
   5. Primeiro, execute o comando de origem para addpath.sh para garantir que todas as ferramentas estejam disponíveis. Execute o comando:
      fonte /opt/addpath.sh
7. Verifique a estrutura da pasta de sequenciamento RNA.
   1. Acesse a pasta de scripts do pipeline transcriptome e garanta que todos os dados do sequenciamento do RNA sejam armazenados dentro da pasta: /home/transcriptome-pipeline/data.
   2. Certifique-se de que todos os resultados obtidos a partir da análise sejam armazenados dentro da pasta do caminho /home/transcriptome-pipeline/resultados.
   3. Certifique-se de que os arquivos de referência de genoma e anotação sejam armazenados dentro da pasta do caminho /home/transcriptome-pipeline/datasets. Esses arquivos ajudarão a suportar todas as análises.
   4. Certifique-se de que todos os scripts sejam armazenados na pasta do caminho /home/transcriptome-pipeline/scripts e separados por cada passo conforme descrito abaixo.
8. Baixe a anotação e o genoma humano.
   1. Acesse a pasta scripts:
      cd /home/transcriptome-pipeline/scripts
   2. Execute este comando para baixar o genoma humano de referência:
      bater downloadGenome.sh
   3. Para baixar a anotação, execute o comando:
      bater downloadAnnotation.sh
9. Alterar a anotação ou a versão do genoma de referência.
   1. Abra downloadAnnotation.sh e downloadGenome.sh para alterar a URL de cada arquivo.
   2. Copie os arquivos downloadAnnotation.sh e downloadGenome.sh para a área de transferência e edite no sistema operacional local.
      cd /home/transcriptome-pipeline/scripts
      cp downloadAnnotation.sh downloadGenome.sh /opt/transferdata
   3. Abra a pasta Host Path , que é selecionada para vincular entre host e contêiner Docker na etapa 1.5.4.
   4. Edite os arquivos usando o software de editor preferido e salve. Finalmente, coloque os arquivos modificados na pasta de script. Execute o comando:
      cd /opt/transferdata
      cp downloadAnnotation.sh downloadGenome.sh /home/transcriptome-pipeline/scripts
      NOTA: Esses arquivos podem ser editados diretamente usando o editor vim ou nano Linux.
10. Em seguida, configure a ferramenta fastq-dump com a linha de comando:
    vdb-config --interativo
    NOTA: Isso permite baixar arquivos de sequenciamento a partir dos dados de exemplo.
    1. Navegue pela página Ferramentas usando a tecla guia e selecione a opção pasta atual. Navegue até a opção Salvar e clique em OK. Em seguida, saia da ferramenta fastq-dump.
11. Inicie o download das leituras do artigo publicado ^{anteriormente7}. É necessário o número de adesão da SRA de cada amostra. Obtenha os números da SRA no site da SRA NCBI (Tabela de Materiais).
    NOTA: Para analisar os dados do RNA-Seq disponíveis em bancos de dados públicos, siga o passo 1.12. Para analisar os dados privados do RNA-seq, siga o passo 1.13.
12. Analisar dados públicos específicos.
    1. Acesse o site do Centro Nacional de Informações de Biotecnologia (NCBI) e busque palavras-chave para um assunto específico.
    2. Clique no link Resultado para BioProject na seção Genomas .
    3. Escolha e clique em um estudo específico. Clique nos Experimentos SRA. Uma nova página é aberta, que mostra todas as amostras disponíveis para este estudo.
    4. Clique no número de adesão "Enviar para:" acima. Na opção "Escolher destino" selecione a opção Arquivo e Formato , selecione ExecutarInfo. Clique em "Criar arquivo" para exportar todas as informações da biblioteca.
    5. Salve o arquivo .csv SraRunInfo no caminho Host definido na etapa 1.5.4 e execute o script de download:
       cp /opt/transferdata/SraRunInfo.csv /home/transcriptome-pipeline/data
       cd /home/transcriptome-pipeline/scripts
       bater downloadAllLibraries.sh
13. Analise dados de sequenciamento privados e inéditos.
    1. Organize os dados de sequenciamento em uma pasta chamada Reads.
       NOTA: Dentro da pasta Leituras , crie uma pasta para cada amostra. Essas pastas devem ter o mesmo nome para cada amostra. Adicione dados de cada amostra dentro de seu diretório. Caso seja um RNA-Seq de ponta emparelhado, cada diretório de amostra deve conter dois arquivos FASTQ, que devem apresentar nomes terminando de acordo com os padrões {sample}_1.fastq.gz e {sample}_2.fastq.gz, sequências para frente e reversa, respectivamente. Por exemplo, uma amostra chamada "Healthy_control" deve ter um diretório com o mesmo nome e arquivos FASTQ chamados Healthy_control_1.fastq.gz e Healthy_control_2.fastq.gz. No entanto, se o sequenciamento da biblioteca for uma estratégia de fim único, apenas um arquivo de leitura deve ser salvo para análise a jusante. Por exemplo, a mesma amostra, "Controle saudável", deve ter um arquivo FASTQ exclusivo chamado Healthy_control.fastq.gz.
    2. Crie um arquivo fenotípico contendo todos os nomes da amostra: Nomeie a primeira coluna como 'Amostra' e a segunda coluna como 'Classe'. Encha a coluna Amostra com nomes de amostra, que devem ser o mesmo nome para os diretórios amostrais e preencha a coluna Classe com o grupo fenotípico de cada amostra (por exemplo, controle ou infectado). Finalmente, salve um arquivo com o nome "metadados.tsv" e envie-o para o /home/transcriptome-pipeline/data/diretório. Confira os metadados.tsv existentes para entender o formato do arquivo fenotípico.
       cp /opt/transferdata/metadata.tsv
       /home/transcriptome-pipeline/data/metadata.tsv
    3. Acesse o diretório Host Path definido na etapa 1.5.4 e copie as novas amostras de diretórios estruturados. Finalmente, mova as amostras de /opt/transferda para o diretório de dados do pipeline.
       cp -rf /opt/transferdata/reads/*
       /home/transcriptome-pipeline/data/reads/
14. Observe que todas as leituras estão armazenadas na pasta /home/transcriptome-pipeline/data/reads.

2. Controle de qualidade dos dados

NOTA: Avalie, graficamente, a probabilidade de erros nas leituras de sequenciamento. Remova todas as sequências técnicas, por exemplo, adaptadores.

1. Acesse a qualidade de sequenciamento das bibliotecas com a ferramenta FastQC.
   1. Para gerar os gráficos de qualidade, execute o programa fastqc. Execute o comando:
      bater FastQC.sh
      NOTA: Os resultados serão salvos na pasta /home/transcriptome-pipeline/results/FastQC. Uma vez que adaptadores de sequência são usados para preparação e sequenciamento da biblioteca, em alguns casos os fragmentos de sequência de adaptadores podem interferir no processo de mapeamento.
2. Remova a sequência do adaptador e as leituras de baixa qualidade. Acesse a pasta Scripts e execute o comando para a ferramenta Trimmomatic:
   cd /home/transcriptome-pipeline/scripts
   bater trimmomatic.sh
   NOTA: Os parâmetros utilizados para o filtro de sequenciamento são: Remover a de baixa qualidade ou 3 bases (abaixo da qualidade 3) (LEADING:3); Remover a baixa qualidade ou 3 bases (abaixo da qualidade 3) (TRAILING:3); Escaneie a leitura com uma janela deslizante de 4 bases, cortando quando a qualidade média por base cair abaixo de 20 (DESLIZANDO:4:20); e Drop lê abaixo as 36 bases de comprimento (MINLEN:36). Esses parâmetros podem ser alterados editando o arquivo de script Trimmomatic.
   1. Certifique-se de que os resultados sejam salvos na seguinte pasta: /home/transcriptome-pipeline/results/trimreads. Execute o comando:
      ls /home/transcriptome-pipeline/results/trimreads

3. Mapeamento e anotação de amostras

NOTA: Após a obtenção das leituras de boa qualidade, estas precisam ser mapeadas para o genoma de referência. Para esta etapa, o mapeador STAR foi usado para mapear as amostras de exemplo. A ferramenta de mapeador STAR requer 32 GB de memória RAM para carregar e executar as leituras e mapeamento do genoma. Para usuários que não possuem 32 GB de memória RAM, leituras já mapeadas podem ser usadas. Nesses casos, pule para o passo 3.3 ou use o mapeador Bowtie2. Esta seção tem scripts para STAR (resultados mostrados em todas as figuras) e Bowtie2 (mapeador necessário de baixa memória).

1. Primeiro indexe o genoma de referência para o processo de mapeamento:
   1. Acesse a pasta Scripts usando a linha de comando:
      cd /home/transcriptome-pipeline/scripts
   2. Para o mapeador STAR, execute:
      bater indexGenome.sh
   3. Para Bowtie mapper, execute:
      bater indexGenomeBowtie2.sh
2. Execute o seguinte comando para mapear leituras filtradas (obtidas a partir da etapa 2) para o genoma de referência (versão GRCh38). Mappers STAR e Bowtie2 são realizados usando parâmetros padrão.
   1. Para o mapeador STAR, execute:
      mapSTAR.sh bater
   2. Para Bowtie2 mapper, execute:
      bater mapBowtie2.sh
      NOTA: Os resultados finais são arquivos BAM (Binary Alignment Map, mapa de alinhamento binário) para cada amostra armazenada em /home/transcriptome-pipeline/results/mapreads.
3. Anotar leituras mapeadas usando a ferramenta FeatureCounts para obter contagem bruta para cada gene. Execute os scripts que anotam as leituras.
   NOTA: A ferramenta FeatureCounts é responsável por atribuir leituras de sequência mapeadas aos recursos genômicos. Os aspectos mais importantes da anotação do genoma que podem ser alterados após a questão biológica incluem, detecção de isoformes, leituras mapeadas múltiplas e junções exon-exon, correspondentes aos parâmetros, GTF.attrType="gene_name" para gene ou não especificar os parâmetros para o nível de meta-recurso, permiteMultiOverlap=TRUE e juncCounts=TRUE, respectivamente.
   1. Acesse a pasta scripts usando a linha de comando:
      cd /home/transcriptome-pipeline/scripts
   2. Para anotar as leituras mapeadas para obter contagens brutas por gene, execute a linha de comando:
      Anotação de Rscript. R
      NOTA: Os parâmetros utilizados para o processo de anotação foram: nome curto de gene de retorno (GTF.attrType="gene_name"); permitir múltiplas sobreposições (permitirMultiOverlap = TRUE); e indicar que a biblioteca é emparelhada (isPairedEnd=TRUE). Para uma estratégia única, use o parâmetro isPairedEnd=FALSE. Os resultados serão salvos na pasta /home/transcriptome-pipeline/countreads.
4. Normalizar a expressão genética.
   NOTA: A normalização da expressão genética é essencial para comparar resultados entre desfechos (por exemplo, amostras saudáveis e infectadas). A normalização também é necessária para realizar a co-expressão e o grau molecular de análises de perturbação.
   1. Acesse a pasta Scripts usando a linha de comando:
      cd /home/transcriptome-pipeline/scripts
   2. Normalize a expressão genética. Para isso, execute a linha de comando:
      Rscript normaliza as séries. R
      NOTA: A expressão da contagem bruta, neste experimento, foi normalizada utilizando os métodos M-values (TMM) e Count Per Million (CPM). Esta etapa visa remover diferenças na expressão genética devido à influência técnica, fazendo a normalização do tamanho da biblioteca. Os resultados serão salvos na pasta /home/transcriptome-pipeline/countreads.

4. Genes expressos diferencialmente e genes co-expressos

1. Identifique genes expressos diferencialmente usando o pacote EdgeR de código aberto. Isso envolve encontrar genes cuja expressão é maior ou menor em comparação com o controle.
   1. Acesse a pasta Scripts usando a linha de comando:
      cd /home/transcriptome-pipeline/scripts
   2. Para identificar o gene expresso diferencialmente, execute o script DEG_edgeR R usando a linha de comando:
      Rscript DEG_edgeR.R.
      NOTA: Os resultados que contenham os genes expressos diferencialmente serão salvos na pasta /home/transcriptome-pipeline/results/degs. Os dados podem ser transferidos para um computador pessoal.
2. Baixe dados do recipiente csblusp/transcriptome.
   1. Transfira dados processados do /home/transcriptome-pipeline para a pasta /opt/transferdata (computador local).
   2. Copie todos os arquivos para o computador local executando a linha de comando:
      cp -rf /home/transcriptome-pipeline/results /opt/transferdata/pipeline
      cp -rf /home/transcriptome-pipeline/data /opt/transferdata/pipeline
      NOTA: Agora, vá para o computador local para garantir que todos os resultados, conjuntos de dados e dados estejam disponíveis para download no Caminho do Host.
3. Identifique módulos de co-expressão.
   1. Acesse o site da Ferramenta de Identificação de Módulos de Co-Expressão (CEMiTool) (Tabela de
      Materiais). Esta ferramenta identifica módulos de co-expressão a partir de conjuntos de dados de expressão fornecidos pelos usuários. Na página principal, clique em Executar no canto superior direito. Isso abrirá uma nova página para carregar o arquivo de expressão.
   2. Clique em Escolher Arquivo abaixo da seção Arquivo de expressão e carregue a matriz de expressão genética normalizada 'tmm_expression.tsv' do Caminho do Host.
      NOTA: Passo 4.4. não é obrigatório.
4. Explore o significado biológico dos módulos de co-expressão.
   1. Clique em Escolher arquivo na seção Fenotipos de amostra e carregue o arquivo com fenótipos de amostra metadata_cemitool.tsv da etapa de dados de download 4.2.2. para realizar uma análise de enriquecimento de conjunto genético (GSEA).
   2. Pressione escolher arquivo na seção Interações genéticas para carregar um arquivo com interações genéticas (cemitool-interactions.tsv). É possível usar o arquivo de interações genéticas fornecidas como exemplo pelo webCEMiTool. As interações podem ser interações proteína-proteína, fatores de transcrição e seus genes transcritos, ou vias metabólicas. Esta etapa produz uma rede de interação para cada módulo de co-expressão.
   3. Clique na seção Escolher arquivo na seção Gene Sets para carregar uma lista de genes relacionados funcionalmente em um arquivo de formato Gene Matrix Transposed (GMT). O arquivo Gene Set permite que a ferramenta realize análises de enriquecimento para cada módulo de co-expressão, ou seja, uma análise de representação excessiva (ORA).
      NOTA: Esta lista de genes pode abranger caminhos, termos GO ou genes alvo de miRNA. O pesquisador pode usar os Módulos de Transcrição sanguínea (BTM) como conjunto genético para esta análise. O arquivo BTM (BTM_for_GSEA.gmt).
5. Defina parâmetros para a realização de análises de co-expressão e obtenha seus resultados.
   1. Em seguida, expanda a seção Parâmetro , clicando no sinal de mais para exibir os parâmetros padrão. Se necessário, mude-os. Verifique a caixa Aplicar VST .
   2. Escreva o e-mail na seção E-mail para receber os resultados como e-mail. Este passo é opcional.
   3. Pressione o botão Executar CEMiTool .
   4. Baixe o relatório completo de análise clicando no Relatório Completo de Download no canto superior direito. Ele vai baixar um arquivo comprimido cemitool_results.zip.
   5. Extrair o conteúdo do cemitool_results.zip com o WinRAR.
      NOTA: A pasta com o conteúdo extraído abrange vários arquivos com todos os resultados da análise e seus parâmetros estabelecidos.

5. Determinação do grau molecular de perturbação das amostras

1. Versão web do Grau Molecular de Perturbação (MDP).
   1. Para executar o MDP, acesse o site do MDP (Tabela de Materiais). MDP calcula distância molecular de cada amostra a partir da referência. Clique no botão Executar .
   2. No link Escolher arquivo , carregue o arquivo de expressão tmm_expression.tsv. Em seguida, carregue o arquivo de dados fenotípico metadata.tsv da etapa de dados do Download 4.2.2. Também é possível enviar um arquivo de anotação de via em formato GMT para calcular o escore de perturbação das vias associadas à doença.
   3. Uma vez que os dados sejam carregados, defina a coluna Class que contém as informações fenotípicas usadas pelo MDP. Em seguida, defina a classe de controle selecionando o rótulo que corresponde à classe de controle.
      NOTA: Existem alguns parâmetros opcionais que afetarão a forma como os escores amostrais são calculados. Se necessário, o usuário é capaz de alterar o método médio estatístico, o desvio padrão e a porcentagem superior dos genes perturbados.
   4. Depois disso, pressione o botão Executar MDP e os resultados do MDP serão mostrados. O usuário pode baixar os números clicando no Gráfico de Download em cada parcela, bem como a pontuação do MDP no botão Download MDP Score File .
      NOTA: Em caso de perguntas sobre como enviar os arquivos ou como o MDP funciona, basta passar pelo Tutorial e Sobre páginas da web.

6. Análise de enriquecimento funcional

1. Crie uma lista de DEGs para baixo e outra de DEGs regulamentadas. Os nomes dos genes devem estar de acordo com os símbolos genéticos de Entrez. Cada gene da lista deve ser colocado em uma linha.
2. Salve as listas genéticas no formato txt ou tsv.
3. Acesse o site do Enrichr (Tabela de Materiais) para realizar a análise funcional.
4. Selecione a lista de genes clicando no Arquivo Escolher. Selecione um dos DEGs e pressione o botão Enviar.
5. Clique em Caminhos no topo da página web para realizar análises de enriquecimento funcional com a abordagem ORA.
6. Escolha um banco de dados de caminhos. O banco de dados de caminhos "Reactome 2016" é amplamente utilizado para obter o significado biológico dos dados humanos.
7. Clique no nome do banco de dados da pathway novamente. Selecione Gráfico de barras e verifique se ele está classificado pelo ranking p-value. Se não, clique no gráfico da barra até que ele seja classificado por p-value. Este gráfico de barras inclui as 10 principais vias de acordo com os valores p.
8. Pressione o botão Configuração e selecione a cor vermelha para a análise de genes para cima ou cor azul para a análise de genes regulados. Salve o gráfico de barras em vários formatos clicando em svg, png e jpg.
9. Selecione Tabela e clique em Exportar entradas para a tabela no canto inferior esquerdo do gráfico da barra para obter os resultados da análise de enriquecimento funcional em um arquivo txt.
   NOTA: Este arquivo de resultados de enriquecimento funcional abrange em cada linha o nome de uma via, o número de genes sobrepostos entre a lista DEG submetida e a via, o valor p, valor p ajustado, razão de odds, escore combinado e o símbolo genético dos genes presentes na lista DEG que participam da via.
10. Repita os mesmos passos com a lista de outros DEGs.
    NOTA: A análise com DEGs não regulados fornece caminhos enriquecidos para genes regulados e a análise com genes up-regulated fornece caminhos enriquecidos para genes up-regulated.

Representative Results

O ambiente de computação para análises de transcriptome foi criado e configurado na plataforma Docker. Essa abordagem permite que usuários iniciantes do Linux usem sistemas de terminais Linux sem conhecimento de gerenciamento a priori. A plataforma Docker usa os recursos do sistema operacional host para criar um contêiner de serviço que inclua ferramentas específicas de usuários (Figura 1B). Um contêiner baseado na distribuição Linux OS Ubuntu 20.04 foi criado e foi totalmente configurado para análises transcriômicas, que é acessível através do terminal de linha de comando. Neste contêiner, há uma estrutura de pasta predefinida para conjuntos de dados e scripts necessários para todas as análises do pipeline (Figura 1C). Um estudo publicado pelo nosso grupo de ^pesquisa7 foi utilizado para análises, e compreendeu 20 amostras de indivíduos saudáveis e 39 amostras de indivíduos infectados agudamente (Figura 1D).

O processo de sequenciamento total do RNA pode gerar erros de leitura, que podem ser causados por um cluster com duas ou mais transcrições ou o esgotamento dos reagentes. As plataformas de sequenciamento retornam um conjunto de arquivos "FASTQ" contendo a sequência (leitura) e a qualidade associada para cada base nucleotídea (Figura 2A). A escala de qualidade phred indica a probabilidade de uma leitura incorreta de cada base (Figura 2B). Leituras de baixa qualidade podem gerar um viés ou expressão genética inadequada, desencadeando sucessivos erros para análises a jusante. Ferramentas como trimmomática foram desenvolvidas para identificar e remover leituras de baixa qualidade das amostras e aumentar a probabilidade de leituras de mapeamento (Figura 2C,D).

O módulo de mapeamento foi pré-configurado com o alinhador STAR e o hospedeiro humano GRCh38 como o genoma de referência. Nesta etapa, as leituras de alta qualidade recuperadas da etapa anterior são usadas como entrada para se alinharem ao genoma de referência humana (Figura 3A). O alinhador STAR produz um alinhamento de leituras mapeadas para um genoma de referência no arquivo formato BAM. Com base nesse alinhamento, a ferramenta FeatureCounts executa a anotação de recursos (genes) dessas leituras alinhadas usando a anotação de referência do host humano no formato de arquivo GTF (Figura 3B). Finalmente, a matriz de expressão com cada nome genético como uma linha, e cada amostra como uma coluna é gerada (Figura 3C). Um arquivo de metadados adicionais contendo os nomes da amostra e os respectivos grupos amostrais também precisa ser fornecido para uma análise mais a jusante. A matriz de expressão genética representa o número de contagens mapeadas para cada gene entre as amostras, que podem ser usadas como entrada EdgeR para identificar DEGs. Além disso, essa matriz de expressão genética foi normalizada utilizando TMM e CPM, a fim de remover a variabilidade técnica e corrigir a medição do RNA-seq, considerando a proporção de genes expressos no tamanho total da biblioteca entre as amostras. Esta matriz foi ainda utilizada como entrada para análises de co-expressão e MDP.

O CEMiTool identifica e analisa os módulos de ^{co-expressão12}. Genes que estão no mesmo módulo são co-expressos, o que significa que eles exibem padrões de expressão semelhantes nas amostras do conjunto de dados. Esta ferramenta também permite a exploração da significância biológica de cada módulo identificado. Para isso, fornece três análises opcionais - análise de enriquecimento funcional pela GSEA, análise de enriquecimento funcional por Análise de Representação Sobre Representação (ORA) e análise de rede. A análise de enriquecimento funcional da GSEA fornece informações sobre a expressão genética de cada módulo em cada fenótipo (Figura 4A). De acordo com isso, permite a identificação dos módulos que são reprimidos ou induzidos a cada fenótipo. A análise ORA mostra as 10 principais funções biológicas significativamente enriquecidas de cada módulo classificadas por valores p ajustados. É possível combinar os resultados GSEA e ORA para identificar processos biológicos prejudicados e se eles estão sendo reprimidos ou induzidos pelo fenótipo de interesse. As análises de rede fornecem um interactome de cada módulo (Figura 4A). Permite a visualização de como os genes de cada módulo interagem. Além disso, a análise da rede fornece informações sobre os genes mais conectados, os hubs, que são identificados por seus nomes na rede. O tamanho dos nódulos representa o grau de conectividade.

Para identificar OSTS, um script interno foi desenvolvido para executar uma análise diferencial de ponta a ponta em uma linha de comando única e concisa. O script executa todas as etapas necessárias para realizar uma análise DEG, comparando diferentes grupos de amostras fornecidos pelo usuário em um arquivo de metadados. Além disso, os resultados do DEG são armazenados em listas separadas de genes não regulados e up-regulated e, em seguida, compilados em uma figura pronta para publicação (Figura 4B) usando pacote EnhancedVolcano R do Bioconductor.

A análise do grau molecular de perturbação realizada pela ferramenta MDP permite identificar amostras perturbadas de indivíduos saudáveis e ^infectados11. O escore de perturbação é calculado considerando todos os genes expressos para cada amostra infectada pelo CHIKV e considerando as amostras saudáveis como o grupo de referência (Figura 5A). O MDP também realiza a análise utilizando apenas os 25% mais altos dos genes mais perturbados dessas amostras (Figura 5B). As amostras podem apresentar uma grande variabilidade, dada a formação genética, idade, sexo ou outras doenças anteriores. Esses fatores podem alterar o perfil do transcriptome. Com base nisso, o MDP sugere quais amostras são potenciais outliers biológicos para removê-las e melhorar os resultados a jusante (Figura 5A,B).

Uma análise de enriquecimento funcional por ORA pode ser realizada utilizando-se do Enrichr para identificar o significado biológico dos DEGs. Os resultados fornecidos com base na lista de genes regulados indicam os processos biológicos reprimidos no fenótipo estudado, enquanto os resultados fornecidos com base na lista de genes regulados apresentam os processos biológicos que são induzidos no fenótipo de interesse. Os processos biológicos mostrados no gráfico de barras gerados pelo Enrichr são os 10 principais conjuntos genéticos enriquecidos com base no ranking p-value (Figura 6).

Figura 1: Meio Ambiente Docker e estudo de exemplo. (A) A plataforma Docker usa os recursos do OS Host para criar "Containers" para o sistema Linux que contém ferramentas para análises de transcriptome. (B) O Docker Container simula um sistema Linux para executar scripts de pipeline. (C) A estrutura da pasta do pipeline de transcriptome foi criada e organizada para armazenar conjuntos de dados e scripts para análise. (D) O estudo do nosso grupo foi utilizado como exemplo de análises de transcriptome. Clique aqui para ver uma versão maior desta figura.

Figura 2: Controle de qualidade do sequenciamento. (A) O arquivo de formato FASTQ é usado para representar a qualidade da base sequencial e nucleotídea. (B) Equação de pontuação de phred, onde a cada 10 aumenta uma base de probabilidade de registro mal interpretada. (C) e (D) O Boxplot representa uma distribuição de qualidade de cada base nucleotídea antes e depois da execução trimmomática, respectivamente. Clique aqui para ver uma versão maior desta figura.

Figura 3: Processo de mapeamento e anotação da sequência à expressão da contagem genética. (A) O mapeamento consiste em alinhar a sequência da transcrição e a sequência do genoma para identificar a localização genômica. (B) As leituras mapeadas do genoma de referência são anotadas com base em sua localização genômica de sobreposição. (C) Com base nas ferramentas de arquivo de mapeamento, como featureCounts, a expressão genética é resumida. Clique aqui para ver uma versão maior desta figura.

Figura 4: Rede de genes co-expressa e análise estatística de DEGs. (A) Módulos de co-expressão baseados na expressão genética e na rede de interações proteína-proteína a partir de genes de módulos. (B) Análise estatística de indivíduos infectados e saudáveis do CHIKV e expressão genética diferencial em vermelho (critérios de p-valor e log2FC), roxo (apenas p-valor), verde (apenas log2FC) e cinza (sem significado). Clique aqui para ver uma versão maior desta figura.

Figura 5: Grau Molecular de Perturbação (MDP) de indivíduos infectados agudamente e saudáveis. (A) Escore de MDP para cada amostra usando todos os genes expressos do transcriptome. (B) Escore de MDP para cada amostra utilizando apenas os 25% mais altos dos genes mais perturbados. Clique aqui para ver uma versão maior desta figura.

Figura 6: A análise funcional dos GENES REGULAdos (A) e (B) regulados foi submetida à ferramenta do site Enrichr para avaliar caminhos biológicos ou conjuntos genéticos representativos. Os valores P foram calculados para cada via e apenas diferenças significativas foram mostradas no gráfico. Clique aqui para ver uma versão maior desta figura.

Discussion

A preparação das bibliotecas de sequenciamento é um passo crucial para responder às questões biológicas da melhor maneira possível. O tipo de transcrição de interesse do estudo orientará qual tipo de biblioteca de sequenciamento será escolhida e conduzirá análises bioinformáticas. Por exemplo, a partir do sequenciamento de um patógeno e interação do host, de acordo com o tipo de sequência, é possível identificar sequências de ambas ou apenas das transcrições do host.

O equipamento de sequenciamento de última geração, por exemplo, a Plataforma Descolumina, mede os escores de qualidade de sequenciamento, o que representa a probabilidade de que uma base seja chamada incorretamente. As análises a jusante são muito sensíveis a sequências de baixa qualidade e levam à expressão genética mal lida ou mal interpretada. Outro obstáculo na realização de análises e interpretação corretas são as sequências de adaptadores. As sequências de adaptadores ajudam na preparação e sequenciamento da biblioteca, e na maioria dos casos, os adaptadores também são sequenciados. Estudos recentes identificaram que o impacto da ferramenta de mapeamento nos resultados finais é ^mínimo13. No entanto, em estudos de patógenos-hospedeiros, o processo de mapeamento pode gerar resultados ligeiramente melhores ao testar diferentes limiares para minimizar o problema de sequências de lócus multi mapeadas.

Os resultados diferenciais da expressão genética devem ser interpretados com certa cautela, especialmente quando o número de amostras por grupo é muito pequeno e as amostras vieram de diferentes ensaios e interferindo por efeitos em lote que os DEGs resultam. Esses resultados são sensíveis a vários fatores: (i) a filtragem de dados aplicada, como a remoção de genes de baixa expressa e o número de amostras a serem mantidas; (ii) desenho do estudo, para comparar apenas entre grupos amostrais ou cada paciente infectado versus todos os pacientes de controle, conforme ilustrado no estudo ^CHIKV7; e (iii) método estatístico utilizado para identificar DEGs. Aqui, ilustramos um exemplo básico com o EdgeR para identificar DEGs assumindo um valor p limiar de 0,05. Também é conhecido na literatura que, em comparação com outros métodos de benchmark, o EdgeR pode ter uma grande variedade de variabilidade na identificação de ^DEGs14. Pode-se considerar a troca entre tais métodos diferentes e levar em conta o número de réplicas disponíveis e a complexidade do design ^{experimental14}.

O CEMiTool realiza análises de módulos de ^{co-expressão12}. Esta ferramenta está disponível através do pacote R no repositório biocondutor e também está disponível em uma versão fácil de usar através do webCEMiTool; Esta última é a versão usada neste protocolo atual. Trata-se de um software alternativo em relação ao ^WGCNA15 apresentando vários benefícios em relação aos ^últimos16, incluindo o fato de ser mais fácil de ^usar17. Além disso, esta ferramenta possui um método automático para filtrar genes, enquanto no WGCNA o usuário deve filtrar os genes antes do uso do WGCNA. Além disso, esta ferramenta tem parâmetros padrão estabelecidos, enquanto no WGCNA o usuário deve selecionar manualmente as análises dos parâmetros. A seleção manual de parâmetros prejudica a reprodutibilidade; portanto, a seleção de parâmetros automáticos garante melhor reprodutibilidade.

Em certos casos, o CEMiTool não é capaz de encontrar um limite macio apropriado, também chamado de valor β. Neste caso, o usuário deve verificar se os dados do RNA-seq apresentam forte dependência de variância média. Se a média apresentar uma forte relação linear com a variância (considerando todos os genes), o usuário deve refazer as análises verificando o parâmetro "Aplicar VST" para remover a dependência média de variância dos dados transcriômicos. É sempre fundamental verificar se há uma forte dependência de variância média nos dados e removê-los quando eles estão presentes.

O CEMiTool tem sido amplamente utilizado para identificar e explorar o significado biológico dos módulos de co-expressão. Um estudo de infecção aguda chikv mostrou um módulo com maior atividade em pacientes após 2 a 4 dias do início dos ^sintomas7. O enriquecimento funcional deste módulo pela ORA apresentou um aumento de monócitos e neutrófilos7. Um estudo de vacinação contra a gripe utilizando transcrição sanguínea da linha de base para o dia 7 pós-vacinação apresentou módulos de co-expressão funcionalmente enriquecidos para processos biológicos relacionados a células T, B e assassinos naturais, monócitos, neutrófilos, respostas de interferon e ativação ^{plaquetária18}.

Considerando a variabilidade dos conjuntos de dados transcriômicos, identificar e quantificar a heterogeneidade dos dados pode ser um desafio, já que muitas variáveis podem influenciar o perfil de expressão ^{genética7,11}. O MDP fornece uma maneira de identificar e quantificar amostras perturbadas de indivíduos saudáveis e infectados, seguindo estas etapas: (i) calcular um método de centralidade (mediana ou média) e desvio padrão das amostras de controle; (ii) utilizar os valores obtidos para calcular o escore z de todos os genes; (iii) estabelecer um limite z-score absoluto superior a 2, indicando desvios representativos das amostras de controle; e (iv) calcular a média dos valores genéticos utilizando os escores filtrados para cada amostra. Apesar de ter algumas limitações para a análise de scRNA-seq, esta ferramenta foi funcional na determinação do escore de perturbação a partir de microarray e dados ^RNA-seq11. Além disso, um estudo anterior usou esta ferramenta para demonstrar o grau molecular de perturbação elevado no transcritor sanguíneo em pacientes com tuberculose e diabetes ^mellitus19. Neste trabalho, mostrou-se a perturbação do controle e das amostras infectadas agudamente pelo CHIKV utilizando indivíduos saudáveis como grupo de referência.

A análise de enriquecimento funcional realizada pelo Enrichr é a ^ORA20,21. ORA é um tipo de análise de enriquecimento funcional em que o usuário deve fornecer a lista de DEGs para a ferramenta. A lista de DEGs é geralmente separada em uma lista DEG regulamentada e em uma lista DEG regulamentada. Existem outras ferramentas para executar o ORA, entre elas, o gProfiler, que está disponível em uma versão web ^amigável22 e o goseq23 que está disponível como um pacote R no Bioconductor. Outro tipo de análise de enriquecimento funcional é o GSEA. Para executar o GSEA, o usuário deve fornecer todos os genes em uma lista ranqueada. Esta lista é geralmente classificada de acordo com a expressão genética na mudança de dobra.

Enrichr sempre fornece os 10 principais conjuntos genéticos enriquecidos com base em seus valores p no resultado gráfico da barra. Portanto, o usuário deve estar atento ao interpretar os resultados, se houver menos de 10 conjuntos genéticos enriquecidos, o gráfico da barra também mostrará processos biológicos não enriquecidos. Para evitar esse erro, o usuário deve estabelecer um limite para o valor p e observar os valores p das vias antes de assumir que todos os conjuntos genéticos do gráfico da barra são enriquecidos. Além disso, o usuário deve estar ciente de que a ordem dos 10 conjuntos genéticos exibidos no gráfico da barra está de acordo com os valores p, não os valores p ajustados. Caso o usuário queira mostrar todas as vias enriquecidas em um gráfico de barras ou até mesmo reordenar de acordo com os valores p ajustados, recomenda-se que o usuário crie seu próprio gráfico de barras usando a tabela baixada. O usuário pode fazer um novo gráfico de barras usando excel ou até mesmo software R.

Materials

Name	Company	Catalog Number	Comments
CEMiTool	Computational Systems Biology Laboratory	1.12.2	Discovery and the analysis of co-expression gene modules in a fully automatic manner, while providing a user-friendly HTML report with high-quality graphs.
EdgeR	Bioconductor (Maintainer: Yunshun Chen [yuchen at wehi.edu.au])	3.30.3	Differential expression analysis of RNA-seq expression profiles with biological replication
EnhancedVolcano	Bioconductor (Maintainer: Kevin Blighe [kevin at clinicalbioinformatics.co.uk])	1.6.0	Publication-ready volcano plots with enhanced colouring and labeling
FastQC	Babraham Bioinformatics	0.11.9	Aims to provide a simple way to do some quality control checks on raw sequence data coming from high throughput sequencing
FeatureCounts	Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research	2.0.0	Assign mapped sequencing reads to specified genomic features
MDP	Computational Systems Biology Laboratory	1.8.0	Molecular Degree of Perturbation calculates scores for transcriptome data samples based on their perturbation from controls
R	R Core Group	4.0.3	Programming language and free software environment for statistical computing and graphics
STAR	Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research	2.7.6a	Aligner designed to specifically address many of the challenges of RNA-seq data mapping using a strategy to account for spliced alignments
Bowtie2	Johns Hopkins University	2.4.2	Ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences
Trimmomatic	THE USADEL LAB	0.39	Trimming adapter sequence tasks for Illumina paired-end and single-ended data
Get Docker	Docker	20.10.2	Create a bioinformatic environment reproducible and predictable (https://docs.docker.com/get-docker/)
WSL2-Kernel	Windows	NA	https://docs.microsoft.com/en-us/windows/wsl/wsl2-kernel
Get Docker Linux	Docker	NA	https://docs.docker.com/engine/install/ubuntu/
Docker Linux Repository	Docker	NA	https://docs.docker.com/engine/install/ubuntu/#install-using-the-repository
MDP Website	Computational Systems Biology Laboratory	NA	https://mdp.sysbio.tools
Enrichr Website	MaayanLab	NA	https://maayanlab.cloud/Enrichr/
webCEMiTool	Computational Systems Biology Laboratory	NA	https://cemitool.sysbio.tools/
gProfiler	Bioinformatics, Algorithmics and Data Mining Group	NA	https://biit.cs.ut.ee/gprofiler/gost
goseq	Bioconductor (Maintainer: Matthew Young [my4 at sanger.ac.uk])	NA	http://bioconductor.org/packages/release/bioc/html/goseq.html
SRA NCBI study	NCBI	NA	https://www.ncbi.nlm.nih.gov/bioproject/PRJNA507472/