Figures_vcmsa.Rmd

---
title: "vcMSA figures"
output: html_document
---



Figures and analysis accompanying
"Vector-clustering Multiple Sequence Alignment: Aligning into the twilight zone of protein sequence similarity with protein language models"
Claire McWhite, Isabel Armour-Garb, & Mona Singh


Set up
```{r}

library(tidyverse)
library(colorspace)
library(cowplot)
library(plotly)
library(ggridges)
library(fs)
theme_set(theme_cowplot(font_size = 8))

#1 = orange, 2 = lightblue, 3 = green, 4  = yellow, 5 = blue, 6 = red, 7 = pink, 8 = gray
palette_OkabeIto <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999")
colorspace::swatchplot(palette_OkabeIto)


```

```{r}
library(seqinr)
read_fasta <- function(fasta_filename, annot = FALSE){
    fasta <- seqinr::read.fasta(fasta_filename, as.string = TRUE)

    # Convert seqinr SeqFastadna object to data.frame
    fasta_df <- fasta %>%
                   sapply(function(x){x[1:length(x)]}) %>%
                   as.data.frame %>%
                   broom::fix_data_frame(newcol = "ID", newnames = "Sequence")

    if(annot == TRUE){
        annot_df <- getAnnot(fasta) %>%
                         sapply(function(x){x[1:length(x)]}) %>%
                         as.data.frame() %>%
                         broom::fix_data_frame(newnames = "Annot")

        fasta_df <- cbind(fasta_df, annot_df)
    }
    return(fasta_df)
}
```



Analysis for Figure 2A & B
```{r}

files <- fs::dir_ls("/scratch/gpfs/ia3026/vcmsa_tests/mnc_final/",glob = "*.csv")# %>% head(200)

dat <- read_csv(files, id="path",  show_col_types = FALSE, col_names = c("fasta", "seq1", "seq2", "aa1", "aa2", "sim"))



dat_formatted_fwd <- dat %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% filter(seq1 == 0 & seq2 == 1 )

dat_formatted_rev <- dat %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% filter(seq1 == 1 & seq2 == 0 ) %>% rename_at(vars(seq1, seq2), ~ if_else(. == "seq1", "seq2", "seq1")) %>% rename_at(vars(aa1, aa2), ~ if_else(. == "aa1", "aa2", "aa1"))



dat_formatted_new <- rbind(dat_formatted_fwd, dat_formatted_rev)

dat_formatted_new %>% full_join(dat_formatted, by = c("protein", "fasta", "seq1", "seq2", "aa1", "aa2")) %>%
  filter(protein == "6PF2K")

refs <- read_csv("/scratch/gpfs/cmcwhite/aln_datasets/quantest2/QuanTest2/Labels/match_positions_both.csv", col_names = c("protein", "aa1", "aa2", "label"))


total_aligned_positions <- refs %>%
  filter(label == 1) %>%
  filter(grepl("^0-", aa1)) %>%
  filter(grepl("^1-", aa2)) %>%
  group_by(protein) %>%
     summarize(aligned_positions = n())


dat_formatted_ref <- dat_formatted %>%
  left_join(refs)


dat_evaluated <- dat_formatted_ref %>% filter(label == 1)

dat_evaluated %>%
  group_by(label) %>%
     summarize(mean = mean(sim))

dat_evaluated %>%
  group_by(label) %>%
     summarize(n = n())

bc_eval <- dat_evaluated %>%
  filter(label ==1) %>%
  group_by(protein, label) %>%
  summarize(n = n()) 
```



Figure 2A and 2B
```{r}
theme_set(theme_cowplot(font_size = 8))
# Figure 2A
correct_pos<- bc_eval %>%
  left_join(total_aligned_positions) %>%
  mutate(prop_correct = n/aligned_positions) %>%
     ggplot(aes(x= prop_correct)) +
        geom_histogram(binwidth = 0.05) +
  scale_y_continuous(expand = c(0,0)) +
  xlim(0,1) +
 
  xlab("Proportion correct positions") +
  ylab("Count of proteins")
correct_pos


# In text value of number of alignments over where >50% of aligned columns are also reciprocal best hits
num_aligns <- nrow(bc_eval)
num_aligns_over_0.5correct <- bc_eval %>%
  left_join(total_aligned_positions) %>%
  mutate(prop_correct = n/aligned_positions) %>%
  filter(prop_correct >= 0.5) %>%
  nrow()


# 93.04
num_aligns_over_0.5correct/num_aligns


homstrad_sims <- read_csv("/scratch/gpfs/cmcwhite/aln_datasets/homstrad_aln_sims.csv") %>%
  mutate(protein = str_replace(file, ".filt.aln", "")) %>% select(-file)

# Figure 2B
sim_correct_plot <- bc_eval %>%
  left_join(total_aligned_positions) %>%
  mutate(prop_correct = n/aligned_positions) %>%
  left_join(homstrad_sims)%>%
     ggplot(aes(x= same, y =  prop_correct, label = protein)) +
        geom_point(color = "grey50", alpha = 0.6, size = 1) +
  ylab("Proportion correct positions") +
  xlab("Reference alignment identity") +
    scale_y_continuous(limits = c(0,1)) +
  background_grid(major = "y")
sim_correct_plot

Figure2 <- plot_grid(correct_pos, sim_correct_plot, labels= c("A", "B"), label_size = 8)
Figure2 %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/Figure2.png", ., height = 2, width = 4,, units = "in", dpi = 600)
Figure2 %>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/Figure2.pdf", ., height = 2, width = 4, units = "in")

```


Remove files where the three reference, and rest of the proteins are substantially different sizes
```{r}
directory_fasta <- "/scratch/gpfs/cmcwhite/aln_datasets/quantest2/QuanTest2/Test/"
fasta_files <- fs::dir_ls(directory_fasta, glob= "*20seqs*fasta")#
length_diffs <- map_dfr(fasta_files, ~read_fasta(paste0(.)), .id = "fasta") %>%
  mutate(prot = basename(fasta)) %>%
  mutate(prot = str_replace(prot, ".vie.20seqs.fasta", "")) %>%
  select(-fasta) %>%
  mutate(ref = case_when(ID %in% c("seq0001", "seq0002", "seq0003") ~ TRUE,
                         TRUE~ FALSE)
                ) %>%
  mutate(seqlen = str_length(Sequence)) %>%
  group_by(prot, ref) %>%
     summarize(mean_len = mean(seqlen)) %>%
  pivot_wider(names_from = ref, values_from = mean_len) %>%
  mutate(diff = `TRUE`/`FALSE`) %>%
  mutate(log2diff = log2(diff)) %>%
  mutate(abs_log2fc = abs(log2diff)) %>% arrange(desc(abs_log2fc))

extreme_length_diffs <- length_diffs %>% filter(abs_log2fc >= 2) %>% pull(prot)

```


