EBPR.Rmd

---
title: "Trait-based Comparative Transcriptomics Analysis"
author: "Joris van Steenbrugge & Ben Oyserman"
output: html_document
---
```{r,results='hide'}
library(TbasCO)
library(magrittr)

```
#Loading and Preprocessing Data


All of the preprocessing is handled by the `Pre_process_input` function.
This function consequtively performs the following tasks:

1. Data Loading
2. Normalization
3. Filtering


### Data loading
The data set should contain the following Elements:

* A column named **Gene**
* Multiple expression data columns that have at least the word **Sample** in it
* A column named **Annotation**
* A column named **Bin**

Additional columns may be present but are not used in the analysis.
Loading of data is done by assigning a file.path to the function's `filepath`
variable.


## Databases
```{r}
#load(file = 'data/kegg_module_2019_07_23.RData"')
database <- Combine_databases(kegg_brite_20191208, kegg_module_20190723)

```



### Normalization
All of the preprocessing is handled by the `Pre_process_input` function.
This function consecutively performs the following tasks:

1. Data Loading
2. Normalization -> Default normalization method between samples and by bin
3. Filtering -> the choice between Mean Absolute Deviation ('MAD') and based on Standard Deviation ('stdev')
```{r Preprocessing}

RNAseq.data <- Pre_process_input(file.path,
                                 database = annotation.db.path,
                                 normalize.method    = T,
                                 filter.method       ='MAD',
                                 filter.low.coverage = T)




genome.taxonomy.phylum <- list('8'  = 'p__Bacteroidetes',
                         '11' = 'p__Proteobacteria',
                         '16' = 'p__Proteobacteria',
                         '17' = '.',
                         '19' = 'p__Bacteroidetes',
                         '20' = '.',
                         '22' = 'p__Proteobacteria',
                         '25' = 'p__Bacteroidetes',
                         '26' = 'p__Actinobacteria',
                         '28' = 'p__Bacteroidetes',
                         '29' = 'p__Proteobacteria',
                         '31' = 'p__Proteobacteria',
                         '32' = 'p__Proteobacteria',
                         '36' = 'p__Actinobacteria',
                         '39' = 'p__Proteobacteria',
                         '42' = 'p__Bacteroidetes',
                         '45' = 'p__Bacteroidetes',
                         '48' = 'p__Proteobacteria',
                         '53' = 'p__Bacteroidetes'
                         )

```

### Define the distance metrics to be used
In this vignette, the Pearson Correlation and the Normalized Rank Euclidean Distance are used.
The user is free to implement different metrics, as long as the function parameters
are identical to the ones in the example PC and NRED functions.
In these examples, rowA and rowB are rows from RNAseq.data$table.
```{r Defining Distance Metrics}
# Calculates the Pearson Correlation
PC <- function(rowA, rowB, RNAseq.features){
  return(cor(as.numeric(rowA[RNAseq.features$sample.columns]),
             as.numeric(rowB[RNAseq.features$sample.columns])
             )
         )
}

# Calculates the Normalized Rank Euclidean Distance
NRED <- function(rowA, rowB, RNAseq.features) {
  r.A <- as.numeric(rowA[ RNAseq.features$rank.columns ])
  r.B <- as.numeric(rowB[ RNAseq.features$rank.columns ])
  return(
    sum((r.A - r.B) * (r.A - r.B))
  )
}

# Combine multiple distance metrics to complement each other.
distance.metrics <- list("NRED" = NRED,
                         "PC"   = PC)
```

### Main Pipeline
```{r Main Analysis}
bkgd.individual         <- Individual_Annotation_Background(RNAseq.data, 
                                                            N       = 5000, 
                                                            metrics = distance.metrics,
                                                            threads = 3)

bkgd.individual.Zscores <- Calc_Z_scores(bkgd.individual, distance.metrics)

bkgd.traits             <- Random_Trait_Background(RNAseq.data, 
                                                   bkgd.individual.Zscores,
                                                   N = 5000, 
                                                   metrics = distance.metrics,
                                                   threads = 4)

pairwise.distances      <- Calc_Pairwise_Annotation_Distance(RNAseq.data,
                                                             RNAseq.data$features$annotation.db,
                                                             distance.metrics,
                                                             bkgd.individual.Zscores,
                                                             show.progress = F,
                                                             threads = 4)

trait.attributes        <- Identify_Trait_Attributes(RNAseq.data = RNAseq.data, 
                                                     pairwise.distances = pairwise.distances,
                                                     threads = 4)

trait.attributes.pruned <- Prune_Trait_Attributes(trait.attributes, bkgd.traits, 
                                                  RNAseq.data,
                                                  p.threshold = 0.05,
                                                  pairwise.distances = pairwise.distances,
                                                  bkgd.individual.Zscores = bkgd.individual.Zscores)

sbs.trait.attributes    <- Traitattributes_To_Sbsmatrix(trait.attributes.pruned, RNAseq.data$features$bins)
```

### Plot background of individual genes
```{r bunch_o_figs, fig.height=8, fig.width=12}
Plot_Background_Individual_Genes(bkgd.individual.Zscores)

Plot_Metric_Comparison(bkgd.individual)
Plot_Redundancy_Traits(RNAseq.data)
```

