<jupyter_start><jupyter_text>GoalJacobo de la Cuesta-Zuluaga.While `PanX` provides a fancy visualization tool, it can be a little restrictive, therefore, I will perform some of the summary statistics and plotting in the present notebook.  Init<jupyter_code># Libraries
library(tidyverse)
library(data.table)
library(UpSetR)
library(stringr)
library(reshape2)
library(RColorBrewer)
library(phangorn)
library(cowplot)
library(forcats)
library(vegan)
# R misc functions
source("/ebio/abt3_projects/vadinCA11/data/V11/R_utils/utils.R")
# Load itol scripts
itol_dir = "/ebio/abt3_projects/vadinCA11/data/V11/itol"
itol_code = file.path(itol_dir, "itol.R") 
source(file = itol_code)<jupyter_output><empty_output><jupyter_text>Var<jupyter_code># Plots dir
plots_dir = "/ebio/abt3_projects/vadinCA11/data/V11/out_plots"
# Read gene cluster table
gene_clust_raw = read.delim("/ebio/abt3_projects/vadinCA11/data/V11/panx_output/tables/gene_clust_summary.txt", 
                               stringsAsFactors=FALSE)

metadata_table = read.table("/ebio/abt3_projects/vadinCA11/data/V11/assemblies_metadata/Assembly_metadata_curated.txt",
                              sep = "\t", dec = ".", header = T, row.names = 1,
                              as.is = T, comment.char = "") %>% 
                            rownames_to_column(var = "Taxonomy")

# InterPro Annotation
# Gene Clusters
Interpro_raw = "/ebio/abt3_projects/vadinCA11/data/V11/cluster_annot/consensus.fasta.tsv" %>%
    read.delim(sep = "\t", na.strings=c("","NA"),  header = F)

colnames(Interpro_raw) = c("Cluster", "MD5", "Len", "Database", 
                           "Accession", "Description", "Start", 
                           "End", "e_Val", "Status", "Date", 
                           "Interpro_accession",  "Interpro_annot")

# Core methanogenesis proteins
Methano_core_annot = "/ebio/abt3_projects/vadinCA11/data/V11/assemblies_metadata/Methanogenesis_arCOGs.txt" %>% 
    read_tsv() %>% 
    column_to_rownames("X1") %>% 
    t() %>% 
    as.data.frame %>%
    select(arCOGs, Distribution, Product)
rownames(Methano_core_annot) = NULL

Methano_core_count = "/ebio/abt3_projects/vadinCA11/data/V11/assemblies_metadata/Methanogenesis_counts.txt" %>%
    read_tsv() %>%
    select(-X1) %>%
    select(Assembly = Genomes, everything())

# Phylogenetic tree 
# Load phylogenomic tree
#phylo_tree_file = file.path("/ebio/abt3_projects/vadinCA11/data/V11/anvio_output/phylogenomic_tree.txt")
phylo_tree_file = file.path("/ebio/abt3_projects/vadinCA11/data/V11/phylophlan/trees/RAxML_bipartitions.methanomassilii.tre")
phylogenomic_tree_raw = read.tree(phylo_tree_file)
phylogenomic_tree_raw$node.label = c() # solves the "Labels are not unique" message


# Root tree and remove root for downstream analyses
outgroup_clades = c("GCA_900090055", "GCA_000195915", "GCA_900176435", "GCA_000152265", "GCA_001402945")
phylogenomic_tree_raw = phylogenomic_tree_raw %>% root(., outgroup = outgroup_clades)
phylogenomic_tree = phylogenomic_tree_raw %>% drop.tip(., outgroup_clades)
# Name and symbol of COG categories
# From http://clovr.org/docs/clusters-of-orthologous-groups-cogs/
COG_name = c("RNA processing and modification" , "Chromatin structure and dynamics" , 
             "Energy production and conversion" , "Cell cycle control, cell division, chromosome partitioning" , 
             "Amino acid transport and metabolism" , "Nucleotide transport and metabolism" , 
             "Carbohydrate transport and metabolism" , "Coenzyme transport and metabolism" , 
             "Lipid transport and metabolism" , "Translation, ribosomal structure and biogenesis" ,
             "Transcription" , "Replication, recombination and repair" ,
             "Cell wall/membrane/envelope biogenesis" , "Cell motility" , 
             "Post-translational modification, protein turnover, and chaperones" , 
             "Inorganic ion transport and metabolism" ,
             "Secondary metabolites biosynthesis, transport, and catabolism" , 
             "Signal transduction mechanisms" , 
             "Intracellular trafficking, secretion, and vesicular transport" , "Defense mechanisms" , 
             "Extracellular structures" , "Nuclear structure" , "Cytoskeleton" ,
             "General function prediction only", "Function unknow" , "No annotation")


COG_symbol = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", 
               "N", "O", "P", "Q", "T", "U", "V", "W", "Y", "Z", "R", "S", "X")

COG_definition = data.frame(Symbol = COG_symbol, Name = COG_name)

eggnog_header = c('query_name', 'seed_eggNOG_ortholog', 'seed_ortholog_evalue', 
                  'seed_ortholog_score', 'predicted_gene_name', 'GO_terms',
                  'KEGG_KOs', 'BiGG_reactions', 'Annotation_tax_scope', 
                  'OGs', 'bestOG|evalue|score', 'COG_cat', 'eggNOG_annot') 

