R/sm_diversity.R

In manuelsmendoza/smmcdr3: CDR3 sequences analysis

sm_diversity <- function(input_file, output_dir = pwd()){
  pkgs <- list("vegan", "iNEXT", "ggplot2", "ggpubr")
  pkg_attach <- function(x){
    library(package = x, character.only = T)
  }
  lapply(pkgs, pkg_attach)
  
  file_name <- unlist(strsplit(x = input_file, split = "_"))
  var_reg <- grep(pattern = "IGHV", x = file_name)
  isotype <- max(grep(pattern = "^C", x = file_name))
  title_text <- paste(file_name[1], file_name[var_reg], file_name[isotype])
  tissue <- unlist(strsplit(x = unlist(strsplit(x = input_file, split = "_"))[1], split = "/"))[length(unlist(strsplit(x = unlist(strsplit(x = input_file, split = "_"))[1], split = "/")))]
  
  data <- read.table(file = input_file, header = TRUE, comment.char = "", sep = "\t")
  data <- data[,c(6:ncol(data))]
  
  shannon_text <- diversity(x = t(data), index = "shannon")
  simpson_text <- diversity(x = t(data), index = "simpson")
  invsimpson_text <- diversity(x = t(data), index = "invsimpson")
  
  diversity_data <- cbind.data.frame(shannon_text, simpson_text, invsimpson_text)
  colnames(diversity_data) <- c("Shannon", "Simpson", "InvSimpson")
  
  shannon_plot <-
    ggplot(data = diversity_data, aes(x = rownames(diversity_data), y = Shannon, fill = rownames(diversity_data))) +
    geom_bar(stat = "identity") +
    ggtitle(label = paste(tissue, file_name[var_reg], file_name[isotype])) +
    theme_pubr() +
    theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90))
  simpson_plot <-
    ggplot(data = diversity_data, aes(x = rownames(diversity_data), y = Simpson, fill = rownames(diversity_data))) +
    geom_bar(stat = "identity") +
    ggtitle(label = paste(tissue, file_name[var_reg], file_name[isotype])) +
    theme_pubr() +
    theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90))
  invsimpson_plot <-
    ggplot(data = diversity_data, aes(x = rownames(diversity_data), y = InvSimpson, fill = rownames(diversity_data))) +
    geom_bar(stat = "identity") +
    ggtitle(label = paste(tissue, file_name[var_reg], file_name[isotype])) +
    theme_pubr() +
    theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90))
  
  diversity_plot <- ggarrange(plotlist = list(shannon_plot, simpson_plot, invsimpson_plot), ncol = 3, nrow = 1)
  
  write.table(x = diversity_data, file = paste0(output_dir, "/", paste0(tissue, file_name[var_reg], file_name[isotype]), "-diversity_indexes.tsv"), append = FALSE, sep = "\t", row.names = TRUE, col.names = TRUE)
  ggsave(plot = diversity_plot, path = output_dir, filename = paste0(paste0(tissue, file_name[var_reg], file_name[isotype]), "-diversity_indexes.png"), device = "png", width = 22, height = 17, units = "cm", dpi = 350)
  print("Diversity indexes exported OK!!!")
}