In currocam/FascinRSCA: R and Python code for SCA with fascin protein

knitr::opts_chunk$set(echo = TRUE)

library("protr")
library(tidyverse)
library(magrittr)
library(FascinRSCA)
# load FASTA files
fasta_file <- protr::readFASTA("../inst/extdata/emboss/NonRedundant_FSCN_Human_Homologous.fasta")
df <- FascinRSCA::freezing_cluster %>%
  filter(fasta_header %in% names(fasta_file))
df$annotation %<>%
  fct_collapse(fascin2 = c("fascin2","fascin2a"), fascin = c("fascin", "singed", "uncharacterized"))

Actualmente tenemos las siguientes secuencias. Nos puede interesar tener el mismo número de las tres clases (por oversampling o undersampling, quizás, se puede implementar maś adelante).

df$annotation %>% fct_count()

A continuación, vamos a buscar un número igual de secuencias que no son fascina. En un principio vamos a buscar 300 secuencias, de las cuales 150 van a ser proteínas que intervienen en procesos similares o en tienen funciones parecidas y que extraeremos de distintas especies. El proceso se muestra a continuación:

Proteínas funcionalmente similares

En primer lugar, cargamos el acceso a la base de datos STRING para humanos y accedemos a las proteínas de fascina 1 y fascina 2 en humano.

library(STRINGdb)
string_db <- STRINGdb$new(version="11", species= 9606,
  score_threshold=200, input_directory = '')
string_proteins <- string_db$get_proteins()
human.fascin <- string_db$mp(c("FSCN1", "FSCN2"))

La anotación funcional de dichas proteínas es la siguiente:

`%!in%` <-  Negate(`%in%`)

human.fascin.category <- string_db$get_annotations(human.fascin) %>%
  filter(category == "Function") %>%
  filter(description %!in%c("binding",
                            "protein binding",
                            "drug binding",
                            "molecular adaptor activity",
                            "protein-containing complex binding",
                            "protein binding, bridging"))
go.terms <- human.fascin.category %>% 
  pull(term_id)

A continuación, descargamos las proteínas que contienen dichos términos y que pertenecen a alguna de nuestras especies en formato .gpad. Descargamos varios archivos porque hay un máximo de 50 organismos por búsqueda. Juntamos los resultados en el siguiente DataFrame.

requestURL <- paste0("https://www.ebi.ac.uk/QuickGO/annotations?taxonId=%20",
  paste(unique(df$taxid)[1:50], collapse = ","),                     
  "%20&taxonUsage=exact&goId=GO:0051015,GO:0008092,GO:0003779&goUsage=exact")

path <- "../inst/extdata/quickGo/fascin_func/"
files_gpad <- purrr::map_chr(dir(path, pattern = "*.gpad"), ~paste0(path, .x))

gpad <- files_gpad %>%
  purrr::map(~read.table(.x, comment = "!",
    header = FALSE, fill = TRUE,
    col.names = c("DB", "DB.ID", "Qualifier", "GO.ID", "DB:Reference",
      "Evidence.Code", "WithFrom", "InteractingTaxon", "Date", "AssignedBy",
      "AnnotationExtension", "AnnotationProperties")) %>%
    select(DB, DB.ID, Qualifier,GO.ID))%>%
  bind_rows(.id = "column_label")

A continuación, sampleamos 150 proteínas, escribimos la lista de identificadores en un archivo y descargamos las proteínas de UniProt.

gpad %>%
  distinct()%>%
  group_by(GO.ID) %>%
  slice_sample(n =50)%>%
  pull(DB.ID) %>%
write("RelateddbID.txt")

Proteínas funcionalmente dispares

A continuación, vamos a obtener 150 proteínas, de nuevo pertenecientes a las especies de nuestro dataset, de forma aleatoria.

fasta.NonRelated <- unique(df$taxid)%>%
  sample(300, replace = TRUE)%>%
  purrr::map_dfr(possibly(otherwise = tibble(seq = NA,
         name = NA),
  function(x){
  fasta_seq <- glue::glue('https://www.uniprot.org/uniprot/?query=reviewed:yes+AND+',
                    'organism:{taxid}&random=yes&format=fasta',
                    taxid = sample(x))%>%
    readFASTA()
  tibble(seq = fasta_seq %>% chuck(1),
         name = names(fasta_seq))
  })) 
fasta.NonRelated %>%
  filter(!is.na(seq))%>%
  distinct(name, .keep_all = TRUE)%>%
  slice_sample(n = 150)%$%
  seqinr::write.fasta(sequences = lapply(seq, seqinr::s2c), names = name,
  file.out = "NonRelated.fasta")