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

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

En este cuaderno se muestra el proceso de obtención de un alineamiento de fascinas 1 que será utilizado en el SCA para la idetificación de sectores. El objetivo es conseguir, mediante el uso de HMM profiles, un alineamiento con mayor número de secuencias efectivas.

En primer lugar, vamos a escribir en un archivo fasta todas las secuencias de fascina 1 del freezing anotadas como fascina 1 y que sean no redundantes.

df <- FascinRSCA::freezing_cluster %>%
  filter(annotation == "fascin1") %>%
  distinct(fasta_seq, .keep_all = TRUE)
dim(df)

A continuación, vamos a escribir un nuevo cabecero para las secuencias que contenga información taxonómica en el formato de pySCA.

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:5]

El formato es el siguiente: fasta_header original (único y con '|' sustituido por '_')|nombre de la entrada| Nombre de la especie | Categorías taxonómicas a las que pertenece en orden jerárquico

A continuación, vamos a escribir dichas secuencias en un nuevo archivo fasta. Para ello usaremos la función multiple_fasta_files

fasta.file <- df %>%
  FascinRSCA::multiple_fasta_files(group_by_cols = "annotation", seq_col = "fasta_seq",
                       header_col = "annotated_header", path ="fascin1", toWrite = TRUE)

La distribución taxonómica de las secuencias que ahora componen el archivo fascin1.fasta es el siguiente:

edges <- df %>%
  select(phylum, class,order, family)%$%
  rbind.data.frame(list(phylum, class) %>% setNames(c("from", "to")),
                   list(class, order) %>% setNames(c("from", "to")),
                   list(order, family) %>% setNames(c("from", "to"))) %>%
  drop_na()
vertices <- df %$%
  rbind.data.frame(fct_count(phylum) %>% setNames(c("nodes", "values")) %>%
                      mutate(level = "phylum"),
                   fct_count(class) %>% setNames(c("nodes", "values"))%>%
                      mutate(level = "class"),
                   fct_count(order) %>% setNames(c("nodes", "values"))%>%
                      mutate(level = "order"),
                   fct_count(family) %>% setNames(c("nodes", "values"))%>%
                      mutate(level = "family"))%>%
  drop_na()
vertices$new_label <- with(vertices, ifelse(values > 4, as.character(nodes), NA))

mygraph <- graph_from_data_frame( edges, vertices=vertices )


ggraph(mygraph, layout = 'circlepack', weight=values) + 
  geom_node_circle(aes(fill = depth)) +
  geom_node_label( aes(label=new_label, fill=depth, size=values),
                   repel = TRUE, color = "white") +
  theme_void() + 
  theme(legend.position="FALSE") + 
  scale_fill_viridis()

T-Coffee

A continuación, vamos a emplear la herramienta t_coffee para realizar el alineamiento. En primer lugar, vamos a seleccionar las secuencias más diversas (el 30% más "informativo") pero forzando a mantener las secuencias en humano.

```{bash, eval = FALSE} t_coffee -other_pg seq_reformat -in fascin1.fasta -action +trim _seq_N30_fNAME HUMAN \ -output fasta_seq > fascin1.trim.fasta

Será necesario proteger el header de nuestras secuencias (al ser tan largo). Vamos a codificarlo:

```{bash, eval = FALSE}
t_coffee -other_pg seq_reformat -in fascin1.trim.fasta -output \
    code_name > fascin1.trim.code_name

t_coffee -other_pg seq_reformat -code fascin1.trim.code_name -in \   
    fascin1.trim.fasta > fascin1.trim.coded.fasta  

A continuación, usamos el servidor tcoffee para realizar un alineamiento con el archivo fascin1.trim.coded.fasta usando el modo accurate.

