## Obtain  summary by otu of the frequencey

### read the table created with the barplot function where we kept all unassigned features of all samples
all_unassigned_table <- read.table("r_figures/quantitative_barplots/Absolute/Absolute/Table/Kingdom_Unclassified_Absolute>20_study_name_noyes_CURML_sample_source_abundancy_table.tsv", sep = "\t", header = T)

### Use summarise to have only once each OTU and the sum of the abundance for each
summarized_unassiged_table <- all_unassigned_table %>%
  group_by(OTU) %>%
  summarise(sum_abudance = sum(Abundance))

### Filter to keep only the ones with  more than 20 reads, to reduce the number of sequences to blast
filtered_unassigned_table <- summarized_unassiged_table %>%
  filter(sum_abudance >= 20)

### Read the fasta file of all sequences denoised by dada
all_fasta_seq <- read_fasta(file_path = "1_dada_denoised/dna-sequences.fasta")

### Transform the rownames as first column, needed to join with deplyr
all_fasta_seq <- as.data.frame(all_fasta_seq)
all_fasta_seq <- rownames_to_column(all_fasta_seq, var = "OTU")

### Keep only the sequences in the filtered table of all unassigned sequences
filtered_seq <- semi_join(all_fasta_seq, filtered_unassigned_table, by = c("OTU"))

### Write in into a fasta format table
df.fasta = dataframe2fas(filtered_seq, file="filtered_seq.fasta")

   
## Rename the sequence column in the all_fasta_seq table
names(all_fasta_seq)[2]<-paste("sequences")

## Calculate the sequences length
sequence_length <- all_fasta_seq %>%
  mutate(length = nchar(as.character(all_fasta_seq$sequences)))

## Get a table with the overall abundance of all features et their assignment
all_sequences_counts <- physeq_df %>%
  group_by(OTU, Kingdom) %>%
  summarise(sum_abundance = sum(Abundance))

## Join the two just created table with all sequences and abundance
joined_sequences_tables <- full_join(all_sequences_counts, sequence_length, by = "OTU")

## Join by length and by kingdom retaining the sum of the abundance for plotting
joined_sequences_tables_length_sum <- joined_sequences_tables %>%
  group_by(length, Kingdom) %>%
  summarise(length_sum_abundance = sum(sum_abundance))

## Create a version without the 0
joined_sequences_no_zero <- joined_sequences_tables_length_sum %>%
  ungroup() %>%
  filter(length_sum_abundance > 0)
  ### Write it as tsv
  write.table(x = joined_sequences_no_zero, file = "length_of_sequences.tsv", append = F, quote = F, sep = "\t", row.names = F, col.names = TRUE)

## Plot the size by Kingdom
read_length_plot <- joined_sequences_tables_length_sum %>%
  ggplot(aes(x = length, y = length_sum_abundance, fill = Kingdom)) +
  geom_col(position = "stack", width = 1) 

##Save this plots 
ggsave(read_length_plot, filename = "amplicons_length_distribution.png", width = 10, height = 5)

## Zoom to lower < values for better visualization

read_length_plot_zoomed <- read_length_plot + coord_cartesian(
xlim = c(250, 450), ylim = c(0, 50000))

##Save this plots 
ggsave(read_length_plot_zoomed, filename = "amplicons_length_distribution_zoomed.png", width = 10, height = 5)

## Function to load in R environnement the taxonomy table, the metadata table, the features counts table, the tree, the Shannon and the PCO objects


load_objects_fct <- function(taxonomy_class_table, replace_empty_tax = TRUE, features_counts_table , metadata_table, tree, shannon, pco) {

    
  ## Qiime2 .qza features counts table
    
    ### Load Qiime2 .qza counts table into R environment
    features_counts_table<-read_qza(features_counts_table)
    
    ### Shows the unique identifier of the file
    print("features counts table")
    print(features_counts_table$uuid)
    
    ### Print information about the object
    names(features_counts_table)
      
    ### Shows the unique identifier of the file
    print("features count table")
    print(features_counts_table$uuid)
    
    ### Show the first 5 samples and features
    print(features_counts_table$data[1:5,1:5]) 
    
    ### Save this table in global environnement
    features_counts_table <<- features_counts_table

    
}