Input files
```{r}
# Chainbreak seems to be less of a problem now
#chainbreak_prot <- read_csv("/home/cmcwhite/aln_results/families_to_exclude.txt", col_names = "prot")

# Gold standard identities
sims <- read_csv("/home/cmcwhite/aln_results/Sims_ref.csv") %>%
  mutate(prot = str_replace(file, ".ref", "")) %>% select(-file) %>% 
  # filter(!prot %in% chainbreak_prot$prot) %>%
    mutate(identbin = round(same, digits = 1)) 

allscores <- read_csv("/scratch/gpfs/cmcwhite/aln_comparison/eval_scores.csv", col_names = FALSE) %>%
  select(score = X1, experiment = X2, score_type = X4, prot = X5) %>%
  filter(!is.na(prot)) #%>%
  #mutate(set = 20) 


allscores_standard <- allscores %>%
  mutate(score_type = str_replace(score_type, "_", "")) %>%
  filter(score_type %in% c("tc", "tcHGI", "tcHGIEB")) %>%
  separate(experiment, into = c("tmp0", "tmp1", 'set', "tmp2", 'layers', 'model', 'padding', 'threshold', 'bc', 'rbh', 'tmp3'), sep = "_", extra = "merge", remove = FALSE) %>%
  separate(prot , into = c("prot", "style", "tmp4", "alg"), extra = "drop", sep = "[.]") %>%
  mutate(layers = case_when(is.na(layers) ~ "None", TRUE ~ layers)) %>%
  mutate(model = case_when(is.na(model) ~ "None", TRUE ~ model))  %>%
  mutate(alg = case_when(alg == "with" ~ tmp4, TRUE ~ alg)) %>%
  mutate(alg = case_when(grepl("layer", experiment) ~ "semantic", TRUE ~ alg))  %>%
  filter(alg != "MAFFT-FINSI") %>% # Not named correctly, this is MAFFT-LINSI
  #mutate(set = "20") %>%
  
  mutate(score= as.numeric(score))# %>%



```


Secondary structure measure 
```{r}

 allscores_ss <- allscores %>%
  filter(!score_type %in% c("sp", "tc", "col")) %>% 
  separate(experiment, into = c('set', NA, 'layers', NA, 'model'), sep = "_", extra = "merge", remove = FALSE) %>%
  separate(prot , into = c("prot",  NA, "alg"), extra = "drop", sep = "[.]") %>%
  mutate(layers = case_when(is.na(layers) ~ "None", TRUE ~ layers)) %>%
  mutate(model = case_when(is.na(model) ~ "None", TRUE ~ model))  %>%
  mutate(alg = case_when(grepl("layer", experiment) ~ "semantic", TRUE ~ alg))  %>%
  filter(experiment != "20seqs_nopca") %>%
  mutate(score= as.numeric(score))# %>%
       #filter(!prot %in% chainbreak_prot$prot)



 #allscores_ss %>%
#    filter(model != "fa_hilo") %>%
#    filter(model != "fa_hilo_tmp") %>%
#  filter(score_type == "col_HGIEB") %>%
#   unique() %>%
#  group_by(experiment, model, alg) %>%
#     summarize(medianscore = median(score), meanscore = mean(score),n = n()) %>%
#  ungroup %>%
#  arrange(desc(medianscore)) %>% View()

```




### Functions for comparing and plotting
```{r}

get_layer_compare <- function(allscores_standard, experiment_name, score_type = "tc", ref = "20seqs_prog") {

layer_compare_sel <- allscores_standard %>%
    filter(experiment %in% c({{experiment_name}}, {{ref}})) %>% # 
  left_join(sims) %>%
  unique %>%
  select(score, prot, alg, score_type, same, sim, identbin)


  print(layer_compare_sel)


layer_compare_semantic <- layer_compare_sel %>% filter(alg == "semantic") %>% select(-alg) %>% rename(semantic = score)
layer_compare_others <- layer_compare_sel %>% filter(alg != "semantic") %>% rename(other = score)
layer_compare_compare <- layer_compare_semantic %>%
  filter(score_type == {{score_type}}) %>%

  left_join(layer_compare_others, by = c("prot", "score_type", "same", "sim", "identbin")) %>%
    filter(!is.na(alg))

print(layer_compare_sel)
prots_per_identbin <- layer_compare_compare %>%
  select(prot, identbin) %>%
  unique() %>%
  group_by(identbin) %>%
   summarize (n= n()) %>%
  ungroup()

perform_plot_setup <- layer_compare_compare %>% mutate(scorediff = semantic - other) %>%
  mutate(winner = case_when(scorediff > 0 ~ "Better",
                            scorediff < 0 ~ "Worse")) %>%
  filter(!is.na(winner)) %>%

  group_by(alg, score_type, winner) %>%
    mutate(totwins = n()) %>%
  ungroup %>%
  group_by(alg, score_type, winner, identbin, totwins) %>%
      summarize(num = n()) %>%
  ungroup()

return(layer_compare_compare)
}

# Figure 3A plotting function
perform_plot_fxn <- function(layer_compare_compare, order2){


print(layer_compare_compare )


perform_plot <-layer_compare_compare %>%
  mutate(alg = paste0(alg," vs. vcMSA" )) %>%
  mutate(win = case_when(vcMSA > other ~ "win", TRUE ~ "loss")) %>%
  arrange(win) %>%
  mutate(win = fct_inorder(win)) %>%
  mutate(new_y = (other + vcMSA )/sqrt(2)) %>%
    mutate(new_x = (vcMSA -other)/sqrt(2)) %>%
  mutate(diff =vcMSA  - other) %>%
  
  ggplot(aes(x = new_x, y = new_y)) +
  geom_point(aes( label = paste(prot, other, vcMSA ), color = same), alpha = 0.7, size = 1) +
 
  geom_abline(intercept = 0, slope = 1, size = 0.1, color = "black") + 
  geom_abline(intercept = 0, slope =-1, size = 0.1, color = "black") + 
      theme(strip.background = element_blank()) +
  xlab("Column score") +
  ylab("vcMSA\nColumn score") + 
   geom_vline(xintercept = 0, color = "red", linetype = "dashed") + 
  scale_color_continuous_sequential(palette = "viridis", name = "Sequence\nidentity", rev = FALSE) +
  scale_y_continuous(limits = c(-.1, 150.1), expand = c(0,0)) +
    scale_x_continuous(limits = c(-70.7,70.7), expand = c(0,0)) +
 facet_wrap(~fct_relevel(alg, order2)) +
  theme(panel.spacing.x =  unit(1.5, "lines")) +
  theme(axis.line = element_blank(), axis.ticks = element_blank(), axis.title = element_blank(), axis.text = element_blank())

return(perform_plot)
}

```