```{bash, eval = FALSE} t_coffee -in=fascin1.trim.coded.fasta -mode=accurate -blast=LOCAL -output=score_html clustalw_aln fasta_aln score_ascii phylip -maxnseq=150 -maxlen=2500 -case=upper -seqnos=off -outorder=input -run_name=result -multi_core=4 -pdb_db=/db/pdb/derived_data_format/blast/2021-11-16/pdb_seqres.fa -protein_db=/db/ncbi/201511/blast/db/nr.fa -quiet=stdout

Por último, decodificamos el resultado final:
```{bash, eval = FALSE}
t_coffee -other_pg seq_reformat -decode fascin1.trim.code_name -in \
    fascin1.trim.coded.fasta_aln > fascin1.trim.fasta_aln

HMMER

En primer lugar, vamos a construir el HMM:

```{bash, hmmbuild, eval = FALSE} hmmbuild fascin1.hmm fascin1.trim.fasta_aln > fascin1.desc

idx name nseq alen mlen eff_nseq re/pos description

---- -------------------- ----- ----- ----- -------- ------ -----------

1 fascin1.trim 26 611 492 0.78 0.588

CPU time: 0.41u 0.00s 00:00:00.41 Elapsed: 00:00:00.42

Como refleja la salida del programa, el el alineamiento `fascin1.trim` está compuesto por 26 secuencias y tiene 611 columnas alineadas. HMMER lo ha convertido en un perfil de 492 posiciones consenso, lo que significa que ha tomado 119 columnas del alineamiento que contenían gaps como inserciones relativas al consenso.

Las 26 secuencias fueron solo contadas como 0.78 secuencias efectivas porque son muy parecidas entre ellas. El perfil terminó con una entropía relativa por posición de 0,588 bits. 

A continuación, usamos `hmmsearch`. Vamos a buscar frente las siguientes bases de datos: SwissProt, RefProteomes y Ensembl. 

Los resultados son los siguientes:

```r
path <- "../inst/extdata/ann_fasta/annotation/hmm/search/"

files_hits <- purrr::map_chr(dir(path, pattern = "*.tsv"), ~paste0(path, .x))
hits.hmmer <- purrr::map_dfr(files_hits, ~.x %>%
  read_tsv() %>%
  mutate(db = .x %>%
           str_match("search/(.*)\\..*$") %>%
           extract(2)))
glimpse(hits.hmmer)

Para interpretar los resultados, vamos a fijarnos primero en el E-value (número esperado de falsos positivos). Escogeremos aquellas secuencias con un valor E-values $< 10^{−3}$. No obstante, para comparar resultados entre bases de datos distintas sería maś adecuado emplear el sequence bit score, puesto que no depende del tamaño de la base de datos.

Únicamente nos fijaremos en el valor de bias cuando sea del mismo orden de magnitud que sequence bit score. Los siguientes valores corresponden al domino que mejor ha puntuado y, en un principio, solo los tendremos en cuenta cuando el single best domain E-value $> 10^{−3}$, cuando no sean significativos, puesto que eso puede significar que sea un homólogo lejano multidominio.

hits.hmmer %>% 
  pull(db) %>% fct_count()
hits.hmmer  <- hits.hmmer %>%
  mutate(fullseq.Evalue = `E-value`)%>%
  group_by(`Target Name`)%>%
  mutate(best.Evalue = min(`Domain Ind. E-value`)) %>%
  ungroup()
hits.hmmer %>%
  select(db, fullseq.Evalue, best.Evalue, Score, Bias) %$%
  split(., db) %>%
  map(summary)

Como podemos observar el valor máximo de E-value de la secuencia completa para las tres búsquedas realizadas es inferior al valor establecido como umbral, $\text{E-value} < 10{-3}$. Esto no ocurre, sin embargo para el best.Evalue, el mejor E-value de los dominios individuales de la proteína. A continuación, vamos a filtrar los casos en que Domain Ind. E-value no es significativo. Estas secuencias, como se comentaba antes, pueden ser secuencias homólogas multidominio o una secuencia repetitiva (puesto que, matemáticamente, la hipótesis nula que se prueba es si se trata de una secuencia aleatoria, y una secuencia repetitiva no es aleatoria). Destacar que escogemos como E-value del dominio el Domain Ind. E-value (el independiente, y no condicional, que supone que la secuencia sí es homóloga).