## Run the function to load all needed object into R

load_objects_fct( features_counts_table =  "1b_vsearch_joining/dereplicate/dereplicated_table.qza")

## Create a phyloseq object from what has been loaded in the previous chunks

  ### Generate a function to create a phyloseq object
  create_phyloseq_fct <- function(physeq_name = "physeq", melted_df_name = "physeq_df", otu_table){

      ### Create the phyloseq object
      physeq_obj<-phyloseq(
      otu_table(otu_table$data, taxa_are_rows = T))
      
      ### Assign the object into the environnement 
      assign(value = physeq_obj, x = physeq_name, envir =  .GlobalEnv)
      
      ### Melt the physeq objet into a dataframe with one row per feature.id and per sample, needed later
      physeq_df <- psmelt(physeq_obj)
      
      ### Assign the object into the environnement 
      assign(value = physeq_df, x = melted_df_name, envir =  .GlobalEnv)
      
    
    }

  ### Use the function
  create_phyloseq_fct(physeq_name = "physeq_vsearch", melted_df_name = "physeq_df_vsearch", otu_table = features_counts_table)

### Read the fasta file of all sequences denoised by dada
all_fasta_seq_vsearch <- read_fasta(file_path = "1b_vsearch_joining/dereplicate/dna-sequences.fasta")

### Transform the rownames as first column, needed to join with deplyr
all_fasta_seq_vsearch <- as.data.frame(all_fasta_seq_vsearch)
all_fasta_seq_vsearch <- rownames_to_column(all_fasta_seq_vsearch, var = "OTU")

## Rename the sequence column in the all_fasta_seq table
names(all_fasta_seq_vsearch)[2]<-paste("sequences")


## Separate the first column where the feature ID and the sample name are merged
all_fasta_seq_vsearch <- separate(all_fasta_seq_vsearch, sep = " ", col = OTU, into = "OTU")

## Calculate the sequences length
sequence_length_vsearch <- all_fasta_seq_vsearch %>%
  mutate(length = nchar(as.character(all_fasta_seq_vsearch$sequences)))

## Get a table with the overall abundance of all features et their assignment
all_sequences_counts_vsearch <- physeq_df_vsearch %>%
  group_by(OTU) %>%
  summarise(sum_abundance = sum(Abundance))

## Join the two just created table with all sequences and abundance
joined_sequences_tables_vsearch <- full_join(all_sequences_counts_vsearch, sequence_length_vsearch, by = "OTU")

### Write it as tsv
  write.table(x = joined_sequences_tables_vsearch, file = "length_of_sequences_vsearch.tsv", append = F, quote = F, sep = "\t", row.names = F, col.names = TRUE)
  
  
## Join by length and by kingdom retaining the sum of the abundance for plotting
joined_sequences_tables_length_sum_vsearch <- joined_sequences_tables_vsearch %>%
  group_by(length) %>%
  summarise(length_sum_abundance = sum(sum_abundance))
  

## Plot the size by Kingdom
read_length_plot <- joined_sequences_tables_length_sum_vsearch %>%
  ggplot(aes(x = length, y = length_sum_abundance)) +
  geom_col(position = "stack", width = 1) 

##Save this plots 
ggsave(read_length_plot, filename = "amplicons_length_distribution_vsearch.png", width = 10, height = 5)

## Zoom to lower < values for better visualization

read_length_plot_zoomed <- read_length_plot + coord_cartesian(
xlim = c(0, 600), ylim = c(0, 50000))

##Save this plots 
ggsave(read_length_plot_zoomed, filename = "amplicons_length_distribution_zoomed_vsearch.png", width = 10, height = 5)