### Figure 3A


delete
```{r}

theme_set(theme_cowplot(font_size = 8))


#score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis<- allscores_standard %>%
#      filter(!prot %in% extreme_length_diffs) %>%
#  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis", score_type = "tc") %>%
#  rename(vcMSA = semantic)
  
#a <- perform_plot_fxn(score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis)


order2 <- c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,  'FAMSA vs. vcMSA','PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "MAFFT-LINSI vs. vcMSA", "UPP vs. vcMSA", "TCOFFEE vs. vcMSA")

#score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis_hfa<- allscores_standard %>%

a <- perform_plot_fxn(score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust , order2)
a


score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis_hfa<- allscores_standard %>%
      filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis_hfa", score_type = "tc") %>%
  rename(vcMSA = semantic)
  
b <- perform_plot_fxn(score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis_hfa, order2)




plot_16_nopad


score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis<- allscores_standard %>%
      filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis", score_type = "tc") %>%
  rename(vcMSA = semantic)
  
b<- perform_plot_fxn(score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis)
plot_16_nopad


score_compare_aln_evaluate_20_nopca_16layer_t5_nopad<- allscores_standard %>%
      filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding0_nobc", score_type = "tc") %>%
  rename(vcMSA = semantic)
  
```


delete
```{r}


  #%>% pull(prot) %>% unique


bind_rows(score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust, add_extra)


order2 <- c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,  'FAMSA vs. vcMSA','PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "MAFFT-LINSI vs. vcMSA", "UPP vs. vcMSA", "TCOFFEE vs. vcMSA")
plot_16_nopad <- perform_plot_fxn(bind_rows(score_compare_aln_evaluate_20_nopca_16layer_t5_nopad, add_extra), order2)
plot_16_nopad


plot_16_nopad %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_nopad.png", ., height = 5, width = 5.8,, units = "in", dpi = 600)
plot_16_nopad%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_nopad.pdf", ., height = 5, width = 5.8, units = "in")


```


```{r}
score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_zcbc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_zcbc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_none_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_none_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)


#scoring
score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_thr70_bc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding0_thr70_bc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)



score_compare_aln_evaluate_50_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_50_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust", score_type = "tc", ref = "50seqs_prog") %>%
  rename(vcMSA = semantic)




```



```{r}

add_extra <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust %>% filter(alg == "CLUSTALO") %>%
    mutate(alg = "extra")

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

order3 <-c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,'PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA",  'FAMSA vs. vcMSA', "UPP vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "TCOFFEE vs. vcMSA")

plot_16_mainref <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust %>% filter(alg != "MAFFT-LINSI") %>%
                                                                   bind_rows(add_extra) %>%
  perform_plot_fxn(order3)


plot_16_mainref

plot_grid(plot_16_main, plot_16_mainref)

plot_16_main %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.png", ., height = 5, width = 5.8,, units = "in", dpi = 600)
plot_16_main%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.pdf", ., height = 5, width = 5.8, units = "in")


```


```{r}

add_extra <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust %>% filter(alg == "CLUSTALO") %>%
    mutate(alg = "extra")

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

order3 <-c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,'PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA",  'FAMSA vs. vcMSA', "UPP vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "TCOFFEE vs. vcMSA")

plot_16_main <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust %>% filter(alg != "MAFFT-LINSI") %>%
                                                                   bind_rows(add_extra) %>%
  perform_plot_fxn(order3)


plot_16_main
plot_16_main %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.png", ., height = 5, width = 5.8,, units = "in", dpi = 600)
plot_16_main%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.pdf", ., height = 5, width = 5.8, units = "in")


```


```{r}

add_extra <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust %>% filter(alg == "CLUSTALO") %>%
    mutate(alg = "extra")

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

order3 <-c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,'PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA",  'FAMSA vs. vcMSA', "UPP vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "TCOFFEE vs. vcMSA")

plot_16_main_ref <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust %>% filter(alg != "MAFFT-LINSI") %>%
                                                                   bind_rows(add_extra) %>%
  perform_plot_fxn(order3)


plot_grid(plot_16_main, plot_16_main_ref , plot_16_main_50, nrow = 1)
plot_16_main %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.png", ., height = 5, width = 5.8,, units = "in", dpi = 600)
plot_16_main%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main.pdf", ., height = 5, width = 5.8, units = "in")


```


```{r}

score_compare_aln_evaluate_50_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust <- allscores_standard %>%
  get_layer_compare(., "aln_evaluate_50_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)


add_extra <- score_compare_aln_evaluate_50_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust %>% filter(alg == "CLUSTALO") %>%
    mutate(alg = "extra")


order3 <-c("extra vs. vcMSA", "MAFFT-FFTNS1 vs. vcMSA" ,'PROBCONS vs. vcMSA', "MUSCLE vs. vcMSA",  'FAMSA vs. vcMSA', "UPP vs. vcMSA", "CLUSTALO vs. vcMSA", "MAFFT-GINSI vs. vcMSA", "TCOFFEE vs. vcMSA")

plot_16_main_50 <- score_compare_aln_evaluate_50_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust %>% filter(alg != "MAFFT-LINSI") %>%
                                                                   bind_rows(add_extra) %>%
  perform_plot_fxn(order3)


plot_16_main_50
plot_16_main_50 %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main_50.png", ., height = 5, width = 5.8,, units = "in", dpi = 600)
plot_16_main_50%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/plot_16_main_50.pdf", ., height = 5, width = 5.8, units = "in")


```



```{r}
score_compare_aln_evaluate_20_nopca_16layer_t5_nopad_HGIEB <- allscores_standard %>%
      filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding0_nobc", score_type = "tc_HGIEB") %>%
  rename(vcMSA = semantic) 
perform_plot_fxn(score_compare_aln_evaluate_20_nopca_16layer_t5_nopad_HGIEB)

score_compare_aln_evaluate_20_nopca_16layer_t5_nopad_HGIEB
score_compare_aln_evaluate_20_nopca_16layer_t5_nopad_HGIEB


```