#Load eggNOG mapper annotations
eggNOG_raw = "/ebio/abt3_projects/vadinCA11/data/V11/cluster_annot/eggNOG/Clusters.emapper.annotations" %>%
        read.delim(., header = F, sep = "\t", comment.char = "#", na.strings=c("","NA"))

colnames(eggNOG_raw) = eggnog_header

# Create a table retaining only one eggNOG annotation and COG category per cluster
eggNOG_clean = eggNOG_raw %>%
        mutate(NOG = factor(str_match(OGs, "[:alnum:]*@NOG"))) %>% 
        select(query_name, NOG, COG_cat, eggNOG_annot) %>%
        separate(COG_cat, into = c("COG_1", "COG_2", "COG_3"), sep = ",")
# Add Clade data to the metadata table and save
# Host Associated, Free Living and EXternal
HA_clade = c("GCA_002504405", "GCA_002506175", "GCA_002504495", "GCA_002506425", 
             "GCA_002503925", "GCA_002494805", "GCA_002498545", "GCA_002497475", 
             "GCA_002508595", "GCA_000800805", "GCA_002506905", "GCA_002496945", 
             "GCA_002495325", "GCA_002503545", "GCA_000300255", "V11_RL001", 
             "GCA_001481295", "GCA_900313075", "GCA_900314325", "GCA_002505245", 
             "GCA_002503785", "GCA_002506865", "GCA_002502545", "GCA_002509415", 
             "GCA_002498765", "GCA_002498785", "GCA_001421175", "GCA_002502965", 
             "GCA_002497155", "GCA_002509425", "GCA_002496785", "GCA_002505275", 
             "GCA_002506325", "GCA_002498605", "GCA_002509405", "GCA_002505345", 
             "GCA_002504645", "GCA_002502765", "GCA_002495665", "GCA_002506565", 
             "GCA_002508625", "GCA_002498425", "GCA_002498805", "GCA_002502925", 
             "GCA_002502465", "GCA_002498365", "GCA_002506995", "GCA_002494585", 
             "GCA_002509465", "GCA_002495495")

FL_clade = c("GCA_003153895", "GCA_003135935", "GCA_002506985", "GCA_002497075", 
             "GCA_000404225", "GCA_000308215", "GCA_002508545", "GCA_002494705", 
             "GCA_001421185", "GCA_002498285", "GCA_002495585", "GCA_002497995", 
             "GCA_002504525", "GCA_002496345", "GCA_002508585", "GCA_002508555", 
             "GCA_002502005", "GCA_002503495", "GCA_002505225", "GCA_002503645")

EX_clade = c("GCA_002496385", "GCA_002499085")

df_HA = data.frame(Assembly = HA_clade, 
                   clade = gsub(pattern = ".*", replacement = "HA", x = HA_clade, perl = T))
df_FL = data.frame(Assembly = FL_clade, 
                   clade = gsub(pattern = ".*", replacement = "FL", x = FL_clade, perl = T))
df_EX = data.frame(Assembly = EX_clade, 
                   clade = gsub(pattern = ".*", replacement = "EX", x = EX_clade, perl = T))

df_Clade = bind_rows(df_HA, df_FL, df_EX)

methanomassilii_data_raw = left_join(metadata_table, df_Clade, by = "Assembly")

methanomassilii_data_raw %>% dfhead
# Given that Interpro returns multiple hits from different DBs
# I will only use the most frequent Interpro ID for a given cluster
# I can always come back to the original Interpro table if needed

# Remove NAs from annotation and determine the most frequent ID
# If ties use the first
best_IPR_annot = Interpro_raw %>%
    filter(!is.na(Interpro_accession)) %>%
    group_by(Cluster, Interpro_accession) %>%
    tally() %>%
    mutate(rank = row_number(-n)) %>%
    filter(rank == 1) %>%
    ungroup %>%
    select(-c(n,rank))

# Create table with the annotation of a given ID
IPR_acc_annot = Interpro_raw %>% 
    select(Interpro_accession, Interpro_annot) %>%
    group_by(Interpro_accession, Interpro_annot) %>%
    filter(!is.na(Interpro_accession) & row_number() == 1) %>%
    ungroup()

# Annotation of each gene cluster
clust_IPR_annot = left_join(IPR_acc_annot, best_IPR_annot, by = "Interpro_accession") %>% 
    filter(!is.na(Cluster))

# Create a single table with all annotations
# Add eggNOG annotations
cluster_table = left_join(gene_clust_raw, eggNOG_clean, by = c("MSA" = "query_name"))
# Add InterPro annotations
cluster_table = left_join(cluster_table, clust_IPR_annot, by = c("MSA" = "Cluster"))


cluster_table = cluster_table %>% 
    select(Gene_ID, MSA, Gene_Name, Annotation, Count, NOG, COG_1, 
    eggNOG_annot, Interpro_accession, Interpro_annot, 
    starts_with("GCA"), starts_with("V11_RL001"))  %>%
    mutate(COG_1 = case_when(is.na(COG_1) ~ "X", TRUE ~ as.character(COG_1))) %>%
    mutate(COG_1 = factor(COG_1, levels = COG_symbol))

