Índice
O reticulate é um pacote do R que realiza uma integração entre sessões do R e Python. Eu recentemente descobri esse pacote e estou muito feliz em executar scripts do Python sem sair do R, é uma verdadeira união do útil ao agradável.
Integrando R e Python com o reticulate
Nesse post eu vou demonstrar o reticulate com dois scripts. Primeiro, eu vou começar uma sessão no R com um script R. Em seguida, vou executar um script Python e manipular o output dele dentro da sessão do R.
O resultado da demonstração será uma tabela (data frame) contendo as coordenadas dos exons de dois genes e as sequências de nucleotídeos desses exons.
O código dessa demonstração está salvo no meu portfólio no GitHub.
Script do R (parte 1/2)
Instalando pacotes
Vou começar listando os pacotes do R que usarei na demonstração. Lembre-se de que pode instalar os pacotes no R com o comando install.packages() ou BiocManager::install() caso sejam da família Bioconductor (por sua vez o BiocManager deve ser instalado com o install.packages() também).
library(here)
library(tidyverse)
library(reticulate)
library(glue)
# Bioconductor
library(GenomicRanges)
library(BSgenome)
library(BSgenome.Hsapiens.UCSC.hg38)
Contéudo do script
Em seguida, usando o pacote here, crio o caminho completo de um arquivo GFF3 e atribuo ao objeto gffPath. O formato GFF3 serve para armazenar propriedades de genes num arquivo de texto para ajudar a representar dados genômicos. O arquivo em questão contém propriedades de todo o genoma humano. No entanto, como é um arquivo grande, não coloquei no meu repositório GitHub, mas você pode baixá-lo no site FTP do NCBI.
gffPath <- here("GCA_000001405.15_GRCh38_full_analysis_set.refseq_annotation.gff")
Edite o caminho do arquivo usando a sintaxe do here caso necessário.
Em seguida, crio um vetor simples no R com o nome de apenas dois genes para deixar a demonstração simples:
genes <- c("HTT", "FMR1")
Apresentando as funções do reticulate
Agora posso executar o script do Python. Caso você utilize ambientes virtuais conda como eu, você pode selecionar o ambiente desejado com a função use_condaenv() do reticulate. Dessa forma, usarei meu ambiente chamado de bioenv:
use_condaenv("bioenv")
Existe uma função análoga para o uso de ambientes virtualenv do Python, a use_virtualenv().
Em seguida, usando o here novamente, configuro o caminho para o script Python e depois o executo com a função source_python() do reticulate:
pythonScript <- here("reticulate_demo_Python_side.py")
Agora vou mostrar o conteúdo desse script.
Script Python
Os módulos Python que usarei nessa demonstração são o pyranges e o pandas. Acredito que o pyranges executa o pandas em segundo plano, mas como diria o Zen do Python, “explícito é melhor que implícito”.
import pyranges as pr
import pandas as pd
Instalando módulos
(Lembre-se de instalar os módulos por meio do conda ou pip):
conda install -c bioconda pyranges
# ou
pip install pyranges
Lendo o arquivo GFF3
Em seguida, carrego o arquivo GFF3 na sessão Python usando a função read_gff3() do pyranges. Já que eu atribuí o caminho do arquivo no objeto gffPath do R, devo usar o prefixo r. para trazê-lo para a sessão do Python, dessa forma:
grch38_gff = pr.read_gff3(r.gffPath)
Convertendo em data frame
Observe que atribuí o conteúdo do arquivo ao objeto grch38_gff. Nesse ponto, posso transformar o objeto num data frame do pandas para facilitar a busca pelos genes que defini anteriormente. Faço isso apenas adicionado o sufixo .df ao nome do objeto:
df = grch38_gff.df
Observe que isso é um método herdado do pyranges, nem todo objeto pode ser convertido num data frame da mesma forma. Consulte os métodos do módulo em questão para saber como fazer tal operação.
Buscando os genes
A vantagem de transformar o objeto pyranges em pandas é que eu posso encontrar os genes desejados por meio do método str.contains(). Portanto, é necessário que eu crie uma string de expressão regular (regex). Farei isso sufixando o caractere $ ao fim de cada gene. Esse caractere, indica, no regex, que a palavra buscada deve estar no fim de uma frase, o que no nosso caso, vai garantir que encontremos o nome completo do gene:
suffix = "$"
search_string = "|".join([gene + suffix for gene in r.genes])
Observe que usei uma list comprehension do Python e o método .join() para criar a string regex e atribuo a ela o objeto search_string. O caractere | (ou pipe) significa “ou” no regex. O resultado é este:
print(search_string)
# HTT$|FMR1$
Filtrando o data frame
Finalmente posso filtrar o data frame para capturar as coordenadas dos genes desejados usando a sintaxe de filtragem do pandas:
df = df[(df["gene"].str.contains(search_string)) & (df["Feature"] == "exon")]
Observe que estou sobrescrevendo o objeto df com o resultado da própria filtragem para facilitar a transferência do objeto para a sessão do R. Veja também que estou trabalhando com duas colunas distintas: a coluna gene (df["gene"]) e Feature (df["Feature"]). A primeira coluna estou filtrando com a string regex criada acima enquanto a segunda coluna estou filtrando para que contenha apenas a característica (feature) exons. Em outras palavras, esse comando selecionará as coordenadas dos exons dos genes desejados.
O comando acima conclui o script Python. Vamos retornar ao script R para terminar a demonstração do reticulate.
Script R: parte 2/2
Nós podemos trazer objetos Python para a sessão do R de maneira similar à apontada para o caminho inverso. No R usamos o prefixo py$. Dessa forma, vou trazer o objeto df da sessão Python e atribuí-lo ao objeto exonsCoordinates do R:
exonsCoordinates <- py$df
Função customizada getSequence()
Antes de continuar, permita-me trazer uma função que criei que irá recuperar as sequências exônicas usando os pacotes GenomicRanges (o pyranges é o análogo Pythonico deste), BSgenome e BSgenome.Hsapiens.UCSC.hg38:
source(here("src", "getSequence.R"))
# Código da função (arquivo getSequence.R):
getSequence <- function(chr, start, end) {
gr <- GenomicRanges::GRanges(glue::glue("{chr}:{start}-{end}"))
refBase <- BSgenome::getSeq(BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38, gr)
refBase <- as.data.frame(refBase)$x[1]
return(refBase)
}
Importante: o primeiro uso do BSgenome/BSgenome.Hsapiens.UCSC.hg38 realizará o download de arquivos de sequência do genoma humano que serão salvos numa localização específica do seu computador. Certifique-se que possui cota de transmissão de dados e espaço suficiente.
Como o nome do terceiro pacote indica, estou utilizando a versão (build) GRCh38 do genoma humano.
Finalizando
Finalmente, usando a sintaxe de pipes (%>%) da família de pacotes tidyverse (mais especificamente o dplyr), crio uma nova coluna no data frame (com a função mutate()) usando as colunas Chromosome, Start e End como argumentos para a função getSequence() mencionada acima. Essa nova coluna conterá as sequências de nucleotídeos dos exons. Observe o uso das funções rowwise() e ungroup(). Já que minha função não é vetorizada, devo usá-la linha-a-linha em vez de coluna-a-coluna como é o comportamento padrão do R. A primeira função garante isso, enquanto a segunda restaura o data frame ao comportamento padrão assim que a getSequence() termina seu trabalho. Em seguida, finalizo selecionando apenas algumas colunas do data frame com a função select():
exonsCoordinatesClean <- exonsCoordinates %>%
rowwise() %>%
mutate(sequence = getSequence(Chromosome, Start, End)) %>%
ungroup() %>%
select(Chromosome, Feature, Start, End, Score, Strand, Frame, ID, gene, sequence)
O resultado é este:
Conclusão
- Demonstrei como executar um script Python com o reticulate sem sair de uma sessão ativa do R;
- Expliquei como recuperar objetos Python em uma sessão R e vice-versa;
- Demonstrei uma função simples que recupera sequências de nucleotídeos a partir do genoma de referência humano (build GRCh38) usando coordenadas cromossômicas.