In text values about combining subalignments
```{r}
has_excluded <- read_csv("/scratch/gpfs/cmcwhite/aln_comparison/aln_evaluate_20_nopca_16layer_t5_padding0_nobc//alignment_files/has_excluded.txt") %>% mutate(set = "excluded")
has_clusters <- read_csv("/scratch/gpfs/cmcwhite/aln_comparison/aln_evaluate_20_nopca_16layer_t5_padding0_nobc//alignment_files/has_seqgroups.txt") %>% mutate(set =  "seqgroups") 

both_stitched <- bind_rows(has_clusters, has_excluded) %>%
  separate(filename, into = c(NA, "filename"), sep= "[/]") %>%
  separate(filename, into = c("prot", NA), sep = "[.]") %>%
  filter(!prot %in% extreme_length_diffs) 


both_stitched %>% group_by(set) %>%
     summarize(n = n())

both_stitched %>% select(prot) %>% unique %>% nrow()



both_stitched %>%
  left_join(sims)


```


score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust <- allscores_standard %>%

     filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)

#### Values in Figure 3A
Win percent
```{r}
score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust%>%
           #filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))


score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust%>%
           filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))


score_compare_aln_evaluate_50_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust%>%
           #filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))


```



delete
```{r}


score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis_hfa %>%
           #filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis_hfa %>%
           filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))



score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis %>%
           filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis %>%
           filter(!prot %in% extreme_length_diffs) %>%
  mutate(win = case_when(vcMSA >= other ~ "Win",
                         TRUE ~ "Lose")) %>%
  group_by(score_type, alg, win) %>%
      summarize(n = n()) %>%
  ungroup %>% 
  pivot_wider(names_from = win, values_from = n) %>%
  arrange(desc(Win))


```

```{r}



theme_set(theme_cowplot(font_size = 8))

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust
score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_thr70_bc_mnc_oldclust

pad_compare <- score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_thr70_bc_mnc_oldclust %>% 
    filter(!prot %in% extreme_length_diffs) %>%
  left_join(score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust, by= c("prot" ,"score_type", "other", "identbin", "same", "sim", "alg")) %>%
  mutate(color = case_when(vcMSA.x >= vcMSA.y ~ "win", TRUE ~ "lose")) %>%
  select(-alg, -other) %>%
  unique


padding_plot <- pad_compare %>%  ggplot(aes(x = vcMSA.x, y = vcMSA.y, label = prot, color = color)) +
  geom_abline() + 
  geom_point(size= 0.8, alpha = 0.8)  +
  xlab("No padding") + 
  ylab("Padding (10X's)") + 
  scale_color_manual(values = palette_OkabeIto) +
  theme(legend.position = "none")

padding_plot

padding_plot %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/padding_plot.png", ., height = 1.4, width = 1.4, units = "in", dpi = 600)
padding_plot%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/padding_plot.pdf", ., height = 1.4, width = 1.4, units = "in")


# On plot value
pad_compare %>%
  group_by(color) %>%
     summarize(n= n())
```


```{r}



theme_set(theme_cowplot(font_size = 8))

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_nobc_mnc_oldclust <- allscores_standard %>%

     #filter(!prot %in% extreme_length_diffs) %>%
  get_layer_compare(., "compare_aln_evaluate_20_nopca_16layer_t5_padding10_thr70_nobc_mnc_oldclust", score_type = "tc") %>%
  rename(vcMSA = semantic)


nobc <- allscores_standard %>% filter(experiment %in% c("aln_evaluate_20_nopca_16layer_t5_padding10_thr70_nobc_mnc_oldclust")) %>% select(prot, score_type,score)
bc <- allscores_standard %>% filter(experiment %in% c("aln_evaluate_20_nopca_16layer_t5_padding10_thr70_refbc_mnc_oldclust")) %>% select(prot, score_type, score) 

batch_compare <- nobc %>% 
  left_join(bc, by= c("prot" ,"score_type")) %>%
    filter(!prot %in% extreme_length_diffs) %>%
  mutate(color = case_when(score.x >= score.y ~ "win", TRUE ~ "lose")) %>%
  filter(score_type == "tc") %>%
  unique



#batch_compare <- nobc %>% full_join(bc, by = c("prot", "score_type")) %>% filter(score_type == "tc")  #%>%
  #ggplot(aes(x = score.x, y = score.y)) +
  #geom_point() +
  #geom_abline(yintercept = 1, slope = 1)




batching_plot <- batch_compare %>%  
  filter(!prot %in% extreme_length_diffs) %>%
  ggplot(aes(x = score.x, y = score.y, label = prot, color = color)) +
  geom_abline() + 
  geom_point(size= 0.8, alpha = 0.8)  +
  xlab("No batch correction\nTotal column score") + 
  ylab("Batch correction\nTotal column score") + 
  scale_color_manual(values = palette_OkabeIto) +
  theme(legend.position = "none")

batching_plot

batching_plot %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/batching_plot.png", ., height = 2, width = 2, units = "in", dpi = 600)
batching_plot%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/batching_plot.pdf", ., height = 2, width = 2, units = "in")


# On plot value
batch_compare %>%
  group_by(color) %>%
     summarize(n= n())
```



```{r}
pad_batch <- plot_grid(padding_plot, batching_plot, ncol = 1, align = "hv")

pad_batch %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/pad_batch.png", ., height = 3, width = 1.6, units = "in", dpi = 600)
pad_batch%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/pad_batch.pdf", ., height = 3, width = 1.6, units = "in")

```

```{r}

score_compare_aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis %>% 
  left_join(score_compare_aln_evaluate_20_nopca_16layer_t5_nopad , by= c("prot" ,"score_type")) %>%
  ggplot(aes(x = vcMSA.x, y = vcMSA.y, label = prot)) + geom_point() + geom_abline()


score_compare_aln_evaluate_20_nopca_16layer_t5_padding0_nobc_mis %>% 
  left_join(score_compare_aln_evaluate_20_nopca_16layer_t5_nopad , by= c("prot" ,"score_type")) %>%
  ggplot(aes(x = vcMSA.x, y = vcMSA.y, label = prot)) + geom_point() + geom_abline()



```


```{r}

score_compare_aln_evaluate_20_nopca_16layer_t5_nopad %>%
           filter(!prot %in% extreme_length_diffs) %>%
   mutate(dif = vcMSA  - other) %>%
  group_by(alg) %>%
     summarize(totdiff = sum(dif)) %>%
   arrange(desc(totdiff))

```


```{r}

score_compare_aln_evaluate_20_nopca_16layer_t5_nopad %>%
  group_by(identbin, alg, score_type) %>%
     summarize(avg_score = median(vcMSA ), avg_other = median(other))


```


