In Jorisvansteenbrugge/TcT: Trait-based Comparative Transcriptomics

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

#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')

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.

# 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

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

Plot_Background_Individual_Genes(bkgd.individual.Zscores)

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

Plot background of individual genes

Plot_Background_Modules(bkgd.traits)

Model_Module Figure

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

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

module_categories <- Get_module_categories(kegg_categories_script)
a <- getMetricDistModule(metric = 'PC', module_categories)

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

Plot_Categories()

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 ))
}

Jorisvansteenbrugge/TcT documentation built on Sept. 26, 2022, 6:50 a.m.