### Plot background of individual genes
```{r bunch_o_figs, fig.height=8, fig.width=12}
Plot_Background_Modules(bkgd.traits)

```


### Model_Module Figure
```{r bunch_o_figs, fig.height=8, fig.width=12}
Module_Names <- RNAseq.data$features$trait_presence_absence[,'39'] %>% which(. == T) %>% names
margins =c(5,23)
Model_Module_List_All <- Model_Module(RNAseq.data, trait.attributes, '39', Module_Names, bkgd.traits)
Plot_Model_Module(Model_Module_List_All, '39', Module_Names, margins, sortbygenome="16")

```

### Figures of functional redundancy and linear regression
```{r bunch_o_figs, fig.height=8, fig.width=8}

TnA_redundancy <- Calc_TnA_redundancy(RNAseq.data)

# Sort based on attributes
Most_redundant_order <- TnA_redundancy[order(TnA_redundancy[,3]),1]



# Color by taxonomy
stc<- string.to.colors(genome.taxonomy.phylum)
# Color by taxonomy in the order based on redundancy
stc_highlights <- stc[match(rev(Most_redundant_order),names(stc))]

par(mfrow=c(4,5),mar=c(2,2,2,2))
#Plot_traits_vs_attributes()


#Figure 2-E
plot(as.numeric(TnA_redundancy[,3])~as.numeric(TnA_redundancy[,2]), col=stc, xlim=c(500,1100), ylim=c(100,400), pch = 19, xlab = "Attributes", ylab = "Traits")

# For some reason, legend is providing error
#legend(x=500, y=400, legend= c("Bacteroidetes","Proteobacteria", "Bacteria","Actinobacteria"), pch=19, col = unique(stc_highlights), cex=0.5)

as.numeric(Most_redundant_order)
j= 1
for (i in rev(Most_redundant_order)) {
  Plot_traits_vs_attributes_highlight(i,stc_highlights[j])
  j = j+1
}


```


Figure 2 C&D
```{r}
module_categories <- Get_module_categories(kegg_categories_script)
a <- getMetricDistModule(metric = 'PC', module_categories)
```

```{r}
attribute_list <- read.csv2("../data/attribute_list_combined.csv", stringsAsFactors = F) 
subset_attribute_list <- attribute_list[attribute_list$L2 == 'KEGG - Secretion Systems', ]
subset_names <- subset_attribute_list[,1]

sbs.trait.attributes.secretion <- Create_Filtered_SBS_Matrix(
                                        trait.attributes.pruned,
                                        RNAseq.data,
                                        subset_names)


Export_EdgeList(sbs.trait.attributes.secretion) %>% write.csv2(file = '/output/path.csv',
                                                               quote = F, row.names = F)

```


## Publication figure creating
```{r}
Plot_Categories()
```


##Benchmark number of genomes
```{r Benchmark number of genomes}
benchmark <- function(RNAseq.data, leave_n_out = 1){
  working_RNAseq.data <- RNAseq.data
  bins <- working_RNAseq.data$features$bins
  to_remove <- sample(bins, leave_n_out, replace = FALSE)
  working_RNAseq.data$features$bins %<>% .[-which(. %in% to_remove)]
  
  trait.attributes        <- Identify_Trait_Attributes(RNAseq.data = working_RNAseq.data, 
                                                     pairwise.distances = pairwise.distances,
                                                     threads = 6)

  trait.attributes.pruned <- Prune_Trait_Attributes(trait.attributes, bkgd.traits, 
                                                  working_RNAseq.data,
                                                  p.threshold = 0.05,
                                                  pairwise.distances = pairwise.distances,
                                                  bkgd.individual.Zscores = bkgd.individual.Zscores)
    
  count <-  trait.attributes.pruned %>% sapply(., length) %>% sum
  
  #average number of attributes per genome
  avg_ta_per_genome <- sapply(working_RNAseq.data$features$bins, function(bin){
    
    totals_current_genome <- sapply(trait.attributes.pruned, function(trait){
      
      has_ta <- sapply(trait, function(TA){
        
        return(bin %in% TA$genomes) 
      }) %>% as.numeric
      return(sum(has_ta))
    }) %>% as.numeric %>% sum
  }) %>% mean
  
  
  genome_lengths <- c()
  for(trait in trait.attributes.pruned){
    for (TA in trait){
      genomes <- TA$genomes
      genome_lengths <- c(genome_lengths, length(genomes))
    }
  }
  #Average number of genomes per attribute
  genomes_per_attributes.avg <- genome_lengths %>% mean
  return(list('count' = count,
              'avg_per_trait' = genomes_per_attributes.avg,
              'avg_per_genome' = avg_ta_per_genome) )
  
}


benchmarks <- c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10)

totals <- matrix(ncol = 4, nrow = 0)

for (leave_n_out in benchmarks){
  output <- benchmark(RNAseq.data, leave_n_out = 3)  
  
  totals <- rbind(totals, c(leave_n_out, output$count, output$avg_per_trait, output$avg_per_genome ))
}
```