### Figure 3B
```{r}
clustalo <- allscores_standard %>%
        # filter(!prot %in% extreme_length_diffs) %>%
  filter(alg == "CLUSTALO") %>%
  filter(set ==20) %>%
  filter(score_type == "tc")

mg <- allscores_standard %>%
        # filter(!prot %in% extreme_length_diffs) %>%
  filter(alg == "MAFFT-GINSI") %>%
  filter(set ==20) %>%
  filter(score_type == "tc")

semantic <- allscores_standard %>%
        # filter(!prot %in% extreme_length_diffs) %>%
  filter(experiment == "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust") %>%
  filter(score_type == "tc") %>%
  mutate(alg = "vcMSA")

both <- bind_rows(mg, clustalo, semantic)  %>% left_join(sims)


ident_plot <- both %>% group_by(alg, identbin) %>%
  mutate(identbin = case_when(identbin >= 0.6 ~ 0.6,
                              TRUE ~ identbin)) %>%
    summarize(medianscore= median(score), n = n()) %>%
  ggplot(aes(x = identbin, y = medianscore, color = fct_relevel(alg, c("vcMSA",  "MAFFT-GINSI" ,"CLUSTALO")))) +
       geom_point(alpha = 1, size = 1.5) +
  
  scale_color_manual(values = c("blue", "#E69F00", "#D55E00") ,name = "Algorithm") +
  background_grid(major = "y") +
  xlab("Sequence identity bin") +
  ylab("Median Total Column score") +
  scale_x_continuous(breaks = c(-0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6), limits = c(0,0.6)) +
    scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100), limits = c(0,103), expand = c(0,0))
ident_plot

ident_plot %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/ident_plot.png", ., height = 2, width = 3,, units = "in", dpi = 600)
ident_plot%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/ident_plot.pdf", ., height = 2, width = 3, units = "in")

```


### Figure 3C
```{r}

#   separate(prot , into = c("prot", "style", NA, "alg"), extra = "drop", sep = "[.]")

sims2 <- sims %>%
  mutate(identbin2 = case_when(identbin == 0.1 ~ 0.1,
                               identbin == 0.2 ~ 0.2,
                               TRUE ~ 0.3))


allseq_scores <- allscores_standard %>%
  filter(!prot %in% extreme_length_diffs) %>%
    filter(score_type == "tc") %>% 
    filter(experiment == "aln_evaluate_20_nopca_16layer_t5_padding10_thr70_bc_mnc_oldclust" | experiment ==  "20seqs_prog") %>%
 left_join(sims2) %>%
  group_by(experiment, model, alg, identbin2) %>%
     summarize(medianscore = median(score), meanscore = mean(score),n = n() ) %>% #, medianlen = median(alnLen), meanlen= mean(alnLen)) %>%
  ungroup %>%

  arrange(desc(medianscore)) 
  
# RERUN!!!!!

struc_scores <- allscores_standard %>%
  #filter(model != "fa_hilo") %>%
  #  filter(model != "fa_hilo_tmp") %>%
  left_join(sims2)%>%
  #filter(same < 0.3) %>%
  filter(!prot %in% extreme_length_diffs) %>%
  filter(score_type == "tcHGIEB") %>% 
   # filter(prot %in% voro$prot) %>%
 # filter(alg %in% c( "semantic") | tmp4 == "TCOFFEE") %>%
  filter(experiment == "aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis" | experiment ==  "20seqs_prog") %>%
  
  group_by(experiment, model, alg, tmp4, identbin2) %>%
     summarize(medianscore_struct = median(score), meanscore_struct= mean(score),n = n() ) %>% #, medianlen = median(alnLen), meanlen= mean(alnLen)) %>%
  ungroup %>%

  arrange(desc(medianscore_struct))  %>%
  mutate(alg = case_when(alg != "semantic" ~ tmp4, TRUE ~ alg)) %>%
  mutate(order  = case_when(identbin2 == 0.1 ~ meanscore_struct, TRUE ~ 100))

theme_set(theme_cowplot(font_size = 8))


plot_ss <- allseq_scores %>%
left_join(struc_scores, by = c("alg", "identbin2")) %>%
  filter(alg != "MAFFT-LINSI") %>%
   # filter(alg %in% c("TCOFFEE", "CLUSTALO", "semantic")) %>%
  mutate(alg = case_when(alg == "semantic"~"vcMSA",
                         TRUE ~ alg)) %>%
  mutate(alg = fct_reorder(alg, order)) %>%
  ggplot(aes(x = meanscore, y= meanscore_struct, color = fct_rev(fct_reorder(alg, order)))) +
  geom_point() +
  xlim(0,100) +
  scale_color_manual(values = c("blue",palette_OkabeIto), name = "Algorithm") +
  theme(strip.background = element_blank()) +
  ylim(0,100) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  facet_wrap(~identbin2) +
  xlab("Overall median Total Score") +
  ylab("Secondary structure\nmedian Total Score")



```


```{r}

fig3bc <- plot_grid(ident_plot, plot_ss, nrow = 1, rel_widths = c(0.43,0.7))


fig3bc %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/fig3bc.png", ., height = 1.5, width = 6, units = "in", dpi = 600)
fig3bc%>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/fig3bc.pdf", ., height = 1.5, width = 6, units = "in")


```