dim(hits.hmmer)
hits.hmmer.filtered <- hits.hmmer %>%
  filter(fullseq.Evalue < 10 ^-3 & best.Evalue < 10 ^-3)
dim(hits.hmmer.filtered)
hits.hmmer.redflag <- hits.hmmer %>%
  filter(fullseq.Evalue < 10 ^-3 & best.Evalue > 10 ^-3)
dim(hits.hmmer.redflag)

Exploramos en profundidad las secuencias que hemos puesto en duda que sean homólogas (aquellas donde el $\text{E-value}$ más bajo de todos los dominios de la secuencia no es significativo, aunque la secuencia completa sí lo sea). Primero mostramos la descripción de dichas secuencias (mejor dicho, de los dominios de dichas secuencias).

En segundo lugar, mostramos un gráfico donde se muestra el $-\log_{10}(\text{E-value})$, tanto para la secuencia completa (coloreados en verde, a la derecha de la línea que indica el umbral de significancia escogido) como para cada uno de los dominios (coloreado en rojo, a la izquierda de la línea). En el eje de la Y están representadas cada una de las secuencias correspondientes a proteínas y la línea gris horizontal une el valor de $-\log_{10}(\text{E-value})$ para la secuencia completa con el del dominio con un mayor valor de signifcancia. Este gráfico permite visualizar

hits.hmmer.redflag %>%
  mutate(Description = fct_lump(Description,n = 10,other_level = "Other"))%>%
  pull(Description)%>%
  fct_count()

hits.hmmer.redflag %>%
  mutate(Description = fct_lump(Description,n = 3,other_level = "Other"))%>%
  arrange(Description)%>%
  ggplot() +
  geom_segment( aes(x=`Target Name`, xend=`Target Name`, y=-log(fullseq.Evalue),
    yend=-log(best.Evalue)), color="grey") +
  geom_point( aes(x=`Target Name`, y=-log(fullseq.Evalue)), color=rgb(0.2,0.7,0.1,0.5), size=2 ) +
  geom_point( aes(x=`Target Name`, y=-log(`Domain Ind. E-value`)), color=rgb(0.7,0.2,0.1,0.5), size=2 ) +
  coord_flip()+
  theme(
    panel.grid.major.x = element_blank(),
    panel.border = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    strip.text.y = element_text(angle = 0),    
    strip.text.x = element_text(angle = 0)
  ) +
  xlab("") +
  ylab("-log(E-value)")+
  geom_hline(yintercept = -log(10^-3), alpha = 0.5)+
  facet_grid(rows = vars(Description), scales = "free", space = "free",
             margins = FALSE,shrink = TRUE)

Por el momento, vamos a excluir estas secuencias, a la espera de profundizar en el significado de estos resultados (Revisión Coral).

A continuación, vamos a construir un DataFrame que contenga la información descargada de HMMER para estas secuencias.

files_fasta <- purrr::map_chr(dir(path, pattern = "*.fa"), ~paste0(path, .x))
fasta_seq <- Biostrings::readAAStringSet(files_fasta)[hits.hmmer.filtered$`Target Name`]

Y escribimos dicho archivo fasta para poder alinearlo a nuestro modelo:

Biostrings::writeXStringSet(fasta_seq, filepath = "hmmer_search.fullseq.fasta")

A continuación, vamos a alinearlo al HMM que habíamos construido y añadiremos, además, el alineamiento del cual se partió:

{bash, eval = FALSE} hmmalign --mapali fascin1.trim.fasta_aln fascin1.hmm hmmer_search.fullseq.fasta > hmmer_search.aln

currocam/FascinRSCA documentation built on March 21, 2022, 6:29 a.m.