cluster_table %>% dfhead
# Export annotated cluster table
cluster_table_file = "/ebio/abt3_projects/vadinCA11/data/V11/cluster_annot/Gene_cluster_annots.tsv"
write_tsv(cluster_table, cluster_table_file)
# Number of annotated clusters
clust_annotations = cluster_table %>%
    select(MSA, matches("annot"), COG_1) %>%
    mutate(Annotation = na_if(Annotation, "hypothetical protein"))

All_clust = clust_annotations %>%
    summarise(n = n()) 

N_prokka = clust_annotations %>%
    filter(!is.na(Annotation)) %>%
    summarise(n = n())

N_eggNOG = clust_annotations %>%
    filter(!is.na(eggNOG_annot)) %>%
    summarise(n = n())

N_InterPro = clust_annotations %>%
    filter(!is.na(Interpro_annot)) %>%
    summarise(n = n())

N_COG = clust_annotations %>%
    filter((!COG_1 %in% c("S", "X"))) %>%
    summarise(n = n())

Annot_clust = clust_annotations %>%
    filter(!is.na(Interpro_annot) | !is.na(eggNOG_annot) | !is.na(Annotation) | (!COG_1 %in% c("S", "X"))) %>%
    summarise(n = n())

Annotator = c("Clusters", "Prokka", "COG", "eggNOG", "Interpro", "Combined") %>%
    factor(x = ., levels = .)

N_annot = bind_rows(All_clust, N_prokka,  N_COG, N_eggNOG, N_InterPro, Annot_clust) %>%
    mutate(Annotator)

N_annot %>%
    ggplot(aes(x = Annotator, y = n)) +
    geom_bar(stat = "identity") +
    labs(y = "Number of annotated clusters")


Annotator = c("Clusters", "Prokka", "eggNOG", "Interpro", "Combined") %>%
    factor(x = ., levels = .)

N_annot = bind_rows(All_clust, N_prokka, N_eggNOG, N_InterPro, Annot_clust) %>%
    mutate(Annotator)

N_annot %>%
    ggplot(aes(x = Annotator, y = n)) +
    geom_bar(stat = "identity")

clust_annotations %>%
    filter(!is.na(Interpro_annot) | !is.na(eggNOG_annot) | !is.na(Annotation)) %>%
    dfhead<jupyter_output>[1] 6745    5<jupyter_text>Colors<jupyter_code># Palette with colors of all qualitative palletes
# Extract all qual palettes
qual_col_pals = brewer.pal.info %>% 
    rownames_to_column() %>% 
    filter(rowname %in% c("Paired", "Set3"))
# Create single vector with all colors
col_vector = map2(qual_col_pals$maxcolors, qual_col_pals$rowname, brewer.pal) %>% 
    unlist %>%
    c(., "#696969", "#000000")
names(col_vector) = as.character(COG_symbol)<jupyter_output><empty_output><jupyter_text>Eukaryote-like proteins  Search clusters annotated as ELPs I will first perform a directed approach, in which I will look into the presence of eukaryote-like proteins and other symbiosis factors present in the different genomes. The factors I will investigate, and the search terms I will include are:* Adhesin-like proteins (ALPs): `Parallel beta-helix repeat`* Sel1 containing proteins (Sel1): `Sel1`* Listeria-Bacteroides repeat containing proteins (List-Bact): `Listeria-Bacteroides`* Tetratrico peptide repeats (TPRs): `Tetratrico` * Ankyrin repeats (ANKs): `Ankyrin`* Leucine-rich repeats (LRRs): `Leucine-rich`* Fibronectin type III (fn3) domains: `Fibronectin type III`* Laminin G domain: `Laminin G`* Bacterial Ig-like domains: `Ig-like`* YadA-like domain (Yersinia adhesin A): `YadA` * TadE-like domain: `TadE` * Invasion protein B (ialB): `Invasion protein B`<jupyter_code>ELPs = c("ALP" = "Parallel beta-helix repeat", "Sel1" = "Sel1", 
         "List_Bact" = "Listeria", "TPR" = "Tetratrico", 
         "ANK" = "Ankyrin", "LRR" = "Leucine-rich", "FN3" = "Fibronectin type III", 
         "Laminin_G" = "Laminin G", "Ig_like" = "Ig-like", 
         "YadA" = "YadA", "TadE" = "TadE", "ialB" = "Invasion protein B")
# Create a table with all clusters annotated as the above ELPs
filter_interpro = function(df, value, prot_name){
    
    clusters = df %>% 
        filter(str_detect(Interpro_annot, value)) %>%
        select(Cluster) %>%
        unique()
    
    clusters = clusters %>% 
        mutate(ELP = rep(prot_name, nrow(clusters)))
        
    clusters
}
# filter_interpro(Interpro_raw, "Sel1")
ELP_cluster_ID = map2_df(ELPs, names(ELPs), function(x, y) filter_interpro(Interpro_raw, x, y))
ELP_cluster_ID %>% dfhead<jupyter_output>[1] 519   2<jupyter_text>Filter gene cluster table<jupyter_code># Create a table with ELPs and cluster names with the counts on each genome
ELP_cluster_table = cluster_table %>% 
    filter(MSA %in% ELP_cluster_ID$Cluster) %>%
    left_join(ELP_cluster_ID, by = c("MSA" = "Cluster")) %>% 
    select(Gene_ID, MSA, ELP, starts_with("GCA"), starts_with("V11"))