```{r}
struc_scores_HGIEB <- allscores_standard %>%

  left_join(sims2)%>%
  filter(!prot %in% extreme_length_diffs) %>%
  filter(score_type == "tcHGIEB") %>% 
   # filter(prot %in% voro$prot) %>%
 # filter(alg %in% c( "semantic") | tmp4 == "TCOFFEE") %>%
  filter(experiment == "aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis" | experiment ==  "20seqs_prog") %>%
  
  group_by(experiment, model, alg, tmp4, identbin2) %>%
     summarize(medianscore_struct = median(score), meanscore_struct= mean(score),n = n() ) %>% #, medianlen = median(alnLen), meanlen= mean(alnLen)) %>%
  ungroup %>%

  arrange(desc(medianscore_struct))  %>%
  mutate(alg = case_when(alg != "semantic" ~ tmp4, TRUE ~ alg)) %>%
  mutate(order  = case_when(identbin2 == 0.1 ~ meanscore_struct, TRUE ~ 100))


struc_scores_TSC <- allscores_standard %>%

  left_join(sims2)%>%
  filter(!prot %in% extreme_length_diffs) %>%
  filter(score_type == "tcTSC") %>% 
   # filter(prot %in% voro$prot) %>%
 # filter(alg %in% c( "semantic") | tmp4 == "TCOFFEE") %>%
  filter(experiment == "aln_evaluate_20_nopca_16layer_t5_padding10_nobc_mis" | experiment ==  "20seqs_prog") %>%
  
  group_by(experiment, model, alg, tmp4, identbin2) %>%
     summarize(medianscore_struct = median(score), meanscore_struct= mean(score),n = n() ) %>% #, medianlen = median(alnLen), meanlen= mean(alnLen)) %>%
  ungroup %>%

  arrange(desc(medianscore_struct))  %>%
  mutate(alg = case_when(alg != "semantic" ~ tmp4, TRUE ~ alg)) %>%
  mutate(order  = case_when(identbin2 == 0.1 ~ meanscore_struct, TRUE ~ 100))

theme_set(theme_cowplot(font_size = 8))


plot_ss <- allseq_scores %>%
left_join(struc_scores, by = c("alg", "identbin2")) %>%
  filter(alg != "MAFFT-LINSI") %>%
   # filter(alg %in% c("TCOFFEE", "CLUSTALO", "semantic")) %>%
  mutate(alg = case_when(alg == "semantic"~"vcMSA",
                         TRUE ~ alg)) %>%
  mutate(alg = fct_reorder(alg, order)) %>%
  ggplot(aes(x = meanscore, y= meanscore_struct, color = fct_rev(fct_reorder(alg, order)))) +
  geom_point() +
  xlim(0,100) +
  scale_color_manual(values = c("blue",palette_OkabeIto), name = "Algorithm") +
  theme(strip.background = element_blank()) +
  ylim(0,100) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  facet_wrap(~identbin2) +
  xlab("Overall median Total Score") +
  ylab("Secondary structure\nmedian Total Score")



```




Supplemental Figure 1A
Illustrate first character effect
```{r}

dat_evaluated_rbh  %>% write_csv("/scratch/gpfs/cmcwhite/aln_comparison/dat_evaluated_rbh.csv")
dat_formatted_ref   %>% write_csv("/scratch/gpfs/cmcwhite/aln_comparison/dat_formatted_ref.csv")


# Conclusion from this... padding only really helps with the alignment is offset ex. 
# Where padding helps
#CODH (0-0-E/1-22-Q), cyclotide (0-0-S/1-3-S), cox (0-0-G/1-15-R), cytochr_c552 (0-0-A/1-6-G)
dat_firstpos <- dat_evaluated_rbh %>% filter(label ==1) %>% filter(grepl("0-0-", aa1))#

dat_firstpos %>%
  group_by(padding ,true_rbh) %>% 
    summarize(n = n())

dat_formatted_ref %>%
  group_by(label, model, padding, bc) %>%
     summarize(mean = mean(sim))


dat_sum <- dat_evaluated_rbh %>%
  filter(label ==1) %>%
  group_by(protein, true_rbh, label,model, padding, bc) %>%
     summarize(n = n()) %>%
  pivot_wider(names_from = true_rbh, values_from = n)

dat_sum

dat_sum %>%
  select(-`NA`) %>%
  pivot_wider(names_from = padding, values_from = true_rbh) %>%
  mutate(diff = padding10-padding0) %>% arrange(desc(diff)) %>% View()

library(ggridges)


dat_evaluated %>%
  mutate(label = replace_na(label, 0)) %>%
  ggplot(aes(x = sim, fill =  paste(layers, padding, bc, label),  y = paste(layers, padding, bc))) +
  geom_density_ridges(alpha =0.5) 

```



Supplemental Figure 1B
Illustrate batch correction
```{r}
                                     #alignment_files_ghf34_3.fasta.thresh0.70.13.t5.bc/ghf34_3.fasta.thresh0.70.13.t5.bc.alignment_group2.prebatch.png.pca.csv
#prebatch <- read_csv("/home/cmcwhite/alignment_files_ghf34_3.fasta.thresh0.70.13.t5.bc/ghf34_3.fasta.thresh0.70.13.t5.bc.alignment_group2.prebatch.png.pca.csv") %>% mutate(set ="Before\nbatch correction")
#postbatch <- read_csv("/home/cmcwhite/alignment_files_ghf34_3.fasta.thresh0.70.13.t5.bc/ghf34_3.fasta.thresh0.70.13.t5.bc.alignment_group2.postbatch.png.pca.csv") %>% mutate(set ="After\nbatch correction")


# This is using the reference batch correction
prebatch <- read_csv("/home/cmcwhite/alignment_files_ghf34_3.fasta.vcmsa.thsh0.70.padding10.16.t5half.pycomabtbc_ref.mnc.4.aln/ghf34_3.fasta.vcmsa.thsh0.70.padding10.16.t5half.pycomabtbc_ref.mnc.4.aln.alignment_group0.prebatch.png.pca.csv") %>% mutate(set ="Before\nbatch correction")
postbatch <- read_csv("/home/cmcwhite/alignment_files_ghf34_3.fasta.vcmsa.thsh0.70.padding10.16.t5half.pycomabtbc_ref.mnc.4.aln/ghf34_3.fasta.vcmsa.thsh0.70.padding10.16.t5half.pycomabtbc_ref.mnc.4.aln.alignment_group0.postbatch.png.pca.csv") %>% mutate(set ="After\nbatch correction")



prepostbatch <- bind_rows(prebatch, postbatch)

bc_plot <- prepostbatch %>%
  mutate(Sequence = case_when(seq == 0 ~ "Protein #1",
                           seq == 1 ~ "Protein #2",
                           seq == 2 ~ "Protein #3")) %>%
  ggplot(aes(x = `0`, y = `1`, color = Sequence)) +
     geom_point(size = 0.3, alpha= 0.7) +
  scale_color_manual(values= palette_OkabeIto, name = "") +
  facet_wrap(~fct_rev(set), nrow = 1, scales = "free") +
  theme(axis.ticks = element_blank(), axis.text = element_blank()) +
  xlab("PC1") +
  ylab("PC2") +
  theme(legend.position = "bottom") +
  theme(strip.background = element_blank())
bc_plot

bc_plot %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/bc_plot.png", ., height = 1.5, width = 2, units = "in", dpi = 600)
bc_plot %>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/bc_plot.pdf", ., height = 1.5, width = 2, units = "in")

```


