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

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(glue)
library(magrittr)

Agrupando secuencias en archivos fasta

A continuación, se van a agrupar las distintas secuencias según su categoría taxonómica, su anotación y el cluster. Para ello se parte del dataset "freezing" redundante, y, por ejemplo, para el archivo fasta de secuencias pertenecientes al cluster 1 de mamíferos primero se filtrará aquellas que cumplan dicha condición y después se eliminarán las secuencias redundantes. De esta manera se minimiza la pérdida de información.

Para ello se hará uso de la función propia multiple_fasta_files.

df <- FascinRSCA::freezing_cluster %>%
  mutate(annotation = fct_collapse(annotation, "fascin2" = c("fascin2", "fascin2a")))# Eliminar o no

Definimos las categorías taxonómicas que nos interesan:

tax_values = c("superkingdom", "kingdom", "phylum","class",
  "order", "suborder","family", "genus")

En primer lugar, vamos a generar los archivos con todas las posibles combinaciones. Comenzamos con aquellas secuencias anotadas como fascinas 1.

df %>%  filter(annotation != "fascin2") %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = tax_values, seq_col = "fasta_seq",
                       header_col = "fasta_header", path ="fascin1", toWrite = TRUE) %>%
  purrr::map(2) %>% unlist() %>% head()

Aquellas secuencias anotadas como fascinas 2.

df %>%  filter(annotation == "fascin2") %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = tax_values, seq_col = "fasta_seq",
                       header_col = "fasta_header", path ="fascin2", toWrite = TRUE) %>%
  purrr::map(2) %>% unlist() %>% head()

Comenzamos con aquellas secuencias que fueron agrupadas en el cluster 1.

df %>%  filter(cluster == 1) %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = tax_values, seq_col = "fasta_seq",
                       header_col = "fasta_header", path ="cluster1", toWrite = TRUE) %>%
  purrr::map(2) %>% unlist() %>% head()

Aquellas secuencias anotadas como fascinas 2.

df %>%  filter(cluster == 2) %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = tax_values, seq_col = "fasta_seq",
                       header_col = "fasta_header", path ="cluster2", toWrite = TRUE) %>%
  purrr::map(2) %>% unlist() %>% head()

Alineamientos

En primer lugar, vamos a escribir los fasta con las secuencias que vamos a emplear para alinear con t_coffee. Comenzaremos con aquellas anotadas como fascina 1.

df %>%
  dplyr::mutate(phylum = forcats::fct_lump_min(phylum, 10,
    other_level = "Other_phyla"),
  class = forcats::fct_lump_min(class, 10,
    other_level = "Other_class"))%>%
  filter(annotation != "fascin2") %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = c("phylum", "class"), 
                                   seq_col = "fasta_seq",
                                   header_col = "fasta_header",
                                   path ="fascin1_lump_min",
                                   toWrite = TRUE) %>%
  purrr::map(2) %>% unlist()

Aquellas anotadas como fascina 2.

df %>%
  dplyr::mutate(phylum = forcats::fct_lump_min(phylum, 10,
    other_level = "Other_phyla"),
  class = forcats::fct_lump_min(class, 10,
    other_level = "Other_class"))%>%
  filter(annotation == "fascin2") %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = c("phylum", "class"), 
                                   seq_col = "fasta_seq",
                                   header_col = "fasta_header",
                                   path ="fascin2_lump_min",
                                   toWrite = TRUE) %>%
  purrr::map(2) %>% unlist()

Aquellos que fueron agrupados en el cluster 1:

df %>%
  dplyr::mutate(phylum = forcats::fct_lump_min(phylum, 10,
    other_level = "Other_phyla"),
  class = forcats::fct_lump_min(class, 10,
    other_level = "Other_class"))%>%
  filter(cluster == 1) %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = c("phylum", "class"), 
                                   seq_col = "fasta_seq",
                                   header_col = "fasta_header",
                                   path ="cluster1_lump_min",
                                   toWrite = TRUE) %>%
  purrr::map(2) %>% unlist()

Aquellos que fueron agrupados en el cluster 2:

df %>%
  dplyr::mutate(phylum = forcats::fct_lump_min(phylum, 10,
    other_level = "Other_phyla"),
  class = forcats::fct_lump_min(class, 10,
    other_level = "Other_class"))%>%
  filter(cluster == 2) %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = c("phylum", "class"), 
                                   seq_col = "fasta_seq",
                                   header_col = "fasta_header",
                                   path ="cluster2_lump_min",
                                   toWrite = TRUE) %>%
  purrr::map(2) %>% unlist()

Anotación de las secuencias

Cambiamos la anotación de los headers de las secuencias para que se adapte a la herramienta de pySCA. Para ello creamos una nueva columna que se llama annotated_header.

df %<>%
  unite("unite_tax", superkingdom:genus, sep = ",", na.rm = TRUE, remove = FALSE)%>%
  mutate(new_fasta_header = str_replace(fasta_header, "[|]", "_"),
         Entry = str_replace(Entry, "[|]", "_"),
         annotated_header = glue("{new_fasta_header}|{Entry}|{species}|{unite_tax}"))

df$annotated_header[[1]]

A continuación, vamos a leer los archivos alineados y a reescribirlos con el nuevo header y en formato fasta.

formato_original <- "clustal"
path <- "../inst/extdata/t_coffee"
input_format <- ".aln" # o "*.aln"
path_aln <- purrr::map_chr(dir(path, pattern = input_format), ~file.path(path, .x))#Obtenemos el path de los distintos archivos
path_aln
purrr::walk(path_aln, function(x){
  aln <- seqinr::read.alignment(file = x, format=formato_original) #Leemos el alineamiento
  names_aln <- df %>%
    filter(fasta_header %in% aln$nam) %>%
    dplyr::pull(annotated_header)

  seqinr::write.fasta(sequences = aln$seq, 
    names = names_aln,
    file.out = paste0(sub(paste0("\\", input_format, "$"), '', x), ".fa"),#Cambiamos la terminación
    nbchar = 80,open = "w")
})