ELP_cluster_table %>% dfhead
# Create table with the count of genes/clusters of each ELP in each genome 
ELP_Sum = ELP_cluster_table %>%
    select(-Gene_ID) %>%
    group_by(ELP) %>%
    summarise_if(is.numeric, .funs = sum) 

# Long form
ELP_Sum_long = ELP_Sum %>%
    gather(key = "Genome", value = "Count", 2:ncol(ELP_Sum))

# Add clade information
ELP_Sum_long = left_join(ELP_Sum_long, df_Clade, by = c("Genome" = "Assembly"))
ELP_Sum_long %>%
    group_by(ELP, clade) %>%
    summarise(mean = mean(Count, na.rm = T))

ELP_Sum_long %>%
    group_by(ELP, clade) %>%
    summarise(sd = sd(Count, na.rm = T))
# Save ELPs table
ELPs_File = "/ebio/abt3_projects/vadinCA11/data/V11/itol/ELPs.tsv"

# Transpose table
ELP_Sum_export = ELP_Sum %>% 
    select(-ELP) %>% 
    t %>% 
    as.data.frame %>% 
    rownames_to_column
    
colnames(ELP_Sum_export) =  c("Assembly", ELP_Sum$ELP)
ELP_Sum_export %>% dfhead

ELP_Sum_export %>% write_delim(., path = ELPs_File, delim = '\t', col_names = T)<jupyter_output>[1] 72 10<jupyter_text>Gene cluster size spectrum  Tables and stats<jupyter_code># Transform counts to presence/absence table
cluster_pres_abs = cluster_table %>%
    select(starts_with("GCA"), starts_with("V11")) %>%
    decostand(., method = "pa")

clust_tmp = cluster_table %>%
    select(-starts_with("GCA"), -starts_with("V11"))

cluster_pres_abs = bind_cols(clust_tmp, cluster_pres_abs)

# Add number of genomes of each clade in each cluster
cluster_pres_abs = cluster_pres_abs %>% 
    select(HA_clade) %>% 
    rowSums() %>%
    mutate(cluster_pres_abs, Count_HA = .)

cluster_pres_abs = cluster_pres_abs %>% 
    select(FL_clade) %>% 
    rowSums() %>%
    mutate(cluster_pres_abs, Count_FL = .)

cluster_pres_abs = cluster_pres_abs %>% 
    select(EX_clade) %>% 
    rowSums() %>%
    mutate(cluster_pres_abs, Count_EX = .)
# Cluster size
cluster_pres_abs %>% 
    count(Count) %>%
    mutate(prop = 100 * (n/sum(n)))<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code># Create a function to plot from the complete cluster table
cluster_spectrum = function(clust_tab, counts, categories){
    table_COGs = clust_tab %>% 
        group_by(eval(parse(text = counts)), eval(parse(text = categories))) %>%
        summarise(n_clusters = n())
    names(table_COGs) = c("Count", "Category", "n_clusters")
    
    table_COGs = table_COGs %>% 
        as.data.frame %>% 
        filter(Count != 0)
    
    max_break = table_COGs %>% pull(Count) %>% max
    
    table_COGs %>%
        ggplot(aes(x = Count, y = (n_clusters), color = Category)) +
            geom_line() +
                labs(x = "Genomes in cluster", y = "Number of clusters", color = "COG category") +
                scale_x_continuous(breaks=seq(0, max_break, 10)) + 
                scale_color_manual(values=col_vector) +
                theme_light() +
                theme(panel.grid.major = element_blank(), 
                      panel.grid.minor = element_blank(), 
                      legend.position = "none")
   
}
# Plots of gene cluster spectrum
# By clade and separating known and unknown clusters
# Mmassilii
Mmassilii_spectrum_u = cluster_pres_abs %>% 
    filter(COG_1 %in% c("X", "S")) %>% 
    cluster_spectrum(., "Count", "COG_1")

Mmassilii_spectrum_k = cluster_pres_abs %>% 
    filter(!(COG_1 %in% c("X", "S"))) %>% 
    cluster_spectrum(., "Count", "COG_1")

# HA
HA_spectrum_u = cluster_pres_abs %>% 
    filter(COG_1 %in% c("X", "S")) %>% 
    cluster_spectrum(., "Count_HA", "COG_1")

HA_spectrum_k = cluster_pres_abs %>% 
    filter(!(COG_1 %in% c("X", "S"))) %>% 
    cluster_spectrum(., "Count_HA", "COG_1")

# FL
FL_spectrum_u = cluster_pres_abs %>% 
    filter(COG_1 %in% c("X", "S")) %>% 
    cluster_spectrum(., "Count_FL", "COG_1")

FL_spectrum_k = cluster_pres_abs %>% 
    filter(!(COG_1 %in% c("X", "S"))) %>% 
    cluster_spectrum(., "Count_FL", "COG_1")

#EX
EX_spectrum_u = cluster_pres_abs %>% 
    filter(COG_1 %in% c("X", "S")) %>% 
    cluster_spectrum(., "Count_EX", "COG_1")

EX_spectrum_k = cluster_pres_abs %>% 
    filter(!(COG_1 %in% c("X", "S"))) %>% 
    cluster_spectrum(., "Count_EX", "COG_1")