Supplemental timings figure
```{r}

metrics <- read_csv("/scratch/gpfs/cmcwhite/aln_comparison//eval_metrics.csv", col_names = FALSE) %>%
  select(score = X1, experiment = X2, alg = X7, prot = X4, eval = X11) %>%
  filter(!is.na(prot)) %>%
  separate(experiment, into = c('set', NA, 'numseqs', NA, 'model'), sep = "_", extra = "merge", remove = FALSE) %>% unique


metrics_prog <- metrics %>% filter(grepl("prog", experiment)) %>%
  # filter(alg == "CLUSTALO") %>%
  pivot_wider(names_from = "eval", values_from = "score") %>%
  filter(set == "100seqs") #%>%
  ggplot(aes(x = set, y = (realtime/1000)/60, group = prot)) +
  geom_line(alpha = 0.1) + 
  geom_point(alpha = 0.1) +
  facet_wrap(~alg)

metrics_sel_nobc <- metrics %>% # filter(prot %in% c("AhpC-TSA", "bromodomain",  "lacI", "Mur_ligase_C", "sh2")) %>%
  filter(alg == "semantic") %>%
  filter(grepl("padding10_thr90_nobc_mnc_oldclust", experiment)) %>%
    filter(!grepl("prev", experiment)) %>%
  pivot_wider(names_from = "eval", values_from = "score")

metrics_sel_refbc <- metrics %>% # filter(prot %in% c("AhpC-TSA", "bromodomain",  "lacI", "Mur_ligase_C", "sh2")) %>%
  filter(alg == "semantic") %>%
  filter(grepl("padding10_thr70_refbc_mnc_oldclust", experiment)) %>%
  filter(!grepl("prev", experiment)) %>%
  pivot_wider(names_from = "eval", values_from = "score")



metrics_200_ref <- metrics_sel_refbc %>% filter(grepl("200", experiment)) %>% sample_n(100)
metrics_200_nobc <- metrics_sel_nobc %>% filter(grepl("200", experiment))# %>% sample_n(50)

metrics_200_both <- metrics_200_ref %>% inner_join(metrics_200_nobc, by = c("prot"))

metrics_200_both
bind_rows(metrics_200_nobc, metrics_200_ref) %>%
 group_by(prot) %>%
     summarize(n = n()) %>%
  filter(n==2)



# In text values
metrics_sel_refbc %>% group_by(experiment) %>%
  
    filter(prot %in% metrics_200_ref$prot) %>%
    summarize(t = median((realtime/1000)/60), min = min((realtime/1000)/60), max = max((realtime/1000)/60))





library(cowplot)
theme_set(theme_cowplot(font_size = 8))
timing_plot1 <- metrics_sel_refbc %>% #select(numseqs, eval, score, prot) %>%
 filter(prot %in% metrics_200_ref$prot) %>%
 #  bind_rows(traces_sel) %>%
   mutate(numseqs = as.numeric(numseqs)) %>%
 # mutate(numseqs = fct_relevel(numseqs, c("20", "50", "100", "200"))) %>%
  ggplot(aes(x = numseqs, y = (realtime/1000)/60, group = prot)) +
  ylab("Minutes") +
  geom_point(alpha=  0.1) +
  geom_line(alpha = 0.1) +
  # scale_y_log10(labels = scales::comma) +
   scale_x_continuous(labels = c(20, 50, 100, 200), breaks = c(20,50,100,200))+ 
   xlab("Number of sequences")

timing_plot1



  
timing_plot1 %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/timing_plot.png", ., height = 2.5, width = 3, units = "in", dpi = 600)
timing_plot1 %>%  ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/timing_plot.pdf", ., height = 2.5, width = 3, units = "in")


metrics_prog 


```



Supplemental true/false positives figure
Demonstrates benefit of maximum noncrossing matching
```{r}
#"mncm" refers to maximum noncrossing matching and "none" refers to no maximum noncrossing matching
files_none <- fs::dir_ls("/scratch/gpfs/ia3026/vcmsa_tests/none_final/",glob = "*.csv")# %>% head(200)

dat_none <- read_csv(files_none, id="path",  show_col_types = FALSE, col_names = c("fasta", "seq1", "seq2", "aa1", "aa2", "sim"))



dat_formatted_fwd_none <- dat_none %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% 
  filter(seq1 == 0 & seq2 == 1 )

dat_formatted_rev_none <- dat_none %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% filter(seq1 == 1 & seq2 == 0 ) %>% 
  rename_at(vars(seq1, seq2), ~ if_else(. == "seq1", "seq2", "seq1")) %>% 
  rename_at(vars(aa1, aa2), ~ if_else(. == "aa1", "aa2", "aa1"))



dat_formatted_none <- rbind(dat_formatted_fwd_none, dat_formatted_rev_none)

rbh_sims <- dat_formatted_none %>% select(protein, seq1, seq2, aa1, aa2, sim)

dat_formatted_ref_none <- dat_formatted_none %>%
  left_join(refs)


none_lab1 <- subset(dat_formatted_ref_none, label == 1)
none_lab1$set = "none_1"
mncm_lab1 <- subset(dat_formatted_ref, label==1) # %>% select(-sim) %>% left_join(rbh_sims)
mncm_lab1$set = "mncm_1"
none_labNA <- subset(dat_formatted_ref_none, is.na(label))
none_labNA$set = "none_NA#"
mncm_labNA <- subset(dat_formatted_ref, is.na(label)) # %>% select(-sim) %>% left_join(rbh_sims)
mncm_labNA$set = "mncm_NA"


none_labNA %>%
  full_join(mncm_labNA, by = c("protein", "fasta", "seq1", "seq2", "aa1", "aa2"))


mncm_1_median <-  median(mncm_lab1$sim, na.rm=TRUE)
mncm_NA_median <-  median(mncm_labNA$sim, na.rm=TRUE)
none_1_median <-  median(none_lab1$sim, na.rm = TRUE)
none_NA_median <-  median(none_labNA$sim, na.rm = TRUE)

alldat <- bind_rows(none_lab1, mncm_lab1) %>% 
  mutate(label = case_when(set == "mncm_1" ~ "MNCM",
                           set == "none_1" ~ "No MNCM"))

alldat_0 <- bind_rows(none_labNA, mncm_labNA) %>%
  mutate(label = case_when(set == "mncm_NA" ~ "MNCM",
                           set == "none_NA" ~ "No MNCM")) 

a <- alldat %>%
  ggplot(aes(sim, fill=label)) +
  geom_histogram(binwidth = 0.01) +
  facet_grid(~fct_rev(label)) + 
  theme(strip.background = element_rect(fill = "white"), legend.position="none") + 
  xlab("Cosine similarity") + 
  ylab("Count") +
  ggtitle("True positive RBHs") +  
  scale_fill_manual(values=c(palette_OkabeIto[1], palette_OkabeIto[2])) + 
  geom_vline(data = alldat[alldat$set == "none_1",], aes(xintercept=none_1_median), colour="black", linetype="dashed") + 
  geom_vline(data = alldat[alldat$set == "mncm_1",], aes(xintercept=mncm_1_median), colour="black", linetype="dashed") +
  xlim(0.5, 1)
a

b <- alldat_0 %>% 
  ggplot(aes(sim, fill=label)) +
  geom_histogram(binwidth = 0.01) +
  facet_grid(~fct_rev(label)) + 
  theme(strip.background = element_rect(fill = "white"),  legend.position="none") + 
  xlab("Cosine similarity") + 
  ylab("Count") +ggtitle("False positive RBHs") +  
  scale_fill_manual(values=c(palette_OkabeIto[1], palette_OkabeIto[2])) + 
  geom_vline(data = alldat_0[alldat_0$set == "none_NA",], aes(xintercept=none_NA_median), colour="black", linetype="dashed") + 
  geom_vline(data = alldat_0[alldat_0$set == "mncm_NA",], aes(xintercept=mncm_NA_median), colour="black", linetype="dashed") +
  xlim(0.5, 1)
b

alldat_0 %>%
  group_by(protein, label) %>%
   summarize(n = n()) %>%
  pivot_wider(names_from = "label", values_from = "n") %>%
  mutate(diff = `No MNCM` - MNCM) %>%
  arrange(desc(diff))

figure <- plot_grid(a, b, ncol = 1, labels = c("A", "B"), align = "v", label_size = 8)
figure

 figure %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures//hist.png", ., height=3, width = 3.5, units = "in", dpi = 600)
 figure %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures//hist.pdf", ., height=3, width = 3.5, units = "in")
```

Supplemental batch correction figure
```{r}
files_bc_combat <- fs::dir_ls("/scratch/gpfs/ia3026/vcmsa_tests/bc_combat_mnc/",glob = "*.csv")

dat_bc_combat <- read_csv(files_bc_combat, id="path",  show_col_types = FALSE, col_names = c("fasta", "seq1", "seq2", "aa1", "aa2", "sim"))

dat_formatted_fwd_bc <- dat_bc_combat %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% 
  filter(seq1 == 0 & seq2 == 1 )

dat_formatted_rev_bc <- dat_bc_combat %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]")  %>% filter(seq1 == 1 & seq2 == 0 ) %>% 
  rename_at(vars(seq1, seq2), ~ if_else(. == "seq1", "seq2", "seq1")) %>% 
  rename_at(vars(aa1, aa2), ~ if_else(. == "aa1", "aa2", "aa1"))



dat_formatted_bc_combat <- rbind(dat_formatted_fwd_bc, dat_formatted_rev_bc)
dat_formatted_ref_bc_combat <- dat_formatted_bc_combat %>%
  left_join(refs)


dat_evaluated_rbh_bc_combat <- dat_formatted_ref_bc_combat %>% filter(label == 1)
  #rowwise() %>%
  #mutate(true_rbh = case_when(label == 1 && rank_aa1 == 1 &&  rank_aa2 == 1 ~ "true_rbh"))


# Check for first position
dat_evaluated_rbh_bc_combat %>%
  group_by(label) %>%
  summarize(mean = mean(sim))

dat_evaluated_rbh_bc_combat %>%
  group_by(label) %>%
  summarize(n = n()) #%>%
  #filter(true_rbh == "true_rbh")


bc_eval_bc_combat <- dat_evaluated_rbh_bc_combat %>%
  filter(label ==1) %>%
  #filter(true_rbh == "true_rbh") %>%
  group_by(protein, label) %>%
  summarize(n = n())


bc_eval_bc_combat$set <- "BC"

counts <- data.frame(table(dat_formatted_ref_bc_combat$protein)) 
colnames(counts) <- c("protein", "freq")
total_df_bc <- left_join(bc_eval_bc_combat, counts)

bc_eval$set <- "MNCM"

counts <- data.frame(table(dat_formatted_ref$protein)) 
colnames(counts) <- c("protein", "freq")
total_df <- left_join(bc_eval, counts)

precision_df <- bind_rows(total_df_bc, total_df)
recall_df <- bind_rows(bc_eval_bc_combat %>%
  left_join(total_aligned_positions), bc_eval %>%
  left_join(total_aligned_positions)) %>% mutate(prop_correct=n/aligned_positions)

#adding in seq sim information before batch correction
files_seqsim <- fs::dir_ls("/scratch/gpfs/ia3026/vcmsa_tests/seqsim/",glob = "*.csv")

dat_seqsim <- read_csv(files_seqsim, id="path",  show_col_types = FALSE)

dat_formatted_seqsim <- dat_seqsim %>%
  mutate(path  = basename(path)) %>%
  separate(path, into = c("protein"), sep= "[.]") %>% filter(seq1 == 0 & seq2 == 1 )

prec_rec_df <- full_join(precision_df, recall_df, by = c('protein', 'set'), suffix = c("_precision", "_recall")) %>%
  merge(dat_formatted_seqsim) %>% mutate(f1 = 2*(prop_correct* prop_correct)/ (prop_correct + prop_correct))


prec_rec_sims <- left_join(prec_rec_df, homstrad_sims) 

f1_supp <- prec_rec_sims%>% filter(set %in% c("MNCM", "BC")) %>% select(protein, f1, set, same, seqsim, sim) %>%
  pivot_wider(names_from = "set", values_from = "f1") %>% mutate(diff = BC - MNCM) %>%
  mutate(above_zero = case_when(diff >= 0 ~ TRUE, diff <0 ~ FALSE)) %>%
  ggplot(aes(x=seqsim, y=diff, label=protein, color = above_zero)) + 
  geom_point(alpha=0.3) + xlab("Cosine similarity") + ylab("Batch corrected - uncorrected \n (F1 score)") +
  geom_vline(xintercept = 0.95, linetype="dashed") +
  geom_hline(yintercept = 0, linetype="dashed") +
  scale_color_manual(values = c(palette_OkabeIto[6], palette_OkabeIto[5])) + 
  theme(legend.position="none") 
f1_supp

f1_supp_df <-  prec_rec_sims %>% 
  filter(set %in% c("MNCM", "BC")) %>% 
  select(protein, f1, set, same, seqsim, sim) %>%
  pivot_wider(names_from = "set", values_from = "f1") %>% mutate(diff = BC - MNCM) %>%
  mutate(above_zero = case_when(diff >= 0 ~ TRUE, diff <0 ~ FALSE))

f1_supp %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/f1_supp.png", ., height=1.5, width = 2, units = "in", dpi = 600)
f1_supp %>% ggsave("/scratch/gpfs/cmcwhite/aln_comparison/figures/f1_supp.pdf", ., height=1.5, width = 2, units = "in")

```



```{r}




```