# Print plot
options(repr.plot.width = 7, repr.plot.height = 3)
Mmassilii_spectrumplt = plot_grid(Mmassilii_spectrum_u, Mmassilii_spectrum_k,
          nrow = 1, labels = c("B", "C"), align = "v")
plot(Mmassilii_spectrumplt)
# Print plot
options(repr.plot.width = 10, repr.plot.height = 4)
plot_grid(HA_spectrum_u, HA_spectrum_k,
          nrow = 1, labels = "AUTO", align = "v")
# Print plot
options(repr.plot.width = 10, repr.plot.height = 4)
plot_grid(FL_spectrum_u, FL_spectrum_k,
          nrow = 1, labels = "AUTO", align = "v")
# Print plot
options(repr.plot.width = 10, repr.plot.height = 4)
plot_grid(EX_spectrum_u, EX_spectrum_k,
          nrow = 1, labels = "AUTO", align = "v")<jupyter_output><empty_output><jupyter_text>Clusters unique to each clade  Core genes of *Methanomassiliicoccales* and each cladeAccording to [Collins et al (2011)](https://doi.org/10.1186/1471-2105-12-S14-S9), the size distribution of gene cluster in bacterial and archaeal genomes follows the same pattern I observed for the complete order and each of the clades. I will focus on the gene clusters located in the right side of the spectrum, which correspond to genes present in most of the analyzed genomes. My reasoning is that these genes would provide better information regarding the differences in adaptation between clades. I will consider genes present in at least 80% of the genomes of a given cluster as the 'soft-core'.<jupyter_code># Determine the core of each clade and the complete order
core_cutoff = 0.8

core_Mmassilii = cluster_pres_abs %>% 
    filter(Count >= (72 * core_cutoff) & Count_FL > 0 & Count_HA > 0 & Count_EX > 0)
    
core_HA = cluster_pres_abs %>% 
    filter(Count_HA >= (length(HA_clade) * core_cutoff))

core_FL = cluster_pres_abs %>% 
    filter(Count_FL >= (length(FL_clade) * core_cutoff))

core_EX = cluster_pres_abs %>% 
    filter(Count_EX >= (length(EX_clade) * core_cutoff))
# Determine unique clusters of each clade
# i.e, clusters present in the core of a clade and absent in the others
unique_HA = core_HA %>%
    filter(Count_FL == 0 & Count_EX == 0)

unique_FL = core_FL %>%
    filter(Count_HA == 0 & Count_EX == 0)

unique_EX = core_EX %>%
    filter(Count_FL == 0 & Count_HA == 0)<jupyter_output><empty_output><jupyter_text>Plot proportion of COG categories by presence/absence<jupyter_code># Function to obtain the fraction labeled as each COG
COG_summary = function(clst, set_name) {
    clst %>% 
        group_by(COG_1) %>% 
        summarise(n = n()) %>%
        mutate(Proportion = n/sum(n), Clade = rep(set_name, nrow(.)))
}
COG_summaries = COG_summary(core_Mmassilii, "Methanomassilii-Shell")

COG_summaries = COG_summary(core_HA, "HA-Shell") %>% bind_rows(., COG_summaries)
COG_summaries = COG_summary(core_FL, "FL-Shell") %>% bind_rows(., COG_summaries)
COG_summaries = COG_summary(core_EX, "EX-Shell") %>% bind_rows(., COG_summaries)

COG_summaries = COG_summary(unique_HA, "HA-Unique") %>% bind_rows(., COG_summaries)
COG_summaries = COG_summary(unique_FL, "FL-Unique") %>% bind_rows(., COG_summaries)
COG_summaries = COG_summary(unique_EX, "EX-Unique") %>% bind_rows(., COG_summaries)
Shell_levels = c("Methanomassilii-Shell", "HA-Shell", "FL-Shell", "EX-Shell", "HA-Unique", "FL-Unique", "EX-Unique")
COG_summaries = COG_summaries %>% 
    mutate(Clade = factor(Clade, levels = Shell_levels))
COG_lvls = c("A", "B", "J", "K", "L", "D", "M", "N", "O", "T", "U", "V", "W", 
             "Y", "Z", "C", "E", "F", "G", "H", "I", "P", "Q", "X", "S") %>% rev
options(repr.plot.width = 10, repr.plot.height = 4)

COG_summaries_plt = COG_summaries %>%
    filter(str_detect(Clade, "Methano") | str_detect(Clade, "-Shell")) %>%
    mutate(Clade = str_remove_all(Clade, "-Shell"), 
           COG_1 = fct_relevel(COG_1, COG_lvls))

cog_barplt_pa = COG_summaries_plt %>%
    #filter(str_detect(Clade, "Methano") | str_detect(Clade, "-Shell")) %>%
    ggplot(., aes(x = Clade, y = Proportion, fill = COG_1)) + 
        geom_bar(stat="identity", position="fill") + 
        scale_fill_manual(values=col_vector) +
        theme_light() +
        labs(y = "Fraction of gene clusters", x = " ", fill = "COG\nCategoty") +
        coord_flip() +
        theme(legend.position="bottom", 
              panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank()) +
        guides(fill = guide_legend(title = "COG category", label.position = "left",
                                   nrow = 2, byrow = TRUE))
plot(cog_barplt_pa)<jupyter_output><empty_output><jupyter_text>Plot proportion of COG categories by gene cluster count<jupyter_code># Methanomassiliicoccales
Genomes_clade = methanomassilii_data_raw %>%
    select(Genome = Assembly, clade)

# Cluster count table in long format
core_Mmassilii_copies = cluster_table %>% 
    filter(MSA %in% core_Mmassilii$MSA) %>%
    select(MSA, COG_1, starts_with("GCA"), V11_RL001) %>% 
    gather(Genome, n, -c(MSA, COG_1))  %>%
    left_join(., Genomes_clade, by = "Genome")


# Total gene clusters in all genomes
total_clusters = core_Mmassilii_copies$n %>% sum

# Sum of gene clusters per category
Mmassilii_clusters_bycat = aggregate(n ~ COG_1, data = core_Mmassilii_copies, sum) %>%
    mutate(Proportion = 100*(n/total_clusters), Clade = "Methanomassilii")

# HA clade
# Cluster count table in long format
core_HA_copies = cluster_table %>% 
    filter(MSA %in% core_HA$MSA) %>%
    select(MSA, COG_1, one_of(HA_clade)) %>% 
    gather(Genome, n, -c(MSA, COG_1))  %>%
    left_join(., Genomes_clade, by = "Genome")

# Total gene clusters in all genomes
total_clusters = core_HA_copies$n %>% sum

# Sum of gene clusters per category
HA_clusters_bycat = aggregate(n ~ COG_1, data = core_HA_copies, sum) %>%
    mutate(Proportion = 100*(n/total_clusters), Clade = "HA")
    

# FL clade
# Cluster count table in long format
core_FL_copies = cluster_table %>% 
    filter(MSA %in% core_FL$MSA) %>%
    select(MSA, COG_1, one_of(FL_clade)) %>% 
    gather(Genome, n, -c(MSA, COG_1))  %>%
    left_join(., Genomes_clade, by = "Genome")

# Total gene clusters in all genomes
total_clusters = core_FL_copies$n %>% sum

# Sum of gene clusters per category
FL_clusters_bycat = aggregate(n ~ COG_1, data = core_FL_copies, sum) %>%
    mutate(Proportion = 100*(n/total_clusters), Clade = "FL")

# EX clade
# Cluster count table in long format
core_EX_copies = cluster_table %>% 
    filter(MSA %in% core_EX$MSA) %>%
    select(MSA, COG_1, one_of(EX_clade)) %>% 
    gather(Genome, n, -c(MSA, COG_1))  %>%
    left_join(., Genomes_clade, by = "Genome")

# Total gene clusters in all genomes
total_clusters = core_EX_copies$n %>% sum

# Sum of gene clusters per category
EX_clusters_bycat = aggregate(n ~ COG_1, data = core_EX_copies, sum) %>%
    mutate(Proportion = 100*(n/total_clusters), Clade = "EX")

# Combine tables
# Note change in clade name
# clade A -> "FL", clade B -> "HA" and clade C -> "EX"
core_proportion_copies = bind_rows(EX_clusters_bycat, 
                                   FL_clusters_bycat, 
                                   HA_clusters_bycat, 
                                   Mmassilii_clusters_bycat) 

core_proportion_copies %>% dfhead
options(repr.plot.width = 10, repr.plot.height = 4)

cog_barplt = core_proportion_copies %>%
    mutate(Clade = fct_relevel(Clade, rev(levels(Clade)))) %>%
    ggplot(., aes(x = Clade, y = Proportion, fill = COG_1)) + 
        geom_bar(stat="identity", position="fill") + 
        scale_fill_manual(values=col_vector) +
        theme_light() +
        labs(y = "Fraction of gene clusters", x = " ", fill = "COG\nCategoty") +
        coord_flip() +
        theme(legend.position="bottom", 
              panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank()) +
        guides(fill = guide_legend(title = "COG category", label.position = "right",
                                   nrow = 2, byrow = TRUE))

plot(cog_barplt)<jupyter_output><empty_output><jupyter_text>COG categories stats<jupyter_code># Mean of each COG category
COG_summaries_stats = core_proportion_copies %>% 
    filter(Clade != "Methanomassilii") %>%
    group_by(COG_1) %>% 
    summarise(mean_cog = mean(Proportion, na.rm = T), 
              sd_cog = sd(Proportion,  na.rm = T)) %>%
    arrange(-mean_cog)

COG_summaries_stats %>% filter(mean_cog >= 0.01)
# Unknowns per clade
core_proportion_copies %>%
    filter(COG_1 %in% c("X", "S")) %>%
    group_by(Clade) %>%
    summarise(Unkown = sum(Proportion))
# Mean COG on complete dataset
genome_cluster_prop %>%
    group_by(COG_1) %>%
    summarise(mean = mean(prop), sd = sd(prop))
# Mean unknown COGs on complete dataset
genome_cluster_prop %>%
    filter(COG_1 %in% c("S", "X")) %>%
    select(-prop) %>%
    spread(COG_1, count) %>%
    mutate(unknowns = 100*((S + X)/total)) %>%
    summarise(mean = mean(unknowns), sd = sd(unknowns))<jupyter_output><empty_output><jupyter_text>Barplot and Size distribution combined plot<jupyter_code>options(repr.plot.width = 10, repr.plot.height = 7)
# first align the top-row plot with the left-most plot of the bottom row
combined_left = align_plots(cog_barplt, Mmassilii_spectrum_u, align = 'v', axis = 'l')


# then build the top row
combined_top = plot_grid(combined_left[[2]], Mmassilii_spectrum_k, rel_widths = c(1.15, 1),
                        labels = c('A', 'B'), label_size = 12)

# then combine with the top row for final plot
combined_plot = plot_grid(combined_top, labels = c('', 'C'), combined_left[[1]],  label_size = 12, ncol = 1)
plot(combined_plot)
# Save plot
plot_file = file.path(plots_dir, "Combined_Cluster_COGs.png")
save_plot(filename = plot_file, plot = combined_plot,
          base_height = 7, base_width = 10, dpi = 300)<jupyter_output><empty_output><jupyter_text>Table with proportion of COG categories in the core of each clade<jupyter_code># Obtain the relative abundance of each COG category in the respective core of
# each genome

# Relabund of each COG category in the core of each genome of the HA clade
core_HA_relabund = cluster_table %>% 
    filter(MSA %in% core_HA$MSA) %>%
    select(MSA, COG_1, one_of(HA_clade)) %>%
    select(COG_1, starts_with("GCA"), V11_RL001) %>% 
    gather(Genome, Count, -COG_1) %>%
    aggregate(Count ~ COG_1 + Genome, ., sum) %>%
    spread(Genome, Count) %>%
    column_to_rownames("COG_1") %>%
    decostand(., method = "total", MARGIN = 2) %>%
    rownames_to_column("COG_1")

# Relabund of each COG category in the core of each genome of the FL clade
core_FL_relabund = cluster_table %>% 
    filter(MSA %in% core_FL$MSA) %>%
    select(MSA, COG_1, one_of(FL_clade)) %>%
    select(COG_1, starts_with("GCA")) %>% 
    gather(Genome, Count, -COG_1) %>%
    aggregate(Count ~ COG_1 + Genome, ., sum) %>%
    spread(Genome, Count) %>%
    column_to_rownames("COG_1") %>%
    decostand(., method = "total", MARGIN = 2) %>%
    rownames_to_column("COG_1")

# Relabund of each COG category in the core of each genome of the EX clade
core_EX_relabund = cluster_table %>% 
    filter(MSA %in% core_EX$MSA) %>%
    select(MSA, COG_1, one_of(EX_clade)) %>%
    select(COG_1, starts_with("GCA")) %>% 
    gather(Genome, Count, -COG_1) %>%
    aggregate(Count ~ COG_1 + Genome, ., sum) %>%
    spread(Genome, Count) %>%
    column_to_rownames("COG_1") %>%
    decostand(., method = "total", MARGIN = 2) %>%
    rownames_to_column("COG_1")

cores_relabund = full_join(core_HA_relabund, core_FL_relabund, by = "COG_1") %>%
    full_join(., core_EX_relabund, by = "COG_1") 

cores_relabund[is.na(cores_relabund)] = 0

# Transpose table
cores_COG_relabund = cores_relabund %>%
    gather(Genome, relabund, -COG_1) %>%
    spread(COG_1, relabund)

clades_Genome = methanomassilii_data_raw %>%
    select(Genome = Assembly, clade)

cores_COG_relabund = full_join(cores_COG_relabund, clades_Genome,  by = "Genome")
# Save tables
make_dir(dir = "/ebio/abt3_projects/vadinCA11/data/V11/COG_tables")
write_tsv(cores_COG_relabund, "/ebio/abt3_projects/vadinCA11/data/V11/COG_tables/COG_cores_relabund.tsv")
write_tsv(core_proportion_copies, "/ebio/abt3_projects/vadinCA11/data/V11/COG_tables/COG_cores_clade.tsv")<jupyter_output>Directory already exists: /ebio/abt3_projects/vadinCA11/data/V11/COG_tables<jupyter_text>Methanogenesis genes<jupyter_code>Methano_core_count %>% dfhead
Methano_core_annot %>% dfhead
options(repr.plot.width = 7, repr.plot.height = 10)
Methano_heatmap_df = Methano_core_count %>%
    gather(COG, Count, -Assembly) %>%
    mutate(COG = factor(COG, levels = Methano_core_annot$arCOGs), 
           Count = factor(Count), 
           Assembly = factor(Assembly, levels = phylogenomic_tree$tip.label))
    
Methano_core_heatmap = Methano_heatmap_df %>%
    ggplot(aes(x = COG, y = Assembly, fill = Count)) + 
    geom_tile() +
    theme_light() +
    #scale_fill_viridis_d() +
    #scale_fill_gradient(low = "white", high = "steelblue") +
    scale_fill_brewer(palette = "Greys") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot(Methano_core_heatmap)
# Itol heatmap
Methano_core_count %>%
    column_to_rownames("Assembly") %>% 
    itol_heatmap(df = ., dataset_label = "Core_methanogenesis", 
                    out_file = "Methanogenesis_heat.txt", out_dir =  itol_dir, 
                 dist_method='euclidean', color_scheme="bw")<jupyter_output>File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/Methanogenesis_heat.txt<jupyter_text>Create ITOL files<jupyter_code>ELP_Sum_export %>% dfhead
# Labels used in tree (Manually created)
itol_labels = file.path(itol_dir, "labels.txt") %>% 
    read.table(file = ., skip = 3, sep = ",", col.names = c("Assembly", "Label"))

methanomassilii_data = methanomassilii_data_raw %>%  
    left_join(., itol_labels, by = "Assembly") %>% 
    left_join(., ELP_Sum_export,  by = "Assembly") %>%
    mutate(GC = GC*100, Len = Len/1000000,
           source_1 = as.factor(source_2),
           source_2 = as.factor(source_2),
           source_3 = as.factor(source_3),
           source_manual = as.factor(source_manual),
           location = as.factor(location))

methanomassilii_data
# Create color strip of microbe source_1
methanomassilii_data %>% 
    column_to_rownames("Assembly") %>%
    select(source_1) %>% 
    itol_colorstrip(df = ., dataset_label = "Environment", 
                    out_file = "source1_colorstrip.txt", out_dir =  itol_dir)

# Create color strip of microbe source_2
methanomassilii_data %>% 
    column_to_rownames("Assembly") %>% 
    select(source_2) %>% 
    itol_colorstrip(df = ., dataset_label = "Source", 
                    out_file = "source2_colorstrip.txt", out_dir =  itol_dir)

# Create color strip of microbe source_3
methanomassilii_data %>% 
    column_to_rownames("Assembly") %>% 
    select(source_3) %>% 
    itol_colorstrip(df = ., dataset_label = "Sample", 
                    out_file = "source3_colorstrip.txt", out_dir =  itol_dir)

# Create color strip of microbe source manual
methanomassilii_data %>% 
    column_to_rownames("Assembly") %>% 
    select(source_manual) %>% 
    itol_colorstrip(df = ., dataset_label = "Isolation_source", 
                    out_file = "sourceM_colorstrip.txt", out_dir =  itol_dir)
# Change Colors in colorstrips

"sed -i -e 's/#FF0000/#d67622/' -e 's/#00FF00/#616163/' -e 's/#0000FF/#c2d8b9/' /ebio/abt3_projects/vadinCA11/data/V11/itol/source1_colorstrip.txt" %>% system

"sed -i -e 's/#FF0000/#A6CEE3/g' -e 's/#FF8B00/#1F78B4/g' -e 's/#E8FF00/#B2DF8A/g' -e 's/#5DFF00/#33A02C/g' -e 's/#00FF2E/#FB9A99/g' -e 's/#00FFB9/#E31A1C/g' -e 's/#00B9FF/#FDBF6F/g' -e 's/#002EFF/#FF7F00/g' -e 's/#5D00FF/#CAB2D6/g' -e 's/#E800FF/#6A3D9A/g' -e 's/#FF008/#FFFF99/g' /ebio/abt3_projects/vadinCA11/data/V11/itol/sourceM_colorstrip.txt" %>% system
methanomassilii_data %>% 
    select(Assembly, ALP, ANK, FN3, Ig_like, List_Bact, LRR, Sel1, TPR, YadA) %>%
    mutate_if(is.numeric, .f = function(x) x/max(x)) %>% 
    column_to_rownames("Assembly") %>% 
    itol_heatmap(df = ., dataset_label = "Adhesion_proteins", 
                 color_scheme = "bw",
                 out_file = "Prots_heatmap.txt", out_dir =  itol_dir)
              
methanomassilii_data %>% 
    select(Assembly, GC, Len, Genes) %>%
    mutate_if(is.numeric, .f = function(x) x/max(x)) %>% 
    column_to_rownames("Assembly") %>% 
    itol_heatmap(df = ., dataset_label = "Genome_Characteristics", 
                 color_scheme = "bw",
                 out_file = "Feats_heatmap.txt", out_dir =  itol_dir)
itol_file = function(Itol_df, to_rownames, to_itol, sufix = "shape", directory){
    var_name = deparse(substitute(to_itol))
    row_name = deparse(substitute(to_rownames))
    Itol_df %>%
    column_to_rownames(row_name) %>%
    select({{to_itol}}) %>%
    itol_externalshape(df = ., 
                       dataset_label = var_name, 
                       out_dir =  itol_dir,
                       out_file = glue::glue("{varname}_{sfx}.txt", 
                                       varname = var_name, sfx = sufix)) 
                       
}
itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = GC, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = Len, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = Genes, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = ALP, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = ANK, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = FN3, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = Ig_like, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = List_Bact, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = LRR, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = Sel1, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = TPR, sufix = "shape", directory = itol_dir)

itol_file(methanomassilii_data, to_rownames = Assembly, 
    to_itol = YadA, sufix = "shape", directory = itol_dir)<jupyter_output>File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/GC_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/Len_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/Genes_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/ALP_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/ANK_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/FN3_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/Ig_like_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/List_Bact_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/LRR_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/Sel1_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/TPR_shape.txt 
File written: /ebio/abt3_projects/vadinCA11/data/V11/itol/YadA_shape.txt<jupyter_text>Session info<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13927463934205203, 'fraction_numerical': 0.033086153431306456, 'mean_word_length': 3.9231381919553963, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 321, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '142874', 'license': 'MIT', 'n_tokens_mistral': 15069, 'n_tokens_neox': 13492, 'n_words': 2735, 'path': 'notebooks/.ipynb_checkpoints/12_Pangenome_statistics_new-checkpoint.ipynb', 'repo_name': 'leylabmpi/Methanomassilii', 'size': '37085'}