marchantia_stress_updated.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"marchantia_stress_updated.ipynb","provenance":[{"file_id":"1CD90kG5BxodVmk3j0wGSx0Bb_g3acifr","timestamp":1656659597588}],"collapsed_sections":["bzX8KZDXIMXV","4qA5pIrH7VZW","exDaka827fIh","RdmYRmdc7oA4","cTKSAfIj7rcq","OZ8K1xcl7dfZ","kAj_ImMzXZwT","_ENdRwbI7iDw","DBM7cYbyxLT1"],"toc_visible":true},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","metadata":{"id":"f4f_0f6taBn9"},"source":["# Code to replicate analysis\n","\n","## Chapters\n","### 1. Miscellaneous preparation steps\n","1. Install non-default modules and upgrade modules\n","1. Mount Google Drive\n","1. Set paths (point to google drive folder to work in)\n","1. Download files (initial set up, skip if continuing)\n","1. Import modules, initialise paths\n","\n","### 2. Generating data for analysis\n","1. Differential Gene Expression (DESeq2)\n","1. Reconstructing the Gene Regulatory Network\n","\n","    1. Supp. Fig. 8A: Responsiveness of DEGs across experiments\n","    1. Reconstruction\n","    1. Supp. Fig. 9: Optimising parameters for filtering GRN\n","    1. Finalising the gene regulatory network; Fig 3A, visualised in cytoscape\n","    1. TF-only GRN: Fig 4A (visualised in Cytoscape) and Supp.Fig 11\n","1. Extract Arabidopsis GO and TFs\n","\n","### 3. Analysis and plotting\n","1. Figure 1 & Supp. Fig 1: Measurements and Student's t-test\n","1. Supp. Fig. 2: QC of RNA-seq data\n","1. Supp. Fig. 3: Volcano plots (DESeq2)\n","1. Supp. Fig. 4: Comparison of DEGs between two controls\n","1. Supp. Fig. 5: Upset plot (up- and down-regulated)\n","1. Supp. Fig. 6, 7: Intersection of DEGs across single and cross stresses (up- and down-regulated)\n","1. Figure 2: Summary of DEGs in Marchantia and inter stress comparisons\n","1. Supp. Fig. 10: Expression of GRN TFs across experiments (clustered)\n","1. Figure 3B: Expression of GRN TFs across experiments\n","1. Figure 3C: Specific expression in GRN TFs\n","1. Supp. Fig. 12: Influence of TFs\n","1. Supp. Fig. 13: Robustly responding second-level MapMan bins across the 7 abiotic stresses\n","1. Figure 4B, C, D; Supp. Fig. 14, 15: Bipartite networks for robustly expressed TFs and biological processes\n","1. Figure 5: Annotation of Ath orthologs with evidence from literature\n","1. Figure 6A, B: Effects of combined stress in terms of significant L2FC\n","1. Figure 6C: Classification of stress interactions\n","1. Figure 6D: Linear regression of all experiments\n","1. Figure 7A-H: Linear regression by stress\n"]},{"cell_type":"markdown","metadata":{"id":"63O3M3hPgwyv"},"source":["# 1. Miscellaneous preparation steps"]},{"cell_type":"markdown","metadata":{"id":"uGwkg7hkhMYY"},"source":["### 1.1 Install non-default modules and upgrade modules"]},{"cell_type":"code","metadata":{"id":"9THYinX_gpJ4"},"source":["# install non-default colab modules\n","# Restart runtime after installation and skip to next step\n","!pip install upsetplot\n","!pip install matplotlib --upgrade"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"DutthK8QhAkJ"},"source":["### 1.2 Mount Google Drive"]},{"cell_type":"code","metadata":{"id":"gLAcF40pZUM4"},"source":["from google.colab import drive\n","drive.mount('/content/gdrive')\n","!rm -rf /content/sample_data"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"50GdOAqr8D6D"},"source":["#@title 1.3 Set path {display-mode: \"form\"}\n","\n","#@markdown Enter the path of the directory you want to work in.\n","\n","drive_path = '/content/gdrive/My Drive/' #@param {type: 'string'}"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"5-ZWhJDzo7e6"},"source":["### 1.4 Download files (first time only)"]},{"cell_type":"code","metadata":{"id":"SFhEwOSCGeF_"},"source":["# Downloads necessary files to perform analyses [only need to be done once]\n","# https://gist.github.com/iamtekeste/3cdfd0366ebfd2c0d805 download raw files directtly from Google Drive\n","!wget --no-check-certificate -r \"https://drive.google.com/uc?id=1cbKgWbEWtstl_2_rb06tI_D-vseprPnT&export=download\" -O marchantia_stress.zip\n","\n","dir_path = drive_path + 'marchantia_stress/'\n","dir_path_safe = dir_path.replace(' ', '\\ ')\n","!mkdir $dir_path_safe\n","!unzip marchantia_stress.zip -d $dir_path_safe"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"u6GBaG6BQBR6"},"source":["###1.5 Import modules, set paths"]},{"cell_type":"code","metadata":{"id":"m3F4Mqz3P__H"},"source":["# import modules\n","import os\n","import string\n","%load_ext rpy2.ipython\n","import pandas as pd\n","import math\n","from matplotlib_venn import venn2\n","from matplotlib import pyplot as plt\n","from ast import literal_eval\n","import seaborn as sns\n","from collections import Counter\n","import random\n","from statsmodels.stats.multitest import multipletests\n","import numpy as np\n","from scipy.cluster.hierarchy import dendrogram, linkage\n","from scipy.spatial.distance import squareform\n","from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n","\n","dir_path = drive_path + 'marchantia_stress/'\n","dir_path_safe = dir_path.replace(' ', '\\ ')"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Dqb8m7ATg19d"},"source":["# 2. Generating data for analysis"]},{"cell_type":"markdown","metadata":{"id":"bzX8KZDXIMXV"},"source":["### 2.1 Differential Gene Expression"]},{"cell_type":"code","metadata":{"id":"DRBwbtjMVsSF"},"source":["# Making necessary directories for outputs\n","mpo_path = dir_path + \"prep_files/mpo/deseq/\"\n","mpo_path_safe = dir_path_safe + \"prep_files/mpo/deseq/\"\n","if not os.path.exists(mpo_path):\n","    !mkdir -p $mpo_path_safe\n","    print(\"Directories made: \" + mpo_path.replace('\\\\', ''))"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"a_uzqH_PUYQ-"},"source":["# Installing DESeq2\n","%%R\n","if (!requireNamespace(\"BiocManager\", quietly = TRUE))\n","    install.packages(\"BiocManager\")\n","\n","BiocManager::install(\"DESeq2\", ask = FALSE)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"rjetqXevenIe"},"source":["# To pull python variables\n","%R -i dir_path\n","%Rget dir_path\n","\n","%R -i dir_path_safe\n","%Rget dir_path_safe\n","\n","%R -i mpo_path\n","%Rget mpo_path"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"0HBke3tVIPk5"},"source":["# adapted from DESeq2_stressonly_phase1n2.R\n","%%R\n","# DESeq2 (Marchantia)\n","library('DESeq2')\n","library('RColorBrewer')\n","\n","sink(paste0(mpo_path, \"phase1n2_sum.txt\"), type=\"output\")\n","\n","raw_counts <- read.table(file = paste0(dir_path, 'prep_files/all_stress_raw.tsv'), sep = '\\t', header = TRUE)\n","raw_counts <- data.frame(raw_counts, row.names = 1)\n","\n","stresses <- c(\"controlH2\", \"controlD2\", \"H\", \"C\", \"M\", \"S\", \"L\", \"D\", \"N\",\n","              \"HS\", \"HM\", \"HN\", \"CS\", \"CM\", \"CN\", \"SM\", \"ML\", \"NL\", \"MN\",\n","              \"SD\", \"MD\", \"ND\", \"HD\", \"CD\", \"CL\", \"LS\", \"SN\")\n","colData = read.csv(paste0(dir_path, 'summary_files/all_stress.txt'), sep = '\\t', row.names=1, header = FALSE)\n","names(colData) <- c('condition')\n","\n","dds = DESeqDataSetFromMatrix(countData=raw_counts,\n","                             colData=colData,\n","                             design=~condition)\n","dds = DESeq(dds)\n","\n","y = 2\n","for (x in 1:2){\n","  for (i in y:length(stresses)){\n","    if (stresses[i] != stresses[x]){\n","      res = results(dds, contrast=c(\"condition\", stresses[i], stresses[x]))\n","      res = res[order(res$pvalue),]\n","      resSig = subset(res, res$padj < 0.05 & abs(res$log2FoldChange) > 1)\n","      resSig = resSig[ order(resSig$padj), ]\n","      print(paste(stresses[i], 'vs', stresses[x]))\n","      summary(res)\n","      summary(resSig)\n","      write.table(as.data.frame(res), file=paste(mpo_path, stresses[i], stresses[x], '_res.tsv', sep = ''),\n","                  quote=FALSE, sep='\\t', col.names = NA)\n","      write.table(as.data.frame(resSig), file=paste(mpo_path, stresses[i], stresses[x], '_resSig.tsv', sep = ''),\n","                  quote=FALSE, sep='\\t', col.names = NA)\n","    }\n","  }\n","  y = y + 1\n","}"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### 2.2 Reconstructing the Gene Regulatory Network"],"metadata":{"id":"K6m_o7t3qqud"}},{"cell_type":"markdown","source":["#### 2.2.1 Supp. Fig. 8A: Responsiveness of DEGs across experiments"],"metadata":{"id":"DBM7cYbyxLT1"}},{"cell_type":"code","source":["\"\"\"\n","DEG count across ALL stresses\n","Separated by TFs and non-TFs\n","\"\"\"\n","\n","import pandas as pd\n","from collections import Counter\n","import matplotlib.pyplot as plt\n","\n","genelist = \"G:/My Drive/Projects/Marchantia_2019/phase1n2/deseq/resSig_compiled.txt\"\n","tflist = \"G:/My Drive/Projects/Marchantia_2019/grn/TF_collate/PlantTFDB_Mpov5r1_prediction_plusTFDB.txt\"\n","\n","genedf = pd.read_csv(genelist, sep=\"\\t\") # 82982 DGEs\n","tflist = [x.split(\"\\t\")[0] for x in open(tflist, \"r\").readlines()] # 397 TFs\n","genedf['type'] = genedf[\"gene\"].apply(lambda x: \"TF\" if x in tflist else \"notTF\")\n","\n","sum(genedf.type == \"TF\") # 1341 occurrences\n","len(set(genedf.gene)) # 12442 genes /19421 genes\n","len(set(genedf.gene[genedf.type == \"TF\"])) # 252 unique TFs / 397 TFs\n","\n","genedf.to_csv(\"G:/My Drive/Projects/Marchantia_2019/grn/resSig_compiled_wType.txt\", sep = \"\\t\", index = False)\n","\n","condcount = Counter(genedf.gene)\n","cdict = {\"genecount\" : [], \"tfcount\": [], \"genes\" : [], \"TFs\" : []}\n","\n","for i in range(max(condcount.values())):\n","\tgenecount = [x for x in condcount if condcount[x] > i]\n","\tTFcount = [x for x in genecount if x in tflist]\n","\tcdict[\"genes\"].append(genecount)\n","\tcdict[\"TFs\"].append(TFcount)\n","\tcdict[\"genecount\"].append(len(genecount))\n","\tcdict[\"tfcount\"].append(len(TFcount))\n","\n","dfcount = pd.DataFrame(cdict)\n","dfcount.index.names = [\"InMoreThan\"]\n","dfcount.to_csv(\"G:/My Drive/Projects/Marchantia_2019/grn/count.txt\", sep = \"\\t\")\n","dfcount.plot.bar(y=[\"genecount\", \"tfcount\"], use_index = True, logy = True)\n","plt.savefig(dir_path + 'figures/' + 'FigS8A.png', dpi=600)"],"metadata":{"id":"v9DfNo0HxP9A"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["#### 2.2.2 Reconstruction"],"metadata":{"id":"cywXxXxY9l7F"}},{"cell_type":"code","source":["%%R\n","# Matrix trimming\n","# Aims:\n","#       1) Trim matrix to only contain genes that are differentially\n","#          expressed in > 5 conditions\n","# Outputs:\n","#       1) Trimmed TPM matrix\n","\n","library(readr)\n","\n","ofile <- file.path(dir_path, \"prep_files/all_stress_mt5.tsv\")\n","ifile <- file.path(dir_path, \"prep_files/all_stress.tsv\")\n","cfile <- file.path(dir_path, \"prep_files/count.txt\")\n","\n","tpm_mat <- as.matrix(read.table(ifile, header=TRUE, sep = \"\\t\",\n","                                row.names = 1, as.is=TRUE))\n","mat_col_ref <- read_tsv(file.path(paste0(dir_path, 'summary_files/all_stress.txt')), col_names = FALSE)\n","# https://cran.r-project.org/web/packages/comprehenr/vignettes/Introduction.html\n","nname <- paste(mat_col_ref$X2, \"_\", to_vec(for(x in str_split(mat_col_ref$X1, \"_\")) x[[2]]), sep = \"\")\n","names(nname) <- mat_col_ref$X1\n","\n","for (i in 1 : length(nname)){\n","  new <- nname[i]\n","  old <- names(nname[i])\n","  colnames(tpm_mat)[colnames(tpm_mat) == old] <- new\n","}\n","\n","count_mat <- as.matrix(read.table(cfile, header=TRUE, sep = \"\\t\",\n","                                row.names = 1, as.is=TRUE))\n","\n","valid_str <- gsub(\"\\\\[|\\\\]\",'', count_mat[6, 3]) #remove '[ and ]'\n","valid_genes <- unlist(strsplit(valid_str, \", \")) #split on ', '\n","pruned_mat <- tpm_mat[valid_genes,]\n","pruned_mat_dum <- rbind(pruned_mat, replicate(length(colnames(pruned_mat)), 0))\n","\n","write.table(pruned_mat_dum, ofile, sep = \"\\t\", col.names=NA)\n","\n","## Note! Remove dummy lines in matrix manually. Updated matrix provided for download."],"metadata":{"id":"obJs3ut7uXyc"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["! cd dir_path + \"/GRN_code\""],"metadata":{"id":"3zWss9LtmynW"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["%%R\n","# GRN Part 3 \n","# Started: 17 February 2022\n","# Aims:\n","#       1) Generate the linear regression models for gene-TF and TF-TF\n","#       2) Regularisations: LASSO/ elastic net\n","#       3) Intermediate outputs per response variable\n","#       4) Final output as a GRN with [goi_name, TF_name, coefficients(beta), gene_name (of TFs), relvar, s, lambda, cv]\n","\n","# Dependencies:\n","# 1) dircreater.r\n","# 2) wrap_elnetv4_Mpo.r\n","# 3) TFelnetv6_Mpo.r\n","# note: cd to /home/qiaowen/Marchantia_2019/scripts/GRN_code/ to access before running code to access the dependencies\n","\n","# Working paths (ChuckNorris)\n","ddir <- file.path(dir_path, \"GRN_code/deps\")\n","wdir <- file.path(dir_path, \"prep_files/Mpo_GRN_models\")\n","mat_path <- file.path(dir_path, \"prep_files/all_stress_mt5_nodum.tsv\")\n","tf_path <- file.path(dir_path, \"prep_files/PlantTFDB_Mpov5r1_prediction_plusTFDB.txt\")\n","\n","## Script to create the elastic net derived grn\n","source('deps/wrap_elnetv3_Mpo.r')\n","#library(edgeR)\n","library(comprehenr)\n","set.seed(2019)\n","options(stringsAsFactors = FALSE)\n","\n","cmode_name <- \"elnet\"\n","\n","extract_matrix <- function(mat, x) {\n","  x_names <- to_vec(for(z in colnames(mat)) if(grepl(x, z, fixed = TRUE) & !grepl('control', z, fixed = TRUE)) z)\n","  new_matrix <- subset(mat, select = x_names)\n","  return(new_matrix)\n","}\n","\n","#stresses <- c(\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all\") # uncomment and change K=3 in TFelnetv6_Mpo.r for 3 fold cross-validation. \n","stresses <- c(\"all\") # K=5 in TFelnetv6_Mpo.r for 5 fold cross-validation\n","\n","for (mat_type in stresses){\n","  print(paste(\"Building models for subset: \", mat_type))\n","  # load gene expression matrix and normalise (log transform followed by z transform)\n","  mat <- as.matrix(read.table(mat_path, header=TRUE, sep = \"\\t\", row.names = 1, as.is=TRUE))\n","  resdir <- file.path(wdir, cmode_name, mat_type)\n","  \n","  if (mat_type != \"all\"){\n","    mat <- extract_matrix(mat, mat_type)\n","  }\n","  \n","  # kick out genes that are completely '0' across all conditions (aftifact of subsetting)\n","  mat <- mat[rowSums(mat) != 0,]\n","  \n","  ###\n","  # Checks if matrix contains zeros\n","  if (sum(mat>0) > 0){\n","    print(\"Warning: '0' present in gene expression matrix.\")\n","    minval <- min(mat[mat > 0])\n","    print(paste(\"Minimum expression value:\", minval))\n","    print(\"Replacing zeros with 1e-12\")\n","    mat[mat == 0] <- 0.000000000001\n","  }\n","  \n","  # Log transforms matrix\n","  log_mat <- log(mat)\n","  #sum(colSums(log_mat == -Inf))\n","  gene_dat <- t(scale(t(log_mat)))\n","  \n","  #import TF annotation - names should be the same as rownames of the read data\n","  TF <- read.delim(tf_path, header = FALSE)\n","  tfs <- TF[, 1]\n","  tfs <- tfs[tfs %in% rownames(gene_dat)] # grab only TFs that are in current dataset\n","  #print(length(tfs))\n","  \n","  #Do elastic net regression analysis for each TF using all other TFs as predictors \n","  #THIS STEP TAKES VERY LONG - RUN ON SERVER OR wITH MORE THEN 8 THREADS\n","  elnet_res <- wrap_elnet(gene_dat, resdir=resdir, thrsh=0, tfs=tfs, parallel=64, cmode=cmode_name)\n","  elnet_all <- elnet_res[,c('Gene.ID', 'predicted', 'rel.coeff')]\n","  colnames(elnet_all) <- c('from', 'to', 'weight' )\n","  elnet_all <- elnet_all[order(abs(elnet_all$weight), decreasing = TRUE),]\n","  save(elnet_all, file=file.path(resdir,'elnet_all.obj'))\n","  write.table(elnet_all, file.path(resdir, 'elnet_all.txt'), sep = \"\\t\", col.names = NA)\n","}"],"metadata":{"id":"2tHtpakrrJYe"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["!cd ~"],"metadata":{"id":"rpolVmqYm7zC"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["#### 2.2.3 Supp. Fig 9: Optimising parameters for filtering GRN"],"metadata":{"id":"EBjT4hOexuUe"}},{"cell_type":"code","source":["%%R\n","# Finding the ideal cutoff to apply to the networks\n","# Aims:\n","#       1) General distribution of relative coefficient\n","#       2) Number of models per R^2 (Fig S9A)\n","#       3) Max coeff. per R^2\n","#       4) Number of nodes and edges per coeff cutoff\n","\n","library(ggplot2)\n","library(RColorBrewer)\n","\n","wdir <- file.path(dir_path, \"prep_files/Mpo_GRN_models\")\n","resname <- \"elnet.txt\"\n","stresses <- c(\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\")\n","cutoffs <- c(0.5, 0.625, 0.75)\n","\n","# load networks\n","read_mat <- function(mat_path){\n","  as.data.frame(read.table(mat_path, header=TRUE, sep = \"\\t\", row.names = 1, as.is=TRUE))\n","}\n","\n","elnet_res <- lapply(stresses, function(x){read_mat(paste(wdir, \"elnet\", x, resname, sep=\"/\"))})\n","names(elnet_res) <- stresses\n","\n","# 1) General relative coefficient distribution\n","for (cf in cutoffs){\n","  # apply relvar cutoff\n","  elnet_cut <- lapply(elnet_res, function(x){subset(x, relvar > cf)})\n","  names(elnet_cut) <- stresses\n","  #multiplot\n","  pdf(file.path(dir_path, \"prep_files\", paste(cf, \"_distribution_hist\", \".pdf\", sep=\"\")), width = 8.5)\n","  par(mfrow=c(3,3))\n","  for (s in stresses){\n","    hist(elnet_cut[[s]][,\"rel.coeff\"], xlab = \"Relative coefficient\", main = paste(cf, s, \"relative coefficient distribution\"))\n","  }\n","  dev.off()\n","}\n","\n","# 2) Number of models per R^2\n","cutoffs <- c(5:10)/10\n","cf_count <- matrix(0, nrow = length(cutoffs), ncol = length(stresses))\n","rownames(cf_count) <- cutoffs\n","colnames(cf_count) <- stresses\n","for (cf in cutoffs){\n","  # apply relvar cutoff\n","  elnet_cut <- lapply(elnet_res, function(x){subset(x, relvar > cf)})\n","  names(elnet_cut) <- stresses\n","  for (s in stresses){\n","    cf_count[as.character(cf), s] = length(unlist((unique(elnet_cut[[s]][\"predicted\"]))))\n","  }\n","}\n","pdf(file.path(dir_path, \"figures\", \"FigS9A.pdf\"), width = 7.5, height = 6)\n","matplot(cf_count, type = \"o\", lwd = 2, col = brewer.pal(length(stresses), \"Dark2\"),\n","        main = \"Number of models at various R^2 cutoffs\", xaxt = \"n\",\n","        lty = 1, xlab=\"R^2 cutoff\", pch=16) # \n","# Add X-axis\n","axis(side=1,at=1:nrow(cf_count),labels=cutoffs)\n","legend(\"topright\", legend = stresses, col=brewer.pal(length(stresses), \"Dark2\"), lty = 1, pch=16, lwd=2)\n","#xlab(\"R^2 cutoff\")\n","dev.off()\n","\n","# 3) Max coeff. per R^2 (absolute)\n","\n","pdf(file.path(dir_path, \"prep_files\", \"coeff_R2_lm.pdf\"))\n","par(mfrow=c(3,3))\n","for (s in stresses){\n","  unique_genes <- unlist(unique(elnet_res[[s]][\"predicted\"]))\n","  R2_val <- sapply(unique_genes, function(x){unlist(unique(subset(elnet_res[[s]], predicted == x, relvar)))})\n","  max_relcoeff <- sapply(unique_genes, function(x){max(abs(subset(elnet_res[[s]], predicted == x, rel.coeff)))})\n","  plot(R2_val, max_relcoeff, main = paste(\"Network\", s), xlab = \"Maximum relative coefficient\", ylab = \"R^2 value\", cex = 0.5)\n","  lm_obj <- lm(max_relcoeff ~ R2_val)\n","  abline(lm_obj, col = \"red\")\n","  mtext(paste(\"R^2:\", round(summary(lm_obj)$r.squared, 2)), 3, adj = 0.02, line = -1, cex = 0.5)\n","}\n","dev.off()\n","\n","# 4) Number of nodes and edges per coeff cutoff (global)\n","get_cutoffs <- function(cf){\n","  quantiles <- seq(0.1, 0.9, 0.1)\n","  node_count <- matrix(0, nrow = length(quantiles), ncol = length(stresses))\n","  edge_count <- matrix(0, nrow = length(quantiles), ncol = length(stresses))\n","  colnames(node_count) <- stresses\n","  rownames(node_count) <- quantiles\n","  colnames(edge_count) <- stresses\n","  rownames(edge_count) <- quantiles\n","  \n","  # apply relvar cutoff\n","  elnet_cut <- lapply(elnet_res, function(x){subset(x, relvar > cf)})\n","  names(elnet_cut) <- stresses\n","  # get quantiles\n","  for (s in stresses){\n","    nw_quantile <- quantile(unlist(abs(elnet_cut[[s]][, \"rel.coeff\"])), quantiles)\n","    for (q in names(nw_quantile)){\n","      quantile_cut <- subset(elnet_cut[[s]], abs(rel.coeff) >= nw_quantile[q])\n","      nodes <- length(unique(append(unique(quantile_cut[, \"predicted\"]), unique(quantile_cut[, \"Gene.ID\"]))))\n","      edges <- nrow(quantile_cut)\n","      rname <- as.character(as.numeric(substr(q, 1, nchar(q)-1))/100)\n","      node_count[rname, s] <- nodes\n","      edge_count[rname, s] <- edges\n","    }\n","  }\n","  # plotting (Node)\n","  matplot(node_count, type = \"l\", lwd = 2, col = brewer.pal(length(stresses), \"Dark2\"),\n","          main = paste(\"Number of nodes, R^2:\", cf), xaxt = \"n\",\n","          lty = 1:length(stresses), xlab=\"coefficient quantiles\", ylab = \"Number of nodes\")\n","  # Add X-axis\n","  axis(side=1,at=1:nrow(node_count),labels=names(nw_quantile))\n","  legend(\"bottomright\", inset=c(0.01 ,0.1), legend = stresses, col=brewer.pal(length(stresses), \"Dark2\"),\n","         lty = 1:length(stresses), cex = 0.8, pt.lwd = 2)\n","\n","  # plotting (Edges)\n","  matplot(edge_count, type = \"l\", lwd = 2, col = brewer.pal(length(stresses), \"Dark2\"),\n","          main = paste(\"Number of edges, R^2:\", cf), xaxt = \"n\",\n","          lty = 1:length(stresses), xlab=\"coefficient quantiles\", ylab = \"Number of edges\")\n","  # Add X-axis\n","  axis(side=1,at=1:nrow(edge_count),labels=names(nw_quantile))\n","  legend(\"topright\", legend = stresses, col=brewer.pal(length(stresses), \"Dark2\"),\n","         lty = 1:length(stresses), cex = 0.8, pt.lwd = 2)\n","}\n","\n","pdf(file.path(dir_path, \"prep_files\", \"network_count.pdf\"), width = 8, height = 11)\n","for (cf in cutoffs){\n","  par(mfrow=c(2,1))\n","  get_cutoffs(cf)\n","}\n","dev.off()"],"metadata":{"id":"sPjZc0ZZx23f"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["\"\"\"\n","Fig S9B -- Part 1\n","Check overlap of AGRIS regulatory network and elastic net\n","\"\"\"\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","from collections import Counter\n","\n","# Initialise Orthofinder output genes: orthogroup\n","OF_path = dir_path + \"/prep_files\" + '/Orthogroups.txt'\n","\n","gene_OF = {}\n","\n","with open(OF_path, 'r') as OF_file:\n","\tfor line in OF_file:\n","\t\tcontent = line.strip(\"\\n\").split(\": \")\n","\t\tog = content[0]\n","\t\tgene_list = content[1].split(\" \")\n","\t\tfor gene in gene_list:\n","\t\t\tgene_OF[gene] = og\n","\n","# initialise agris\n","agris_path = dir_path + 'prep_files/AtRegNet_confirmed.txt'\n","agris_out_path = dir_path + 'prep_files/AtRegNet_confirmed_OG.txt'\n","agris_og_nodup = dir_path + 'prep_files/AtRegNet_OG_nw.txt'\n","# TFLocus, TargetLocus: 1, 4\n","agris_df = pd.read_csv(agris_path, sep=\"\\t\", header=0, index_col=None)\n","agris_df[\"TFOG\"] = [gene_OF[x] if x in gene_OF else pd.NA for x in agris_df.TFLocus]\n","agris_df[\"TargetOG\"] = [gene_OF[x] if x in gene_OF else pd.NA for x in agris_df.TargetLocus]\n","agris_df.dropna(inplace=True, subset=[\"TargetOG\", \"TFOG\"])\n","agris_df.to_csv(agris_out_path, index=False, sep=\"\\t\")\n","\n","no_dup = agris_df.drop_duplicates(subset=[\"TargetOG\", \"TFOG\"])[[\"TargetOG\", \"TFOG\"]]\n","no_dup.to_csv(agris_og_nodup, index=False, sep=\"\\t\")\n","\n","no_dup_og = list(set(no_dup.TargetOG.to_list() + no_dup.TFOG.to_list()))\n","\n","# initialise elnet\n","elnet_dir = dir_path + \"/prep_files/Mpo_GRN_models/elnet/\"\n","stresses = [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\"]\n","elnet = \"elnet.txt\"\n","quantile = [i/10 for i in range(0,10)] # change\n","\n","\n","for i in stresses:\n","\t# add OG to main network file\n","\telnet_file = pd.read_csv(elnet_dir + i + \"/\" + elnet, header=0, sep=\"\\t\", index_col=0)\n","\telnet_file[\"predicted_GO\"] = [gene_OF[x] for x in elnet_file.predicted]\n","\telnet_file[\"Gene_GO\"] = [gene_OF[x] for x in elnet_file[\"Gene.ID\"]]\n","\telnet_file.to_csv(elnet_dir + i + \"/elnet_GO.txt\", index=False, sep=\"\\t\")\n","\t\n","\t# GO terms are found in the agris netreg network\n","\telnet_in_agris = elnet_file[(elnet_file.predicted_GO.isin(no_dup_og)) & (elnet_file.Gene_GO.isin(no_dup_og))]\n","\telnet_in_agris.to_csv(elnet_dir + i + \"/elnet_in_agris.txt\", index=False, sep=\"\\t\")\n","\t\n","\t# remove duplicated interactions\n","\telnet_nodup = elnet_in_agris.drop_duplicates(subset=[\"predicted_GO\", \"Gene_GO\"])[[\"predicted_GO\", \"Gene_GO\"]]\n","\t#elnet_nodup.to_csv(elnet_dir + i + \"/elnet_in_agris_GOnw.txt\", index=False, sep=\"\\t\")\n","\t\t\n","# The actual comparison\n","\n","no_dup[\"target_TF\"] = no_dup.TargetOG + \"_\" + no_dup.TFOG\n","nw_real_ji = pd.DataFrame(0, columns=quantile, index=stresses+[\"union\"])\n","\n","for j in quantile:\n","\tunion_interactions = []\n","\tfor i in stresses:\n","\t\telnet_no_dup = pd.read_csv(elnet_dir + i + \"/elnet_in_agris.txt\", header=0, index_col=None, sep=\"\\t\")\n","\t\t# R^2 > 0.8 cutoff\n","\t\telnet_no_dup = elnet_no_dup[elnet_no_dup.relvar > 0.8]\n","\t\telnet_no_dup[\"target_TF\"] = elnet_no_dup.predicted_GO + \"_\" + elnet_no_dup.Gene_GO\n","\t\t# Get coefficient cutoff for corresponding quantile\n","\t\tcutoff = abs(elnet_no_dup[\"rel.coeff\"]).quantile(j)\n","\t\telnet_no_dup = elnet_no_dup[abs(elnet_no_dup[\"rel.coeff\"]) >= cutoff] # change\n","\t\tif i != stresses[-1]:\n","\t\t\tunion_interactions.extend(elnet_no_dup.target_TF)\n","\t\t# Calculations\n","\t\tintersect = list(set(no_dup.target_TF.to_list()) & set(elnet_no_dup.target_TF.to_list()))\n","\t\tunion = list(set(no_dup.target_TF.to_list()) | set(elnet_no_dup.target_TF.to_list()))\n","\t\tnw_real_ji.loc[i, j] = len(intersect)/len(union)\n","\n","\t# for union network filtered at R^2 0.8 and corresponding quantiles\n","\tintersect = list(set(no_dup.target_TF.to_list()) & set(union_interactions))\n","\tunion = list(set(no_dup.target_TF.to_list()) | set(union_interactions))\n","\tnw_real_ji.loc[\"union\", j] = len(intersect)/len(union)\n","\n","# normal jaccard (GRN&AGRIS/GRN)\n","#nw_real_ji.plot(title=\"Overlap of AGRIS with GRN at various coefficient cutoffs (quantile)\",\n","#\t\t\t\t\t\t  xlabel = \"Networks\", ylabel = \"Jaccard Index\", style='.-').legend(bbox_to_anchor=(1, 1))\n","#plt.xticks(ticks = [i for i in range(9)] ,labels = [\"Light\", \"Dark\", \"Heat\", \"Cold\", \"Salt\", \"Mannitol\", \"Nitrogen deficiency\", \"All\", \"Union\"], rotation = 'vertical')\n","#plt.savefig(dir_path + 'figures/AGRIS_overlap_chart_nJI_quantile0.png', dpi = 600, bbox_inches='tight')\n","nw_real_ji.to_csv(dir_path + 'prep_files/all_quantile_nJI_quantile0.txt', sep=\"\\t\")\n","\n","# Just some summary stats\n","best_quantile_nJI = nw_real_ji.idxmax(axis=1)\n","Counter(nw_real_ji.idxmax(axis=1))\n","best_quantile_nJI.to_csv(dir_path + 'prep_files/best_quantile_nJI_quantile_0.txt', sep=\"\\t\", header=None)\n","\n","# Proceeding with custom cutoffs\n","best_union_interactions_nJI = []\n","best_quantile_overlap_nJI = nw_real_ji.iloc[:-1,:].max(axis=1).to_list()\n","union_df_nJI = []\n","\n","for i in stresses:\n","\telnet_no_dup = pd.read_csv(elnet_dir + i + \"/elnet_in_agris.txt\", header=0, index_col=None, sep=\"\\t\")\n","\t# R^2 > 0.8 cutoff\n","\telnet_no_dup = elnet_no_dup[elnet_no_dup.relvar > 0.8]\n","\telnet_no_dup[\"target_TF\"] = elnet_no_dup.predicted_GO + \"_\" + elnet_no_dup.Gene_GO\n","\t# Get coefficient cutoff for corresponding quantile\n","\tcutoff = abs(elnet_no_dup[\"rel.coeff\"]).quantile(best_quantile_nJI[i])\n","\telnet_no_dup = elnet_no_dup[abs(elnet_no_dup[\"rel.coeff\"]) >= cutoff]\n","\telnet_no_dup[\"nw_source\"] = [i for x in range(len(elnet_no_dup))]\n","\tif i != stresses[-1]:\n","\t\tbest_union_interactions_nJI.extend(elnet_no_dup.target_TF)\n","\t\tunion_df_nJI.append(elnet_no_dup)\n","intersect = list(set(no_dup.target_TF.to_list()) & set(best_union_interactions_nJI))\n","best_quantile_overlap_nJI.append(len(intersect)/len(set(best_union_interactions_nJI) | set(no_dup.target_TF.to_list()))) # [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\", \"union\"]\n","\n","new_networks = [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\", \"union\"]\n","with open(dir_path + 'prep_files/overlap_best_quantile_0.txt', \"w+\") as bf:\n","\tbf.write(\"\\t\" + \"\\t\".join(new_networks) + \"\\n\")\n","\tbf.write(\"quantile\\t\" + \"\\t\".join([str(x) for x in best_quantile[:-1]]) + \"\\tna\\n\")\n","\tbf.write(\"custom_ratio\\t\" + \"\\t\".join([str(x) for x in best_quantile_overlap]) + \"\\n\")\n","\tbf.write(\"quantile\\t\" + \"\\t\".join([str(x) for x in best_quantile_nJI[:-1]]) + \"\\tna\\n\")\n","\tbf.write(\"JI_ratio\\t\" + \"\\t\".join([str(x) for x in best_quantile_overlap_nJI]) + \"\\n\")\n","\n","# normal JI ratio\n","union_all_nJI = pd.concat(union_df_nJI)\n","union_all_nJI.to_csv(elnet_dir + '/union_nw_nJI_quantile0.txt', sep=\"\\t\", index=False)"],"metadata":{"id":"fijXwIgYyHmX"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["\"\"\"\n","Supp. Fig. S9B -- Part 2\n","\"\"\"\n","\n","import pandas as pd\n","import random\n","import seaborn as sns\n","from collections import defaultdict\n","import matplotlib.pyplot as plt\n","\n","agris_path = dir_path + 'prep_files/AtRegNet_confirmed_OG.txt'\n","elnet_dir = dir_path + 'prep_files/Mpo_GRN_models/elnet/'\n","elnet_name = '/elnet_in_agris.txt'\n","union_path = elnet_dir + 'union_nw_quantile0.txt'\n","union_path_nJI = elnet_dir + 'union_nw_nJI_quantile0.txt'\n","\n","# initialise agris\n","agris = pd.read_csv(agris_path, header=0, index_col=None, sep=\"\\t\")\n","agris_OG_nodup = agris.drop_duplicates(subset=[\"TFOG\", \"TargetOG\"])[[\"TargetOG\", \"TFOG\"]]\n","\n","agris_TFOG = agris_OG_nodup.TFOG.to_list()\n","agris_TargetOG = agris_OG_nodup.TargetOG.to_list()\n","agris_targetTF_OG = agris_OG_nodup.TargetOG +\"_\" + agris_OG_nodup.TFOG.to_list()\n","\n","\n","new_networks = [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\", \"union\"]\n","union_networks = [\"union2\", \"union3\", \"union4\", \"union5\", \"union6\", \"union7\"]\n","all_networks = [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\", \"all_cv5\", \"union1\", \"union2\", \"union3\", \"union4\", \"union5\", \"union6\", \"union7\"]\n","quantile_path = dir_path + 'prep_files/overlap_best_quantile_0.txt'\n","with open(quantile_path, \"r\") as q:\n","\tqcon = q.readlines()\n","\tbest_quantile = qcon[1].strip(\"\\n\").split(\"\\t\")[1:-1]\n","\texpected_ratio = [float(x) for x in qcon[2].strip(\"\\n\").split(\"\\t\")[1:]]\n","\tnJI_quantile = qcon[3].strip(\"\\n\").split(\"\\t\")[1:-1]\n","\texpected_JI = [float(x) for x in qcon[4].strip(\"\\n\").split(\"\\t\")[1:]]\n","\n","def get_grn_OG(nw_path, nw, index, cut_ref):\n","\telnet_no_dup = pd.read_csv(nw_path, header=0, index_col=None, sep=\"\\t\")\n","\tif nw != \"union\":\n","\t\t# R^2 > 0.8 cutoff\n","\t\telnet_no_dup = elnet_no_dup[elnet_no_dup.relvar > 0.8]\n","\t\telnet_no_dup[\"target_TF\"] = elnet_no_dup.predicted_GO + \"_\" + elnet_no_dup.Gene_GO\n","\t\t# Get coefficient cutoff for corresponding quantile\n","\t\tcutoff = abs(elnet_no_dup[\"rel.coeff\"]).quantile(float(cut_ref[index]))\n","\t\telnet_no_dup = elnet_no_dup[abs(elnet_no_dup[\"rel.coeff\"]) >= cutoff]\n","\treturn(list(set(elnet_no_dup.target_TF)))\n","\n","def get_grn_int(nw_path, nw, index, cut_ref):\n","\telnet_no_dup = pd.read_csv(nw_path, header=0, index_col=None, sep=\"\\t\")\n","\tif nw != \"union\":\n","\t\t# R^2 > 0.8 cutoff\n","\t\telnet_no_dup = elnet_no_dup[elnet_no_dup.relvar > 0.8]\n","\t\tcutoff = abs(elnet_no_dup[\"rel.coeff\"]).quantile(float(cut_ref[index]))\n","\t\telnet_no_dup = elnet_no_dup[abs(elnet_no_dup[\"rel.coeff\"]) >= cutoff]\n","\t\tgene_TF = elnet_no_dup[\"predicted\"] + \"_\" + elnet_no_dup[\"Gene.ID\"]\n","\telse:\n","\t\tunique_gene_pairs_df = elnet_no_dup.drop_duplicates(subset=[\"predicted\", \"Gene.ID\"])[[\"predicted\", \"Gene.ID\"]]\n","\t\tgene_TF = unique_gene_pairs_df[\"predicted\"] + \"_\" + unique_gene_pairs_df[\"Gene.ID\"]\n","\treturn(gene_TF.to_list())\n","\n","def get_ucut_OG(nw_path, nw):\n","\t# for union netork, presence in X network\n","\telnet_no_dup = pd.read_csv(nw_path, header=0, index_col=None, sep=\"\\t\")\n","\t# Get coefficient cutoff for corresponding quantile \n","\tcutoff = int(nw[-1])\n","\tunique_gene_pairs_df = elnet_no_dup.drop_duplicates(subset=[\"predicted\", \"Gene.ID\"])[[\"predicted\", \"Gene.ID\"]]\n","\tunique_gene_pairs = [[row[\"predicted\"], row[\"Gene.ID\"]] for idx, row in unique_gene_pairs_df.iterrows()]\n","\tOGpairs = []\n","\tfor pair in unique_gene_pairs:\n","\t\tsubset = elnet_no_dup[(elnet_no_dup[\"predicted\"] == pair[0]) & (elnet_no_dup[\"Gene.ID\"] == pair[1])]\n","\t\tif len(subset) >= cutoff:\n","\t\t\tOGpairs.append(subset.target_TF.to_list()[0])\n","\treturn(list(set(OGpairs)))\n","\n","def get_ucut_int(nw_path, nw):\n","\telnet_no_dup = pd.read_csv(nw_path, header=0, index_col=None, sep=\"\\t\")\n","\tcutoff = int(nw[-1])\n","\tunique_gene_pairs_df = elnet_no_dup.drop_duplicates(subset=[\"predicted\", \"Gene.ID\"])[[\"predicted\", \"Gene.ID\"]] #target_TF\n","\tunique_gene_pairs = [[row[\"predicted\"], row[\"Gene.ID\"]] for idx, row in unique_gene_pairs_df.iterrows()]\n","\tintpairs = []\n","\tfor pair in unique_gene_pairs:\n","\t\tsubset = elnet_no_dup[(elnet_no_dup[\"predicted\"] == pair[0]) & (elnet_no_dup[\"Gene.ID\"] == pair[1])]\n","\t\tif len(subset) >= cutoff:\n","\t\t\tintpairs.append(pair[0] + \"_\" + pair[1])\n","\treturn(intpairs)\n","\n","\n","# Permutation test\n","nw_OG_int_nJI = []\n","nw_int_nJI = []\n","for idx, network in enumerate(new_networks):\n","\tif network == \"union\":\n","\t\tfpath = union_path_nJI\n","\telse:\n","\t\tfpath = elnet_dir + network + elnet_name\n","\tnw_OG_int_nJI.append(get_grn_OG(fpath, network, idx, nJI_quantile))\n","\tnw_int_nJI.append(get_grn_int(fpath, network, idx, nJI_quantile))\n","\n","for idx, network in enumerate(union_networks):\n","\tfpath = union_path_nJI\n","\tOG_pairs = get_ucut_OG(fpath, network)\n","\tnw_OG_int_nJI.append(OG_pairs)\n","\tnw_int_nJI.append(get_ucut_int(fpath, network))\n","\texpected_JI.append(len(set(agris_targetTF_OG) & set(OG_pairs))/len(set(OG_pairs) | set(agris_targetTF_OG)))\n","\n","dicto2 = defaultdict(list)\n","for i in range(1000):\n","\trandom.shuffle(agris_TFOG)\n","\tint_pairs = [x + '_' + agris_TFOG[j] for j, x in enumerate(agris_TargetOG)]\n","\tfor k, nw in enumerate(all_networks):\n","\t\tintersection = len(set(int_pairs) & set(nw_OG_int_nJI[k]))\n","\t\tdicto2[\"network\"].append(nw)\n","\t\tdicto2[\"ratio\"].append(intersection/len(set(nw_OG_int_nJI[k]) | set(int_pairs)))\n","ratio_coll_nJI = pd.DataFrame.from_dict(dicto2)\n","pval_nJI = [sum(ratio_coll_nJI[ratio_coll_nJI.network == x].ratio >= expected_JI[i])/1000 for i, x in enumerate(all_networks)]\n","sig_nw = [x for i, x in enumerate(all_networks) if pval_nJI[i] < 0.05]\n","sig_height = [ratio_coll_nJI[ratio_coll_nJI.network == all_networks[i]].ratio.max() + 0.001 for i, x in enumerate(pval_nJI) if x < 0.05]\n","\n","# Plot result of permutuation test\n","fig, ax = plt.subplots()\n","sns.violinplot(x = \"network\", y = \"ratio\", data = ratio_coll_nJI, ax=ax, color=\"silver\", scale='width')\n","ax.set_xticklabels([\"Light\", \"Dark\", \"Heat\", \"Cold\", \"Salt\", \"Mannitol\", \"Nitrogen deficiency\", \"All\", \"Union (1)\", \"Union (2)\", \"Union (3)\", \"Union (4)\", \"Union (5)\", \"Union (6)\", \"Union (7)\"], rotation=90)\n","ax.set_xlabel(\"Network\")\n","ax.set_ylabel(\"Jaccard Index\")\n","ax.plot(all_networks, expected_JI, 'ko')\n","ax.plot(sig_nw, sig_height, 'k*')\n","plt.savefig(dir_path + 'figures/FigS9B.png', dpi = 600, bbox_inches='tight')"],"metadata":{"id":"bziEYoCgyUId"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["#### 2.2.4 Finalising the gene regulatory network; Fig 3A, visualised in Cytoscape"],"metadata":{"id":"uG1L2rpJz-yv"}},{"cell_type":"code","source":["# Get union network (including those not in AGRIS), cutoff R2 > 0.8\n","elnet_name = '/elnet.txt'\n","\n","u_stresses = [\"L\", \"D\", \"H\", \"C\", \"S\",\"M\", \"N\"]\n"," \n","def read_nw(nwpath, nw):\n","\tnwdf = pd.read_csv(nwpath, header=0, index_col=0, sep=\"\\t\")\n","\tnwdf_r2cf = nwdf[nwdf.relvar > 0.8]\n","\tnwdf_r2cf[\"network\"] = nw\n","\treturn(nwdf_r2cf)\n","\n","nw_dfs = [read_nw(elnet_dir + x + elnet_name, x) for x in u_stresses]\n","# concatenate union\n","union_raw = pd.concat(nw_dfs)\n","union_raw.to_csv(elnet_dir + \"union_0.8_ignoreAGRIS_full.txt\", index=False, sep=\"\\t\")"],"metadata":{"id":"T3Erx3UV0QPv"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Aim: To retrieve top TF (highest absolute coefficient) for each gene among all TFs\n","# Fig 3A, visualised in Cytoscape\n","# get union network (including not in AGRIS)\n","\n","\n","# comment next 2 lines if running for the first time\n","union_raw_p = elnet_dir + 'union_0.8_ignoreAGRIS_full.txt'\n","union_raw = pd.read_csv(union_raw_p, header=0, sep=\"\\t\")\n","\n","unique_targets = list(set(union_raw.predicted.to_list())) # 5878\n","top_int_df = pd.DataFrame(columns = [\"predicted\", \"Gene.ID\", \"coeffs\", \"present_in\", \"top_coeff\", \"top_stress\", \"TF_stat\"])\n","len_ut = len(unique_targets)\n","\n","def get_top_TF(target):\n","\tsubset = union_raw[union_raw.predicted == target]\n","\tsorted_subset = subset.reindex(subset[\"rel.coeff\"].abs().sort_values(ascending=False).index)\n","\ttop_TF = sorted_subset.iloc[0,:][\"Gene.ID\"]\n","\ttop_coeff = sorted_subset.iloc[0,:][\"rel.coeff\"]\n","\ttop_nw = sorted_subset.iloc[0,:][\"network\"]\n","\tTF_subset = sorted_subset[sorted_subset[\"Gene.ID\"] == top_TF]\n","\tnw_list = TF_subset.network.to_list()\n","\tcoeff_list = TF_subset[\"rel.coeff\"].to_list()\n","\tTF_stat_list = [\"P\" if x > 0 else \"N\" for x in coeff_list]\n","\tstat_sum = set(TF_stat_list)\n","\tif len(stat_sum) == 1:\n","\t\tif TF_stat_list[0] == \"P\":\n","\t\t\tTF_stat = \"Activator\"\n","\t\telif TF_stat_list[0] == \"N\":\n","\t\t\tTF_stat = \"Repressor\"\n","\telse:\n","\t\tTF_stat = \"Ambiguous\"\n","\ttop_int_df.loc[unique_targets.index(target)] = [target, top_TF, coeff_list, nw_list, top_coeff, top_nw, TF_stat]\n","\n","for i, target in enumerate(unique_targets):\n","\tget_top_TF(target)\n","\tif i % 1000 == 0:\n","\t\tprint(\"Processed: \" + str(i) + \"/\" + str(len_ut))\n","\n","top_int_df.to_csv(elnet_dir + \"union_0.8_ignoreAGRIS_topTF.txt\", index=False, sep=\"\\t\")\n","\n","from collections import Counter\n","import matplotlib.pyplot as plt\n","stat_count = Counter(top_int_df.TF_stat.to_list())\n","plt.pie(stat_count.values(), labels = stat_count.keys(), autopct = lambda x: '{:.0f}'.format(x*sum(stat_count.values())/100))\n","nw_count = Counter(top_int_df.top_stress.to_list()) # Counter({'M': 536, 'C': 816, 'S': 961, 'D': 750, 'H': 816, 'N': 1086, 'L': 913})\n","plt.pie(nw_count.values(), labels = nw_count.keys(), autopct = lambda x: '{:.0f}'.format(x*sum(nw_count.values())/100))"],"metadata":{"id":"_OohCQoNuvCw"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["#### 2.2.5. TF-only GRN, Fig 4A (visualised in cytoscape) and Supp. Fig 11"],"metadata":{"id":"xAAY08O09Xan"}},{"cell_type":"code","source":["import pandas as pd\n","import matplotlib.pyplot as plt\n","\n","# Paths\n","elnet_path = elnet_dir + '/union_0.8_ignoreAGRIS_full.txt'\n","tf_path = dir_path + 'prep_files/PlantTFDB_Mpov5r1_prediction_plusTFDB.txt'\n","\n","elnet_df = pd.read_csv(elnet_path, header=0, sep=\"\\t\")\n","tf_list = [x.split(\"\\t\")[0] for x in open(tf_path, \"r\").readlines()]\n","\n","elnet_tf_df = elnet_df[(elnet_df.predicted.isin(tf_list)) & (elnet_df[\"Gene.ID\"].isin(tf_list))] # network containing TFs only\n","\n","def overall_spec(df_row):\n","    # retrieve status from specificity file\n","\tstatus = [\"UP\" if x > 0 else \"DOWN\" for x in list(df_row) if x != 0]\n","\tstatus = list(set(status))\n","\tif status == [\"NC\"]:\n","\t\treturn \"NS\"\n","\telif status == [\"UP\"]:\n","\t\treturn \"UP\"\n","\telif status == [\"DOWN\"]:\n","\t\treturn \"DOWN\"\n","\telse:\n","\t\treturn \"MIXED\"\n","\n","# specificity file, spec > 0.7\t\n","desc_path = elnet_dir + 'TF_Scond_specificity_07.txt'\n","desc_file = pd.read_csv(desc_path, header=0, sep=\"\\t\")\n","desc_file.columns = [\"desc\"] + desc_file.columns.to_list()[1:]\n","desc_file[\"gene\"] = desc_file.desc.apply(lambda x: x.split(\" \")[0])\n","desc_file[\"overall_spec\"] = desc_file.iloc[:,1:8].apply(lambda x: overall_spec(x), axis=1)\n","# annotation file for network when viewed in cytoscape\n","desc_file.to_csv(dir_path + 'prep_files/network_anno.txt', index=False, sep=\"\\t\")\n","\n","# for each gene, get the top DEG status in stress X containing network\n","elnet_tf_df[\"predicted_stat\"] = elnet_tf_df.apply(lambda x: desc_file[desc_file.gene == x.predicted].overall_spec.values[0] if len(desc_file[desc_file.gene == x.predicted].overall_spec.values)!= 0 else \"NS\", axis=1)\n","elnet_tf_df[\"geneid_stat\"] = elnet_tf_df.apply(lambda x: desc_file[desc_file.gene == x[\"Gene.ID\"]].overall_spec.values[0] if len(desc_file[desc_file.gene == x[\"Gene.ID\"]].overall_spec.values)!= 0 else \"NS\", axis=1)\n","\n","def check_outcome(pstat, gstat, relcoef):\n","\tif pstat == \"NS\" and gstat == \"NS\":\n","\t\t\treturn \"unexpected\"\n","\telif pstat == \"NS\" or gstat == \"NS\":\n","\t\t\treturn \"unexpected\"\n","\telif pstat == \"MIXED\" and gstat == \"MIXED\":\n","\t\t\treturn \"unexpected\"\n","\telif pstat == \"MIXED\" or gstat == \"MIXED\":\n","\t\t\treturn \"unexpected\"\n","\telif len(set([pstat, gstat])) == 1:\n","\t\tif relcoef > 0:\n","\t\t\treturn \"expected\"\n","\t\telse:\n","\t\t\treturn \"unexpected\"\n","\telse:\n","\t\tif relcoef < 0:\n","\t\t\treturn \"expected\"\n","\t\telse:\n","\t\t\treturn \"unexpected\"\n","elnet_tf_df[\"grn_tally\"] = elnet_tf_df.apply(lambda x: check_outcome(x.predicted_stat, x.geneid_stat, x[\"rel.coeff\"]), axis=1)\n","\n","# filter for unique target-TF pairs with top relative coefficient\n","gene_pairs = [list(x) for i, x in elnet_tf_df[[\"predicted\", \"Gene.ID\"]].drop_duplicates().iterrows()]\n","elnet_tf_unique = pd.DataFrame(columns=elnet_tf_df.columns.to_list())\n","for g in gene_pairs:\n","\tsub = elnet_tf_df[(elnet_tf_df.predicted == g[0]) & (elnet_tf_df[\"Gene.ID\"] == g[1])]\n","\telnet_tf_unique.loc[len(elnet_tf_unique)] = sub.loc[abs(sub[\"rel.coeff\"]).idxmax()].to_list()\n","\n","def plot_stats(df, ytype, start=0, end=55):\n","    # plots number of expected and unexpected edges for each coefficient cutoff 0 - 0.54, step 0.01\n","\tcutoff_count = [df[abs(df[\"rel.coeff\"]) > i/100].grn_tally.value_counts() for i in range(start, end)]\n","\tcutoff_ratio = [x.expected/sum(x) for x in cutoff_count]\n","\tif ytype == \"edges\":\n","\t\tcutoff_size = [sum(x) for x in cutoff_count]\n","\t\tylab = \"Number of edges\"\n","\telif ytype == \"nodes\":\n","\t\tcutoff_size = [len(list(set(x.predicted.to_list() + x[\"Gene.ID\"].to_list()))) for x in [df[abs(df[\"rel.coeff\"]) > i/100] for i in range(start, end)]]\n","\t\tylab = \"Number of nodes\"\n","\t\n","\tfig, ax1 = plt.subplots()\n","\tcolor = 'royalblue'\n","\tax1.set_xlabel('Absolute relative coefficient cutoff')\n","\tax1.set_ylabel('Ratio of expected versus total interactions', color=color)\n","\tax1.plot([i/100 for i in range(start,end)], cutoff_ratio, color=color)\n","\tax1.tick_params(axis='y', labelcolor=color)\n","\t\n","\tax2 = ax1.twinx()\n","\tcolor = 'firebrick'\n","\tax2.set_ylabel(ylab, color=color)  # we already labellled the x-label in ax1\n","\tax2.plot([i/100 for i in range(start,end)], cutoff_size, color=color)\n","\tax2.tick_params(axis='y', labelcolor=color)\n","\t\n","\tfig.tight_layout()\n","\tplt.savefig(dir_path + 'figures/' + 'FigS11_{}.png'.format(ytype), dpi=600)\n","\tplt.show()\n","\treturn[[i/100 for i in range(start, end)], cutoff_ratio, cutoff_size]\n","\n","node_details = plot_stats(elnet_tf_unique, \"nodes\")\n","edge_details = plot_stats(elnet_tf_unique, \"edges\")\n","\n","# Find optimal cutoff and output network to file\n","optimal_cutoff = node_details[0][node_details[1].index(max(node_details[1]))] #0.22\n","filtered_elnet = elnet_tf_unique[abs(elnet_tf_unique[\"rel.coeff\"]) > optimal_cutoff]\n","filtered_elnet[\"abs_relcoeff\"] = filtered_elnet[\"rel.coeff\"].apply(abs)\n","filtered_elnet.to_csv(dir_path + 'prep_files/network_cf22_spec.txt', index=False, sep=\"\\t\")"],"metadata":{"id":"8pZr30ab9cD6"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### 2.3 Extract Arabidopsis GO and TFs"],"metadata":{"id":"rhp4-bP6Ju_X"}},{"cell_type":"code","source":["\"\"\"\n","Flatten GO file to 1 line per gene\n","\"\"\"\n","\n","import pandas as pd\n","\n","GOpath = dir_path + 'prep_files/ATH_GO_GOSLIM.txt'\n","TFpath = dir_path + 'prep_files/Ath_TF_list.txt'\n","opath = dir_path + 'prep_files/'\n","\n","GO_df = pd.read_csv(GOpath, skiprows = 4, header=None, sep=\"\\t\")\n","GO_df.columns = ['locus_name', 'TAIR_accession', 'object_name',\n","\t\t\t\t 'relationship_type', 'GO_term', 'GO_ID', 'TAIR_Keyword_ID',\n","\t\t\t\t 'Aspect', 'GOslim_term', 'Evidence_code', 'Evidence_description',\n","\t\t\t\t 'Evidence_with', 'Reference','Annotator', 'Date_annotated']\n","TF_df = pd.read_csv(TFpath, header=0, sep=\"\\t\")\n","TF_genes = list(set(TF_df.Gene_ID))\n","all_genes = GO_df.locus_name.unique().tolist()\n","\n","# define codes relating to experimental evidence and high throughput experiments\n","exp_codes = ['EXP', 'IDA', 'IPI', 'IMP', 'IGI', 'IEP']\n","exp_codes_htp = ['HTP', 'HDA', 'HMP', 'HGI', 'HEP']\n","GO_df_EXP = GO_df[GO_df.Evidence_code.isin(exp_codes)]\n","EXP_genes = list(GO_df_EXP.locus_name.unique())\n","EXP_TF = list(set(EXP_genes) & set(TF_genes))\n","GO_df_HTP = GO_df[GO_df.Evidence_code.isin(exp_codes_htp)]\n","HTP_genes = [x for x in list(GO_df_HTP.locus_name.unique()) if x not in EXP_TF]\n","HTP_TF = list(set(HTP_genes) & set(TF_genes))\n","EXP_HTP_TF = list(set(EXP_TF) | set(HTP_TF))\n","\n","dicto = {'locus': [], 'MolFunc': [], 'BioProc': [], 'type': []}\n","for i, gene in enumerate(TF_genes):\n","\tif i+1 % 100 == 0:\n","\t\tprint(\"{}/1000 genes\".format(i+1))\n","\tMolF = list(GO_df[(GO_df.locus_name == gene) & (GO_df.Aspect == \"F\")].GO_term.unique())\n","\tBioP = list(GO_df[(GO_df.locus_name == gene) & (GO_df.Aspect == \"P\")].GO_term.unique())\n","\tdicto['locus'].append(gene)\n","\tdicto['MolFunc'].append(MolF)\n","\tdicto['BioProc'].append(BioP)\n","\t\n","\tif gene in EXP_TF:\n","\t\tdicto['type'].append('EXP')\n","\telif gene in HTP_TF:\n","\t\tdicto['type'].append('HTP')\n","\telse:\n","\t\tdicto['type'].append('OTHERS')\n","\t\t\n","EXP_TF_aspect_df = pd.DataFrame.from_dict(dicto)\n","EXP_TF_aspect_df.to_csv(opath + 'Ath_TFall_GOanno.txt',\n","\t\t\t\t\t\tsep=\"\\t\", index=False)"],"metadata":{"id":"_EthsbGMJ0T-"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["\"\"\"\n","Map Mpo differentially expressed TFs to Ath\n","\"\"\"\n","\n","# get differentially expressed Mpo TF\n","mpo_tf_path = dir_path + \"prep_files/PlantTFDB_Mpov5r1_prediction_plusTFDB.txt\"\n","mpo_exp_path = dir_path + \"prep_files/Mpo_GRN_models/all_stress_mt5_nodum.tsv\"\n","\n","mpo_tf_list = pd.read_csv(mpo_tf_path, sep=\"\\t\").iloc[:,0]\n","mpo_exp_tf_list = pd.read_csv(mpo_exp_path, header=0, sep=\"\\t\").iloc[:,0]\n","\n","mpo_nw_tf = list(set(mpo_tf_list) & set(mpo_exp_tf_list))\n","\n","# orthogroup file\n","OF_path = dir_path + 'prep_files/Orthogroups.txt'\n","gene_OF = {}\n","\n","with open(OF_path, 'r') as OF_file:\n","\tfor line in OF_file:\n","\t\tcontent = line.strip(\"\\n\").split(\": \")\n","\t\tog = content[0]\n","\t\tgene_list = content[1].split(\" \")\n","\t\tfor gene in gene_list:\n","\t\t\tgene_OF[gene] = og\n","\n","# Ath TF GO file (flat)\n","ath_go_path = dir_path + 'prep_files/Ath_TFall_GOanno.txt'\n","ath_go_df = pd.read_csv(ath_go_path, header=0, index_col=0, sep=\"\\t\")\n","ath_go_df[\"OG\"] = [gene_OF[x] for x in ath_go_df.index.to_list()]\n","\n","# conversion\n","mpo_nw_tf_og = [gene_OF[x] for x in mpo_nw_tf]\n","\n","unique_og = list(set(mpo_nw_tf_og))\n","ath_genes = {og: ath_go_df[ath_go_df.OG == og].index.to_list() for og in unique_og}\n","ath_desc = {og: [list(x)[:-1] for i, x in ath_go_df[ath_go_df.OG == og].iterrows()] for og in unique_og}\n","\n","compilation = pd.DataFrame.from_dict({\"Mpo_genes\": mpo_nw_tf, \"OG\": mpo_nw_tf_og})\n","compilation[\"Ath_genes\"] = compilation.OG.apply(lambda x: ath_genes[x])\n","compilation[\"Ath_desc\"] = compilation.OG.apply(lambda x: ath_desc[x])\n","compilation_nona_expOnly = compilation[[\"EXP\" in [z for y in x for z in y] for x in compilation.Ath_desc]]\n","compilation_nona_expOnly.to_csv(dir_path + 'prep_files/MpoDEGTF_AthAllTF_noNA_exp.txt', index=False, sep=\"\\t\")"],"metadata":{"id":"d4NPWKbTKe87"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"V9HG7DByg6Pg"},"source":["# 3. Analysis and plotting"]},{"cell_type":"markdown","metadata":{"id":"4qA5pIrH7VZW"},"source":["### Figure 1 & Supp. Fig 1: Measurements and Student's t-test"]},{"cell_type":"code","metadata":{"id":"DXpM_NKa743D"},"source":["# adated from measurements_forsupp.py\n","wdir = dir_path + 'prep_files/'\n","odir = dir_path + 'figures/'\n","odir_safe = dir_path_safe + 'figures/'\n","if not os.path.exists(odir):\n","    !mkdir $odir_safe\n","infile = 'phase1n2_measurements_nooutliers.txt'\n","\n","measurements = pd.read_csv(wdir + infile, sep='\\t')\n","# Single letter to full single stress description\n","singled = {'C':'Cold',\n","\t\t   'H':'Heat',\n","\t\t   'S':'Salt',\n","\t\t   'M':'Mannitol',\n","\t\t   'L':'Light',\n","\t\t   'D':'Dark',\n","\t\t   'N':'Nitrogen'\n","\t\t   }\n","\n","# For plotting all controls\n","areatype = ['Parea', 'Earea']\n","titletype = ['15', '21']\n","for i in range(0,2):\n","\tarea = areatype[i]\n","\ttitle = titletype[i]\n","\tcontrol_m = measurements[measurements.Stress == 'Control'][[\"Batch\", area]].groupby('Batch', sort = False).mean()\n","\tcontrol_s = measurements[measurements.Stress == 'Control'][[\"Batch\", area]].groupby('Batch', sort = False).std()\n","\tcontrol_m.plot.bar(yerr=[list(control_s[area]), list(control_s[area])[::-1]], legend=False, title='Control (Day '+ title + ')', capsize=4)\n","\tplt.savefig(odir + 'Control_Day' + title + '.png', dpi = 600, bbox_inches='tight')\n","\tplt.show()\n","\n","# df with only single stress measurements\n","ss_meas = measurements[(measurements.Condition != 'None') & (measurements.Condition != 'mixed')]\n","# df with only crossed stress measurements\n","cs_meas = measurements[measurements.Condition == 'mixed']\n","# df with only controls\n","control_meas = measurements[measurements.Stress == 'Control']\n","\n","controltype = ['Stress']\n","controltitle = ['_merged']\n","\n","xaxislabel = {'Heat': 'Temperature (\\u00B0C)',\n","\t\t\t  'Cold': 'Temperature (\\u00B0C)',\n","\t\t\t  'Mannitol': 'Mannitol (mM)',\n","\t\t\t  'Salt': 'NaCl (mM)',\n","\t\t\t  'Light': 'Light intensity (\\u03bcEm\\u207b\\u00b2s\\u207b\\u00b9)',\n","\t\t\t  'Dark': 'Days',\n","\t\t\t  'Nitrogen': 'KNO\\u2083 (%)'}\n","\n","# Supp. Fig. 1\n","# t-test (control as b, following test, a)\n","tout = open(wdir + 'ttest.txt', 'w+') \n","from scipy import stats as st\n","\n","for ss in list(ss_meas.Stress.unique()):\n","\tfor i, c in enumerate(controltype):\n","\t\tcontrol_batches = ss_meas[ss_meas.Stress == ss].Batch.unique()\n","\t\tcontrol_mean = control_meas[control_meas.Batch.isin(control_batches)].groupby(c, sort = False).mean()\n","\t\tcontrol_std = control_meas[control_meas.Batch.isin(control_batches)].groupby(c, sort = False).std()\n","\t\tstress_mean = ss_meas[ss_meas.Stress == ss].groupby('Condition', sort=False).mean()\n","\t\tstress_std = ss_meas[ss_meas.Stress == ss].groupby('Condition', sort=False).std()\n","\t\t\n","\t\tif c == 'Stress': # t-test\n","\t\t\tcontrol_df = control_meas[control_meas.Batch.isin(control_batches)][['Parea', 'Earea']]\n","\t\t\tstress_conds = ss_meas[ss_meas.Stress == ss].Condition.unique()\n","\t\t\tfor k, a in enumerate(areatype):\n","\t\t\t\tfor scond in stress_conds:\n","\t\t\t\t\tstress_df = ss_meas[(ss_meas.Stress == ss) & (ss_meas.Condition == scond)]\n","\t\t\t\t\ttstat, pval = st.ttest_ind(stress_df[a], control_df[a])\n","\t\t\t\t\ttout.write(('\\t').join(['Day '+ titletype[k], ss + '_' + scond, 'control_merged', str(tstat), str(pval)]) + \"\\n\")\n","\t\t\n","\t\tlabels = control_mean.index.to_list() + stress_mean.index.to_list()\n","\t\tfor j, a in enumerate(areatype): # Day 15 or 21 area\n","\t\t\tcoll_mean = list(control_mean[a]) + list(stress_mean[a])\n","\t\t\tcoll_std = list(control_std[a]) + list(stress_std[a])\n","\t\t\tplt.bar(labels, coll_mean, yerr = coll_std, capsize=4)\n","\t\t\tplt.title(ss + ' (Day ' + titletype[j] + ')')\n","\t\t\tplt.xlabel(xaxislabel[ss]) \n","\t\t\tplt.ylabel('Area (mm\\u00b2)')\n","\t\t\tplt.savefig(odir + ss + '_Day' + titletype[j] + controltitle[i] + '.png', dpi = 600, bbox_inches='tight')\n","\t\t\tplt.show()\n","\n","# cross_stress plot\n","for i, c in enumerate(controltype):\n","\tcs_control_batches = cs_meas.Batch.unique()\n","\tcs_control_mean = control_meas[control_meas.Batch.isin(cs_control_batches)].groupby(c, sort = False).mean()\n","\tcs_control_std = control_meas[control_meas.Batch.isin(cs_control_batches)].groupby(c, sort = False).std()\n","\tcs_stress_mean = cs_meas.groupby('Stress', sort=False).mean()\n","\tcs_stress_std = cs_meas.groupby('Stress', sort=False).std()\n","\tcs_labels = cs_control_mean.index.to_list() + cs_stress_mean.index.to_list()\n","\t\n","\tif c == 'Stress': #t-test\n","\t\tcontrol_df = control_meas[control_meas.Batch.isin(cs_control_batches)][['Parea', 'Earea']]\n","\t\tfor k, a in enumerate(areatype):\n","\t\t\tfor ss in list(cs_meas.Stress.unique()):\n","\t\t\t\tstress_df = cs_meas[(cs_meas.Stress == ss)]\n","\t\t\t\ttstat, pval = st.ttest_ind(stress_df[a], control_df[a])\n","\t\t\t\ttout.write(('\\t').join(['Day '+ titletype[k], ss + '_' + 'mixed', 'control_merged', str(tstat), str(pval)]) + \"\\n\")\n","\t\t\t\t\t\n","\tfor j, a in enumerate(areatype): # Day 15 or 21 area\n","\t\t\tcoll_mean = list(cs_control_mean[a]) + list(cs_stress_mean[a])\n","\t\t\tcoll_std = list(cs_control_std[a]) + list(cs_stress_std[a])\n","\t\t\tplt.bar(cs_labels, coll_mean, yerr = coll_std, capsize=4)\n","\t\t\tplt.title('Cross stress (Day ' + titletype[j] + ')')\n","\t\t\tplt.xticks(rotation=90)\n","\t\t\tplt.xlabel('Experiment') \n","\t\t\tplt.ylabel('Area (mm\\u00b2)')\n","\t\t\tplt.savefig(odir + 'Cross_stress_Day' + titletype[j] + controltitle[i] + '.png', dpi = 600, bbox_inches='tight')\n","\t\t\tplt.show()\n","\t\t\t\n","# single stress reps and cross stress (control - merged)\n","def ss_grab(stress, condition):\n","\t\"\"\"\n","\tSlice the relevant condition for \n","\tParameters\n","\t----------\n","\tstress : string\n","\t\tStress of interest.\n","\tcondition : string\n","\t\tCondition of interest.\n","\n","\tReturns\n","\t-------\n","\tsssub : dataframe\n","\t\tdf of single stress.\n","\n","\t\"\"\"\n","\tsssub = measurements[(measurements.Stress == stress) & (measurements.Condition == condition)]\n","\treturn sssub\n","\n","srep_keys = [['Cold', '3'],\n","\t\t\t ['Heat', '33'],\n","\t\t\t ['Salt', '40'],\n","\t\t\t ['Mannitol', '100'],\n","\t\t\t ['Light', '435'],\n","\t\t\t ['Dark', '3'],\n","\t\t\t ['Nitrogen', '0']]\n","srepdf = measurements[(measurements.Stress == 'Cold') & (measurements.Condition == '3')]\n","\n","for s, c in srep_keys[1:]:\n","\tsrepdf = pd.concat([srepdf, ss_grab(s, c)])\n","\t\n","s_cs_meas = pd.concat([srepdf, cs_meas])\n","\t\n","s_cs_control_batches = list(cs_meas.Batch.unique()) + list(srepdf.Batch.unique())\n","s_cs_control_mean = control_meas[control_meas.Batch.isin(s_cs_control_batches)].groupby('Stress', sort = False).mean()\n","s_cs_control_std = control_meas[control_meas.Batch.isin(s_cs_control_batches)].groupby('Stress', sort = False).std()\n","s_cs_stress_mean = s_cs_meas.groupby('Stress', sort=False).mean()\n","s_cs_stress_std = s_cs_meas.groupby('Stress', sort=False).std()\n","\n","s_cs_meas_label = [x + ' (' + x[0] + ')' if len(x) > 2 else x for x in s_cs_stress_mean.index]\n","s_cs_meas_label[s_cs_meas_label.index('Light (L)')] = 'High light (L)'\n","s_cs_meas_label[s_cs_meas_label.index('Dark (D)')] = 'Darkness (D)'\n","s_cs_labels = s_cs_control_mean.index.to_list() + s_cs_meas_label\n","\n","# Fig1\n","#t-test\n","control_df = control_meas[control_meas.Batch.isin(s_cs_control_batches)][['Parea', 'Earea']]\n","## cross-stress\n","for k, a in enumerate(areatype):\n","\tfor ss in list(cs_meas.Stress.unique()):\n","\t\tstress_df = cs_meas[(cs_meas.Stress == ss)]\n","\t\tfor singleS in ss:\n","\t\t\tsinglecontrol = srepdf[srepdf.Stress == singled[singleS]][a]\n","\t\t\ttstat, pval = st.ttest_ind(stress_df[a], singlecontrol)\n","\t\t\ttout.write(('\\t').join(['Day '+ titletype[k], ss + '_' + 'mixed', 'control_' + singled[singleS], str(tstat), str(pval)]) + \"\\n\")\n","## single stress\n","for k, a in enumerate(areatype):\n","\tfor ss in list(srepdf.Stress.unique()):\n","\t\tcond = srepdf[(srepdf.Stress == ss)].Condition.unique()[0]\n","\t\tstress_df = srepdf[(srepdf.Stress == ss)][a]\n","\t\ttstat, pval = st.ttest_ind(stress_df, control_df[a])\n","\t\ttout.write(('\\t').join(['Day '+ titletype[k], ss + '_' + cond, 'control', str(tstat), str(pval)]) + \"\\n\")\n","\t\t\n","tout.close()\t\n","# plotting\n","colour_seq = ['tomato'] +  ['mediumseagreen']*7 + ['cornflowerblue']*20\n","for j, a in enumerate(areatype): # Day 15 or 21 area\n","\t\tcoll_mean = list(s_cs_control_mean[a]) + list(s_cs_stress_mean[a])\n","\t\tcoll_std = list(s_cs_control_std[a]) + list(s_cs_stress_std[a])\n","\t\tplt.bar(s_cs_labels, coll_mean, yerr = coll_std, capsize=4, color = colour_seq)\n","\t\tplt.title('Area (Day ' + titletype[j] + ')')\n","\t\tplt.xticks(rotation=90)\n","\t\tplt.xlabel('Experiment') \n","\t\tplt.ylabel('Area (mm\\u00b2)')\n","\t\tplt.savefig(odir + 'fig1_Day' + titletype[j] + '.png', dpi = 600, bbox_inches='tight')\n","\t\tplt.show()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"RdmYRmdc7oA4"},"source":["### Supp. Fig 2: QC of RNA-seq data"]},{"cell_type":"code","metadata":{"id":"ZUjfF_uuYTGW"},"source":["# adapted from QC_scaled_updated.py\n","from sklearn.preprocessing import StandardScaler\n","from scipy.stats import pearsonr\n","\n","o_dir = dir_path + 'figures/'\n","sumdir = dir_path + 'summary_files/'\n","\n","expdesc = ['all_stress', 'diurnal_exp', 'single_stress']\n","targetexp = expdesc[0]\n","targetp = sumdir + targetexp + '.txt'\n","expmatp = dir_path + 'prep_files/' + targetexp + '.tsv'\n","exps = [x.split(\"\\t\")[0] for x in open(targetp, \"r\").readlines()]\n","labels = [x.strip().split(\"\\t\")[1] + '_' + x.split(\"\\t\")[0].split('_')[1] for x in open(targetp, \"r\").readlines()]\n","\n","df = pd.read_csv(expmatp, index_col = 0, sep = \"\\t\", header = 0)\n","df.columns = labels\n","\n","# Standard Scaling\n","scaled_features = StandardScaler().fit_transform(df.values)\n","df_scaled = pd.DataFrame(scaled_features, index = df.index, columns = df.columns)\n","\n","# plot cluster map\n","sns.set(font_scale=1.6)\n","\n","methods = \"average\"\n","\n","g1 = sns.clustermap(df_scaled.corr(),\n","\t\t\t\t method = methods,\n","\t\t\t\t figsize=(20,20),\n","\t\t\t\t xticklabels=True,\n","\t\t\t\t yticklabels=True)\n","plt.title(\"All stress (scaled): \" + methods)\n","plt.savefig(o_dir + \"SuppFig2\" + '.png')\n","\n","\n","# PCC of experiments\n","pcc_out = open(dir_path + \"prep_files/mpo/all_stress_PCC.txt\", \"w+\")\n","pcc_out.write(\"exp1\\texp2\\tpcc_val\\tp_value\\n\")\n","exps = list(df_scaled.columns)\n","for exp1 in range(len(exps)):\n","\tfor exp2 in range(exp1):\n","\t\tif exps[exp1].split(\"_\")[0] == exps[exp2].split(\"_\")[0]:\n","\t\t\tpcc_val, p_value = pearsonr(df_scaled[exps[exp1]], df_scaled[exps[exp2]])\n","\t\t\tpcc_out.write(exps[exp1] + \"\\t\" + exps[exp2] + \"\\t\" + str(pcc_val) + \"\\t\" + str(p_value) + \"\\n\")\n","pcc_out.close()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"cTKSAfIj7rcq"},"source":["### Supp. Fig 3: Volcano plots (DESeq2)"]},{"cell_type":"code","metadata":{"id":"5Cp2Q3EKrCD3"},"source":["# adapted from deseq_volcano.py\n","\n","wdir = dir_path + 'prep_files/mpo/deseq/'\n","odir = wdir + 'volcano/'\n","deseqouts = [x for x in os.listdir(wdir) if \"res.tsv\" in x] # controlD2controlH2_res.tsv\n","control = 'controlD2controlH2_res.tsv'\n","deseqouts.pop(deseqouts.index(control))\n","all_stress = list(set([x.split('control')[0] for i, x in enumerate(deseqouts)]))\n","all_stress.sort()\n","s_stress = [x for x in all_stress if len(x) == 1]\n","c_stress = [x for x in all_stress if len(x) == 2]\n","all_stress = s_stress + c_stress\n","\n","# plot control\n","controls = pd.read_csv(wdir + control,\n","\t\t\t  sep = \"\\t\", header = 0, index_col = 0)\n","sns.scatterplot(x = controls['log2FoldChange'],\n","\t\t\t\t y = -np.log10(controls[\"padj\"]),\n","\t\t\t\t #ax = axs[axcord[0][0], axcord[0][1]],\n","\t\t\t\t alpha = 0.2,\n","\t\t\t\t marker = '.',\n","\t\t\t\t legend = False,\n","\t\t\t\t edgecolor = \"none\",\n","\t\t\t\t hue = np.logical_and(abs(controls['log2FoldChange']) > 1,\n","\t\t\t\t\t\t  -np.log10(controls[\"padj\"]) > -np.log10(0.05)))\n","plt.title(\"control \" + control.split(\"control\")[1] + \" vs control H2\")\n","\n","# plot everything else\n","xlen = 5\n","ylen = math.ceil(len(deseqouts)/5)\n","fig, axs = plt.subplots(ylen, xlen, figsize=(30, 37.5), sharex='col', sharey='row')\n","#sns.set(font_scale=1.6)\n","axcord = []\n","for a in range(ylen):\n","\tfor b in range(xlen):\n","\t\taxcord.append([a, b])\n","\n","for i, z in enumerate(all_stress):\n","\tfiles = [x for x in deseqouts if x.startswith(z+'control')]\n","\tfiles.sort()\n","\tfileD2 = pd.read_csv(wdir + files[0],\n","\t\t\t\t\t  sep = \"\\t\", header = 0, index_col = 0)\n","\tfileH2 = pd.read_csv(wdir + files[1],\n","\t\t\t\t\t  sep = \"\\t\", header = 0, index_col = 0)\n","\tD2ax = int(((i*2)//10)*10 + ((i*2)%10)/2)\n","\tH2ax = int(D2ax + 5)\n","\t# Volcano plots\n","\t# against control D2\n","\tsns.scatterplot(x = fileD2['log2FoldChange'],\n","\t\t\t\t y = -np.log10(fileD2[\"padj\"]),\n","\t\t\t\t ax = axs[axcord[D2ax][0], axcord[D2ax][1]],\n","\t\t\t\t alpha = 0.2,\n","\t\t\t\t marker = '.',\n","\t\t\t\t legend = False,\n","\t\t\t\t edgecolor = \"none\",\n","\t\t\t\t hue = np.logical_and(abs(fileD2['log2FoldChange']) > 1,\n","\t\t\t\t\t\t  -np.log10(fileD2[\"padj\"]) > -np.log10(0.05)))\n","\taxs[axcord[D2ax][0], axcord[D2ax][1]].set_title(z + \" vs control D2\")\n","\t\n","\t# against control H2\n","\tsns.scatterplot(x= fileH2['log2FoldChange'],\n","\t\t\t\t y = -np.log10(fileH2[\"padj\"]),\n","\t\t\t\t ax = axs[axcord[H2ax][0], axcord[H2ax][1]],\n","\t\t\t\t alpha = 0.2,\n","\t\t\t\t marker = '.',\n","\t\t\t\t legend = False,\n","\t\t\t\t edgecolor = \"none\",\n","\t\t\t\t hue = np.logical_and(abs(fileH2['log2FoldChange']) > 1,\n","\t\t\t\t\t\t  -np.log10(fileH2[\"padj\"]) > -np.log10(0.05)))\n","\taxs[axcord[H2ax][0], axcord[H2ax][1]].set_title(z + \" vs control H2\")\n","\t\n","for ax in axs.flat:\n","    ax.set(xlabel='log2FoldChange', ylabel='-log10 padj')\n","\n","# Hide x labels and tick labels for top plots and y ticks for right plots.\n","for ax in axs.flat:\n","    ax.label_outer()\n","\n","plt.savefig(dir_path + \"figures/SuppFig3.png\", dpi = 600)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"OZ8K1xcl7dfZ"},"source":["### Supp. Fig 4: Comparison of DEGs between two controls"]},{"cell_type":"code","metadata":{"id":"-Bz6wm_ABLin"},"source":["wdir = mpo_path\n","deseqouts = [x for x in os.listdir(wdir) if \"resSig.tsv\" in x]\n","deseqouts.remove('controlD2controlH2_resSig.tsv')\n","\n","# create subplots\n","xlen = 4\n","ylen = math.ceil(len(deseqouts)/8)\n","figw = xlen * 4\n","figh = ylen * 2.5\n","\n","stress_list = ['H', 'C', 'HM', 'CM', 'M', 'CL', 'ML', 'L', 'HS', 'CS', 'SM', 'LS', 'S',\n","\t\t\t   'HN', 'CN', 'MN', 'NL', 'SN', 'N', 'HD', 'CD', 'MD', 'SD', 'ND', 'D']\n","f_axes = string.ascii_uppercase[:len(stress_list)]\n","axd = plt.figure(constrained_layout=True,\n","\t\t\t\t figsize=(figw, figh)).subplot_mosaic(\n","\t\t\t\t\t \"\"\"\n","\t\t\t\t\t A......\n","\t\t\t\t\t .B.....\n","\t\t\t\t\t CDE....\n","\t\t\t\t\t .FGH...\n","\t\t\t\t\t IJKLM..\n","\t\t\t\t\t NOPQRS.\n","\t\t\t\t\t TUV.WXY\n","\t\t\t\t\t \"\"\",\n","\t\t\t\t\t gridspec_kw = {'hspace' : 0.3}\n","\t\t\t\t\t )\n","\n","counter = 0\n","for i in stress_list:\n","\tfiles = [x for x in deseqouts if x.startswith(i+'control')]\n","\tfiles.sort()\n","\tfileD2 = pd.read_csv(wdir + files[0],\n","\t\t\t\t\t  sep = \"\\t\", header = 0, index_col = 0)\n","\tfileH2 = pd.read_csv(wdir + files[1],\n","\t\t\t\t\t  sep = \"\\t\", header = 0, index_col = 0)\n","\tD2index, H2index = set(fileD2.index.tolist()), set(fileH2.index.tolist())\n","\tstatus = []\n","\t\n","\t# Create sets\n","\tD2only = D2index - H2index\n","\tH2only = H2index - D2index\n","\tD2H2 = D2index & H2index\n","\t\n","\t# Subsets of df for D2only, H2only and D2H2\n","\tD2onlydf = fileD2[fileD2.index.isin(list(D2only))]\n","\tH2onlydf = fileH2[fileH2.index.isin(list(H2only))]\n","\tD2H2df = fileD2.append(fileH2)\n","\tD2H2df = D2H2df[D2H2df.index.isin(list(D2H2))]\n","\tD2H2df.sort_index(inplace=True)\n","\t\n","\t# Create list of lists of all differentially expressed genes with corresponding status\n","\tstress = \"Mpo_\" + files[0].split(\"controlD2\")[0]\n","\tfor j in D2only:\n","\t\tstatus.append([j, stress, 'D2', str(fileD2.loc[j, \"log2FoldChange\"]), str(fileD2.loc[j, \"padj\"]), \"N/A\", \"N/A\"])\n","\tfor k in H2only:\n","\t\tstatus.append([k, stress, 'H2', \"N/A\", \"N/A\", str(fileH2.loc[k, \"log2FoldChange\"]), str(fileH2.loc[k, \"padj\"])])\n","\tfor m in D2H2:\n","\t\tstatus.append([m, stress, 'D2H2', str(fileD2.loc[m, \"log2FoldChange\"]), str(fileD2.loc[m, \"padj\"]), str(fileH2.loc[m, \"log2FoldChange\"]), str(fileH2.loc[m, \"padj\"])])\n","\tstatus.sort()\n","\t\n","\t# Output file with genes and status: D2, H2 or D2H2\n","\twith open(wdir + \"sets/results/\" + stress + \".txt\", \"w+\") as filo:\n","\t\tfilo.write(\"gene\\tstress\\tstatus\\tL2FC_D2\\tpadj_D2\\tL2FC_H2\\tpadj_H2\\n\")\n","\t\tfor n in status:\n","\t\t\tfilo.write(\"\\t\".join(n) + \"\\n\")\n","\t\t\t\n","\t# Venn diagram\n","\tsp_ax = axd[f_axes[stress_list.index(stress.split('_')[1])]]\n","\tvenn2(subsets=(len(D2only), len(H2only), len(D2H2)),\n","\t   set_labels = ('D2', 'H2'),\n","\t   ax = sp_ax)\n","\tsp_ax.set_title(stress.split('_')[1], size=14)\n","\tplt.savefig(dir_path + \"figures/\" + \"suppfig4.png\", dpi = 600)\n","\t\n","\t# increase counter\n","\tcounter += 1"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"SFx5kliR9PQE"},"source":["# Sort genes according to whether they are the same in both controls\n","# adapted from compile_sigGenes_phase1n2.py\n","deseqdir = mpo_path\n","setdir = deseqdir + 'sets/results/'\n","setdir_safe = setdir.replace(' ', '\\ ')\n","if not os.path.exists(setdir):\n","\t!mkdir -p $setdir_safe\n","\n","merfile = open(dir_path + 'mercator/MpoProt.results.txt', 'r')\n","ofile = open(deseqdir + 'resSig_compiled.txt', 'w+')\n","efile = open(deseqdir + 'resSig_failed.txt', 'w+')\n","setfiles = [x for x in os.listdir(setdir) if '.txt' in x]\n","\n","ofile.write(\"\\t\".join(['gene', 'stress', 'L2FC_D2', 'L2FC_H2', 'annotation']) + \"\\n\")\n","efile.write(\"\\t\".join(['gene', 'stress', 'L2FC_D2', 'L2FC_H2', 'annotation']) + \"\\n\")\n","\n","def get_anno(gene):\n","\tgene = gene.lower()\n","\treturn meranno[gene]\n","\n","def up_down(val):\n","\tif val < 0:\n","\t\tstat = \"DOWN\"\n","\telif val > 0:\n","\t\tstat = \"UP\"\n","\telif math.isnan(val):\n","\t\tstat = \"NaN\"\n","\treturn stat\n","\n","meranno = {}\n","\n","for line in merfile:\n","\tif len(line.rstrip().split(\"\\t\")) == 5:\n","\t\tbincode, name, identifier, desc, ptype = line.rstrip().replace(\"'\", \"\").split(\"\\t\")\n","\t\tif identifier not in meranno:\n","\t\t\tmeranno[identifier] = [[bincode, desc]]\n","\t\telse:\n","\t\t\tmeranno[identifier].append([bincode, desc])\n","\n","for i in setfiles:\n","\tcontent = pd.read_csv(setdir + i, sep = \"\\t\", header = 0)\n","\tsigGenes = content[content['status'] == \"D2H2\"]\n","\tsigGenes['annotation'] = sigGenes['gene'].apply(get_anno)\n","\tsigGenes['L2FC_D2'] = sigGenes['L2FC_D2'].apply(up_down)\n","\tsigGenes['L2FC_H2'] = sigGenes['L2FC_H2'].apply(up_down)\n","\tsigGenes = sigGenes.drop(columns = ['status', 'padj_D2', 'padj_H2'])\n","\tfor index, row in sigGenes.iterrows():\n","\t\tif row['L2FC_D2'] != row['L2FC_H2']:\n","\t\t\tefile.write(\"\\t\".join([str(z) for z in row]) + \"\\n\")\n","\t\telse:\n","\t\t\tofile.write(\"\\t\".join([str(z) for z in row]) + \"\\n\")\n","\n","ofile.close()\n","efile.close()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Supp. Figure 5: Upset plot (up- and down-regulated)"],"metadata":{"id":"kAj_ImMzXZwT"}},{"cell_type":"code","metadata":{"id":"0NoknTLri0z9"},"source":["# Fig S5 A & B (adapted from upset.py)\n","import upsetplot\n","from collections import defaultdict\n","\n","wdir = dir_path + 'prep_files/mpo/deseq/'\n","odir = wdir + 'upset/'\n","odir_safe = dir_path_safe + 'prep_files/mpo/deseq/' + 'upset/'\n","\n","if not os.path.exists(odir):\n","\t!mkdir $odir_safe\n","\n","data = pd.read_csv(wdir + 'resSig_compiled.txt', sep = '\\t')\n","\n","### FUNCTIONS ###\n","\n","def upset_matrix(set_dict, stress_types):\n","\tupset_data_sub = upsetplot.from_contents({k: v for k, v in set_dict.items() if k in stress_types})\n","\treturn upset_data_sub # , fig=None\n","\n","\n","def plot_selected(cond_dict, cond_list, filename, title, orient = \"horizontal\", cutoff=50):\n","\tdf_set = upset_matrix(cond_dict, cond_list)\n","\tdf_set = df_set.sort_index()\n","\t\n","\t# preparation to output all data\n","\tindex_names = list(df_set.index.names)\n","\tindex_list = df_set.index.to_list()\n","\tset_count = Counter(index_list)\n","\tcounter_list = [[k, v] for k, v in set_count.items()]\n","\tcounter_list.sort(key = lambda x: x[1], reverse=True)\n","\t\n","\t# writing output to file\n","\twith open(filename + \"_matrix.txt\", \"w+\") as ofile:\n","\t\tofile.write(\"\\t\".join(index_names + ['count', 'genes']) + \"\\n\")\n","\t\tfor i in counter_list:\n","\t\t\tglist = df_set.loc[i[0],:].id.to_list()\n","\t\t\tofile.write(\"\\t\".join([str(int(x)) for x in i[0]] + [str(i[1]), str(glist)]) + \"\\n\")\n","\t\n","\t# writing top 50 to file\n","\tset_cutoff = set_count.most_common(cutoff)\n","\tselection = [x[0] for x in set_cutoff]\n","\twith open(filename + \"_top50.txt\", \"w+\") as cfile:\n","\t\tcfile.write(\"\\t\".join([str(index_names), 'count', 'genes']) + \"\\n\")\n","\t\tfor i in range(len(set_cutoff)):\n","\t\t\tglist = df_set.loc[set_cutoff[i][0],:].id.to_list()\n","\t\t\tcfile.write(\"\\t\".join([str(set_cutoff[i][0]), str(set_cutoff[i][1]), str(glist)]) + \"\\n\")\n","\t# selection for plotting\n","\tsel_matrix = df_set.loc[selection[0], :]\t\n","\tfor i in range(1, len(selection)):\n","\t\tsel_matrix = sel_matrix + df_set.loc[selection[i],:]\n","\tupsetplot.plot(sel_matrix, orientation = orient, sort_by = 'cardinality')\n","\tplt.title(title, size=20)\n","\tif \"upreg\" in filename:\n","\t\tfigname = 'c'\n","\telse:\n","\t\tfigname = 'd'\n","\tplt.savefig(dir_path + 'figures/FigS5' + figname + '.png', dpi=600)\n","\n","### END ###\n","\n","# Reshape data to have for every category,\n","cond_dict_U = defaultdict(list) # genres_movies\n","cond_dict_D = defaultdict(list)\n","for index, row in data.iterrows():\n","\tif row['L2FC_D2'] == 'UP':\n","\t\tcond_dict_U[row['stress'].split(\"_\")[1]].append(row['gene'])\n","\telif row['L2FC_D2'] == 'DOWN':\n","\t\tcond_dict_D[row['stress'].split(\"_\")[1]].append(row['gene'])\n","\t\t\n","all_stress_list = [x.split('_')[1] for x in data.stress.unique()]\n","\n","# initialise dictionaries of up and downregulated genes for each condition\n","cond_dict_U_set = dict()\n","cond_dict_D_set = dict()\n","for k, v in cond_dict_U.items():\n","    cond_dict_U_set[k] = set(v)\n","for k, v in cond_dict_D.items():\n","    cond_dict_D_set[k] = set(v)\n","\n","# Plot horizontal (default)\n","plot_selected(cond_dict = cond_dict_D_set,\n","\t\t\t  cond_list = all_stress_list,\n","\t\t\t  filename = odir + \"all_downreg\",\n","\t\t\t  title = \"Downregulated genes\")\n","\n","plot_selected(cond_dict = cond_dict_U_set,\n","\t\t\t  cond_list = all_stress_list,\n","\t\t\t  filename = odir + \"all_upreg\",\n","\t\t\t  title = \"Upregulated genes\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"_ENdRwbI7iDw"},"source":["### Supp Figs 6 & 7: Intersection of DEGs across single and cross stresses (up- and down-regulated)"]},{"cell_type":"code","metadata":{"id":"2DTxA1RKpwvR"},"source":["# Supp. figs 6 & 7 (adapted from indivenn_hm.py)\n","from collections import defaultdict\n","from matplotlib_venn import venn3\n","\n","wdir = dir_path + 'prep_files/mpo/deseq/'\n","odir = dir_path + 'figures/'\n","data = pd.read_csv(wdir + 'resSig_compiled.txt', sep = '\\t')\n","\n","# Mercator bin conversion\n","dicto = literal_eval(open(dir_path + 'prep_files/merdict.txt', 'r').read())\n","\n","all_s = [x.split(\"_\")[1] for x in data.stress.unique()]\n","single = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 1]\n","cross = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 2]\n","\n","data.annotation = data.annotation.apply(literal_eval)\n","data[\"mername\"] = data.annotation.apply(lambda x: dicto[int(x[0][0].split('.')[0])])\n","\n","dict_A = defaultdict(list)\n","dict_U = defaultdict(list)\n","dict_D = defaultdict(list)\n","\n","def sum_to_dict(dicto, stress, reg):\n","\tif reg == \"ALL\":\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress)]\n","\telse:\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress) & (data.L2FC_D2 == reg)]\n","\tdicto[stress].append(set(subset.gene.to_list()))\n","\tdicto[stress].append(subset.mername.to_list())\n","\n","def dict_to_df(dicto):\n","\tdf = pd.DataFrame.from_dict(dicto, orient='index', columns=[\"gene\", \"mername\"])\n","\treturn df\n","\n","for s in all_s:\n","\t#sum_to_dict(dict_A, s, \"ALL\")\n","\tsum_to_dict(dict_U, s, \"UP\")\n","\tsum_to_dict(dict_D, s, \"DOWN\")\n","\n","#df_A = dict_to_df(dict_A)\n","df_U = dict_to_df(dict_U)\n","df_D = dict_to_df(dict_D)\n","\n","def plot_venn(df, s1, s2, c1, title, axis):\n","\tvenn3([df.loc[s1].gene, df.loc[s2].gene, df.loc[c1].gene],\n","\t   (s1, s2, c1),\n","\t   ax = axis)\n","\taxis.set_title(title, size=20)\n","\n","# create subplots\n","xlen = 4\n","ylen = math.ceil(len(all_s)/4)\n","figw = xlen * 4\n","figh = ylen * 3.5\n","\n","a_axes = string.ascii_uppercase[:len(all_s)]\n","def plot_subplot(df, title_ext):\n","\taxa = plt.figure(constrained_layout=True,\n","\t\t\t\t figsize=(figw, figh)).subplot_mosaic(\n","\t\t\t\t\t \"\"\"\n","\t\t\t\t\t ABCD\n","\t\t\t\t\t EFGH\n","\t\t\t\t\t IJKL\n","\t\t\t\t\t MNOP\n","\t\t\t\t\t QR..\n","\t\t\t\t\t \"\"\"\n","\t\t\t\t\t )\n","\n","\tfor c in range(len(cross)):\n","\t\tst = cross[c]\n","\t\tplot_venn(df, st[0], st[1], st, st + title_ext, axa[a_axes[c]])\n","\tplt.savefig(odir+'supp_fig6or7' + title_ext +'.png', dpi=600)\n","\n","\n","#df_col = [[df_A, ''], [df_U, \"_upregulated\"], [df_D, \"_downregulated\"]]\n","df_col = [[df_U, \"_upregulated\"], [df_D, \"_downregulated\"]]\n","\n","for x in df_col:\n","\tdf_type, ext = x[0], x[1]\n","\tplot_subplot(df_type, ext)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"exDaka827fIh"},"source":["### Figure 2: Summary of DEGs in Marchantia and inter stress comparisons\n","\n"]},{"cell_type":"code","metadata":{"id":"A2YBe6KzehDL"},"source":["# Fig 2A and B (adapted from DGE_count_sizecorr.py)\n","cross = pd.read_csv(dir_path + 'prep_files/mpo/deseq/resSig_compiled.txt', sep='\\t')\n","\n","cross.stress = [x.split('_')[1] for x in list(cross.stress)]\n","stress_l = list(cross.stress.unique())\n","stress_l.sort(key=lambda x: len(x))\n","c_dgecount = cross.groupby(['stress', 'L2FC_D2']).count().gene.to_frame(name='count')\n","\n","unstacked = c_dgecount.unstack().reindex(stress_l)\n","ax = unstacked.plot.bar(figsize=(7,3), stacked=True, ylabel='Number of DEGs', color=['navy', 'firebrick'])\n","\n","handles, labels = ax.get_legend_handles_labels()\n","ax.legend(handles=handles[::-1], labels=[x.split(', ')[1].split(')')[0] for x in labels][::-1])\n","plt.savefig(dir_path + 'figures/fig4a.png', dpi=600)\n","\n","wdir = dir_path + 'prep_files/'\n","infile = 'phase1n2_measurements_nooutliers.txt'\n","\n","measurements = pd.read_csv(wdir + infile, sep='\\t')\n","\n","def ss_grab(stress, condition):\n","\t\"\"\"\n","\tSlice the relevant condition for \n","\tParameters\n","\t----------\n","\tstress : string\n","\t\tStress of interest.\n","\tcondition : string\n","\t\tCondition of interest.\n","\n","\tReturns\n","\t-------\n","\tsssub : dataframe\n","\t\tdf of single stress.\n","\n","\t\"\"\"\n","\tsssub = measurements[(measurements.Stress == stress) & (measurements.Condition == condition)]\n","\treturn sssub\n","\n","srep_keys = [['Cold', '3'],\n","\t\t\t ['Heat', '33'],\n","\t\t\t ['Salt', '40'],\n","\t\t\t ['Mannitol', '100'],\n","\t\t\t ['Light', '435'],\n","\t\t\t ['Dark', '3'],\n","\t\t\t ['Nitrogen', '0']]\n","srepdf = measurements[(measurements.Stress == 'Cold') & (measurements.Condition == '3')]\n","\n","for s, c in srep_keys[1:]:\n","\tsrepdf = pd.concat([srepdf, ss_grab(s, c)])\n","\n","crepdf = measurements[measurements.Condition == 'mixed']\n","m_nocon = pd.concat([srepdf, crepdf])\n","m_nocon.Stress = [x[0] if len(x) > 2 else x for x in m_nocon.Stress]\n","noHL = m_nocon[m_nocon.Stress != 'HL']\n","\n","avg_meas = noHL.groupby('Stress').mean()[['Parea', 'Earea']]\n","avg_meas.reindex(stress_l)\n","\n","totaldeg = cross.groupby(['stress']).count()['L2FC_D2']\n","totaldeg.reindex(stress_l)\n","avg_meas['totaldeg'] = totaldeg\n","avg_meas = avg_meas[['totaldeg', 'Parea', 'Earea']]\n","\n","# size plots by df\n","ax = avg_meas.plot.scatter(x='totaldeg', y='Parea', color='orange', label='Area (Day 15)')\n","\n","for ind, dat in avg_meas.iterrows():\n","    ax.annotate(ind, (dat['totaldeg'], dat['Parea']),\n","\t\t\t\txytext=(-4,-12), textcoords='offset points')\n","\n","avg_meas.plot.scatter(x='totaldeg', y='Earea', color='navy', label='Area (Day 21)', ax=ax)\n","for ind, dat in avg_meas.iterrows():\n","    ax.annotate(ind, (dat['totaldeg'], dat['Earea']),\n","\t\t\t\txytext=(-4,-12), textcoords='offset points')\n","\t\n","# size plot with regression\n","from scipy import stats\n","\n","q_colnames = avg_meas.columns.to_list()\n","dicto = {'Parea' : 'Day 15', 'Earea' : 'Day 21'}\n","\n","def plot_reg(df, title):\n","\tlabels = []\n","\tcol = q_colnames[0]\n","\tfor col2 in q_colnames[1:]:\n","\t\tplt.scatter(col, col2, data=df)\n","\t\t#m, c = np.polyfit(df[col], df[col2], 1)\n","\t\tm, c, r_value, p_value, std_err = stats.linregress(df[col], df[col2])\n","\t\tfor ind, dat in avg_meas.iterrows():\n","\t\t\tplt.annotate(ind, (dat[col], dat[col2]),\n","\t\t\t  xytext=(-4,-12), textcoords='offset points')\n","\t\tplt.plot(df[col], m*df[col] + c)\n","\t\tlabels.append(dicto[col2] + ' (R\\u00b2: ' + str(round(r_value**2,2)) + ', p: ' + str(round(p_value, 2)) + ')')\n","\tplt.legend(labels)\n","\tplt.xlabel('DEG count')\n","\tplt.ylabel('Size (mm\\u00b2)')\n","\tplt.savefig(dir_path + 'figures/fig4b.png', dpi=600)\n","\t\n","plot_reg(avg_meas, 'Number of DEGs vs Size')"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"8EibR5DSslOq"},"source":["# Figure 2 C-F (adapted from plot_venn_sum.py)\n","\n","from collections import defaultdict, Counter\n","from matplotlib_venn import venn3, venn3_circles\n","import random\n","from scipy import stats\n","\n","wdir = dir_path + 'prep_files/mpo/deseq/'\n","odir = wdir + 'indivenn_hm/'\n","data = pd.read_csv(wdir + 'resSig_compiled.txt', sep = '\\t')\n","\n","# Mercator\n","### DICTIONARY OF MERCATOR BINS ###\n","mfile = dir_path + 'mercator/MpoProt.results.txt'\n","\n","meranno = defaultdict(list)\n","merbin = defaultdict(list)\n","map2anno = {}\n","\n","merfile = open(mfile, 'r')\n","merfile.readline()\n","for line in merfile:\n","\tlinecon = line.rstrip().replace(\"'\", \"\").split(\"\\t\")\n","\tif len(linecon) == 5:\n","\t\tbincode, name, identifier, desc, ptype = linecon\n","\t\tmeranno[identifier].append(dicto[int(bincode.split('.')[0])])\n","\t\tmerbin[identifier].append('.'.join(bincode.split('.')[:2]))\n","\tif len(linecon[0].split('.')) == 2:\n","\t\tmap2anno[linecon[0]] = linecon[1]\n","\t\t\n","\n","all_s = [x.split(\"_\")[1] for x in data.stress.unique()]\n","single = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 1]\n","cross = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 2]\n","\n","data.annotation = data.annotation.apply(literal_eval)\n","data[\"mername\"] = data.annotation.apply(lambda x: [dicto[int(y[0].split('.')[0])] for y in x]) # different from cell above, hence the repetitive code\n","\n","dict_A = defaultdict(list)\n","dict_U = defaultdict(list)\n","dict_D = defaultdict(list)\n","\n","def sum_to_dict(dicto, stress, reg):\n","\tif reg == \"ALL\":\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress)]\n","\telse:\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress) & (data.L2FC_D2 == reg)]\n","\tdicto[stress].append(set(subset.gene.to_list()))\n","\tdicto[stress].append([y for x in subset.mername.to_list() for y in x])\n","\tdicto[stress].append(['.'.join(y[0].split('.')[:2]) for x in subset.annotation.to_list() for y in x])\n","\n","def dict_to_df(dicto):\n","\tdf = pd.DataFrame.from_dict(dicto, orient='index', columns=[\"gene\", \"mername\", \"mapbin2\"])\n","\treturn df\n","\n","for s in all_s:\n","\tsum_to_dict(dict_A, s, \"ALL\")\n","\tsum_to_dict(dict_U, s, \"UP\")\n","\tsum_to_dict(dict_D, s, \"DOWN\")\n","\n","df_A = dict_to_df(dict_A)\n","df_U = dict_to_df(dict_U)\n","df_D = dict_to_df(dict_D)\n","\n","# =============================================================================\n","# \n","# # Summary of stress response\n","# \n","# =============================================================================\n","\n","# Q1 : ji_cal(a, b) [%]\n","def ji_cal(a, b):\n","\t# jaccard index calculation\n","\treturn len(a&b) / len(a|b)\n","# Q2: |(A − AB)/A − (B − AB)/B| [% difference]\n","def suppInX(a, b, ab):\n","\treturn len((a-ab))/len(a) - len((b-ab))/len(b)\n","\n","# Q3: (AB - A - B) / AB [%]\n","def novel(a, b, ab):\n","\treturn len(ab - a - b) / len(ab)\n","\n","def q_col(df, colnames):\n","\t\"\"\"\n","\tCollates the params for each cross stress and output in df\n","\t\n","\tParameters\n","\t----------\n","\tdf : dataframe\n","\t\tdataframe to use (all genes, upreg/downreg only).\n","\n","\tReturns\n","\t-------\n","\tq_df : dataframe\n","\t\tdatafram containing JI of all cross stress.\n","\n","\t\"\"\"\n","\tq_dict = {}\n","\tfor c in range(len(cross)):\n","\t\tst = cross[c]\n","\t\ta = df.loc[st[0]].gene\n","\t\tb = df.loc[st[1]].gene\n","\t\tab = df.loc[st].gene\n","\t\tq_dict[st] = [\n","\t\t\tji_cal(a,b),\n","\t\t\t#perXInAB(a,ab),\n","\t\t\t#perXInAB(b,ab),\n","\t\t\tsuppInX(a,b,ab),\n","\t\t\t#suppInAB(a,b,ab),\n","\t\t\tnovel(a,b,ab)\n","\t\t\t\t]\n","\tq_df = pd.DataFrame.from_dict(q_dict, orient=\"index\", columns = colnames)\n","\treturn q_df\n","\n","q_colnames = [\"similarity\", \"suppression\", \"novel interaction\"]\n","\n","q_A, q_U, q_D = [q_col(df_A, q_colnames), q_col(df_U, q_colnames), q_col(df_D, q_colnames)]\n","qdf_col = [[q_U, \"Upregulated DEGs\"], [q_D, \"Downregulated DEGs\"]]\n","\n","def plot_q_subplots(df, outerax):\n","\tax = [axe[x] for x in outerax]\n","\t#plt.suptitle(title, fontsize=14)\n","\tfor i, axis in enumerate(ax):\n","\t\tif q_colnames[i] == \"suppression\":\n","\t\t\tsns.heatmap(df[q_colnames[i]].to_frame().transpose(), cmap='coolwarm', ax=axis)\n","\t\telse:\n","\t\t\tsns.heatmap(df[q_colnames[i]].to_frame().transpose(), cmap='Blues', ax=axis)\n","\t\taxis.set_yticklabels([q_colnames[i]], rotation=0)\n","\t\tcbar = axis.collections[0].colorbar\n","\t\tminval = round(df[q_colnames[i]].min(),2)\n","\t\tmaxval = round(df[q_colnames[i]].max(),2)\n","\t\n","\t\twhile round(minval*100,2) % 5 != 0:\n","\t\t\tminval += 0.01\n","\t\twhile round(maxval*100,2) % 5 != 0:\n","\t\t\tmaxval -= 0.01\n","\t\tcbar.set_ticks([minval, maxval])\n","\n","def plot_reg(df, outerax):\n","\tlabels = []\n","\tstatscol = []\n","\tfor i, col in enumerate(q_colnames[:-1]):\n","\t\tfor j, col2 in enumerate(q_colnames[i+1:]):\n","\t\t\touterax.scatter(col, col2, data=df)\n","\t\t\tm, c, r_value, p_value, std_err = stats.linregress(df[col], df[col2])\n","\t\t\tstatscol.append([m, c, r_value, p_value, std_err])\n","\t\t\t# m, c = np.polyfit(df[col], df[col2], 1)\n","\t\t\touterax.plot(df[col], m*df[col] + c)\n","\t\t\tlabels.append(col[:3] + ' v ' + col2[:3] + ' ($\\mathregular{R^{2}}$: '+str(round(r_value**2,1))+', p: ' + str('{:.2f}'.format(round(p_value,2))+')'))\n","\touterax.legend(labels, fontsize=\"x-small\")\n","\treturn statscol\n","\t\n","def dum_venn(a_b, c_a, c_b, c_ab, ax, col, title, ac=20, bc=20, cc=20):\n","\t'''\n","\n","\tParameters\n","\t----------\n","\ta_b : int\n","\t\tSize of A&B.\n","\tc_a : int\n","\t\tSize of A&C-B.\n","\tc_b : int\n","\t\tSize of B&C-A.\n","\tc_ab : int\n","\t\tSzie of C&(A&B).\n","\tax : axis handle\n","\t\tAxis handle of subplot to plot into.\n","\tcol : list\n","\t\tList containing lists of patch id and corresponding colour.\n","\ttitle:\n","\t\t\n","\tac : int, optional\n","\t\tSize of set a. The default is 20.\n","\tbc : int, optional\n","\t\tSize of set b. The default is 20.\n","\tcc : int, optional\n","\t\tSize of set c. The default is 20.\n","\n","\tReturns\n","\t-------\n","\tNone.\n","\n","\t'''\n","\t\n","# =============================================================================\n","# \ta_b = 8 # A&B\n","# \tc_a_b = 4 # C&A-B/ C&B-A\n","# \tc_ab = 3   # C&(A&B)\n","# =============================================================================\n","\tdum = list(string.ascii_uppercase + string.ascii_lowercase)\n","\trandom.shuffle(dum)\n","\ta = set(dum[:20])\n","\tb = set(list(a)[:a_b] + [x for x in dum if x not in a][:bc-a_b])\n","\tc = set(list(a-b)[:c_a] +\n","\t\t list(b-a)[:c_b] + list(a&b)[:c_ab] +\n","\t\t [x for x in dum if x not in a and x not in b][:cc-c_a-c_b-c_ab])\n","\n","\tv = venn3([a, b, c],\n","\t   ('A', 'B', 'AB'),\n","\t   ax = ax) # ax = axis\n","\tvenn3_circles([a, b, c], linewidth=1, color='k', ax=ax)\n","\tfor i in col:\n","\t\tv.get_patch_by_id(i[0]).set_color(i[1])\n","\tfor idx, subset in enumerate(v.subset_labels):\n","\t\tv.subset_labels[idx].set_visible(False)\n","\tax.set_title(title, fontsize=16)\n","\n","# =============================================================================\n","# #\n","# # Initialising subplot\n","# #\n","# =============================================================================\n","#figsize=(figw, figh)\n","\n","top_mosaic = [[\"v1\", \"v2\", \"v3\"]]\n","eq_mosaic = [\t[\"e1\", \"e2\", \"e3\"]\t]\n","middle_mosaic = [\n","\t[\"u1\", \"d1\"],\n","\t[\"u2\", \"d2\"],\n","\t[\"u3\", \"d3\"]\n","]\n","bottom_mosaic = [[\"r1\", \"r2\"]]\n","\n","figw, figh = 11, 9\n","fig = plt.figure(figsize=(figw, figh))\n","axc = fig.subplot_mosaic(\n","\ttop_mosaic,\n","\tgridspec_kw={\n","\t\t\"bottom\": 0.75,\n","\t\t\"top\": 1,\n","\t\t#\"wspace\": 0.5,\n","\t\t#\"hspace\": 0.5,\n","\t\t}\n","\t)\n","axd = fig.subplot_mosaic(\n","\teq_mosaic,\n","\tgridspec_kw={\n","\t\t\"bottom\": 0.55,\n","\t\t\"top\": 0.8,\n","\t\t#\"wspace\": 0.5,\n","\t\t#\"hspace\": 0.5,\n","\t\t}\n","\t)\n","\n","axe = fig.subplot_mosaic(\n","\tmiddle_mosaic,\n","\tgridspec_kw={\n","\t\t\"bottom\": 0.38,\n","\t\t\"top\": 0.6,\n","\t\t#\"wspace\": 0.5,\n","\t\t\"hspace\": 0.2,\n","\t\t}\n","\t)\n","axf = fig.subplot_mosaic(\n","\tbottom_mosaic,\n","\tgridspec_kw={\n","\t\t\"bottom\": 0,\n","\t\t\"top\": 0.3,\n","\t\t#\"wspace\": 0.5,\n","\t\t#\"hspace\": 0.5,\n","\t\t}\n","\t)\n","\n","for axy in ['e1', 'e2', 'e3']:\n","\taxd[axy].axis('off')\n","for axy in [\"v1\", \"v2\", \"v3\"]:\n","\taxc[axy].set_anchor('N')\n","\n","axd['e1'].text(0.39, 0.45, r\"$\\frac{A \\cap B}{A \\cup B}$\", fontsize=20)\n","axd['e2'].text(0.04, 0.45, r\"$\\frac{A-B-AB}{A}-\\frac{B-A-AB}{B}$\", fontsize=20)\n","axd['e3'].text(0.29, 0.45, r\"$\\frac{AB-A-B}{AB}$\", fontsize=20)\n","\n","for seq, x in enumerate(qdf_col):\n","\tdf_type, title = x[0], x[1]\n","\tplot_q_subplots(df_type, [x[seq] for x in middle_mosaic])\n","for axy in ['u1', 'u2', 'd1', 'd2']:\n","\taxe[axy].set_xticklabels([])\n","\taxe[axy].xaxis.set_visible(False)\n","for axy in ['d1', 'd2', 'd3']:\n","\taxe[axy].set_yticklabels([])\n","\taxe[axy].yaxis.set_visible(False)\n","\n","v1col = [['100', 'white'], ['110', 'limegreen'],\n","\t\t ['101', 'white'], ['111', 'limegreen'],\n","\t\t ['010', 'white'], ['011', 'white'],\n","\t\t ['001', 'white']]\n","v2col = [['100', 'red'], ['110', 'white'],\n","\t\t ['101', 'white'], ['111', 'white'],\n","\t\t ['010', 'cornflowerblue'], ['011', 'white'],\n","\t\t ['001', 'white']]\n","v3col = [['100', 'white'], ['110', 'white'],\n","\t\t ['101', 'white'], ['111', 'white'],\n","\t\t ['010', 'white'], ['011', 'white'],\n","\t\t ['001', 'darkorchid']]\n","dum_venn(a_b=8,c_a=4, c_b = 4,c_ab=3,\n","\t\t ax=axc['v1'], col=v1col, title='Similarity')\n","dum_venn(a_b=8, c_a=10, c_b=6, c_ab=3,\n","\t\t ax=axc['v2'], col=v2col, title='Suppression')\n","dum_venn(a_b=8, c_a=4, c_b=4, c_ab=3,\n","\t\t ax=axc['v3'], col=v3col, title='Novel interaction')\n","\n","reg_stats = []\n","for i, x in enumerate(qdf_col):\n","\tdf_type, title = x[0], x[1]\n","\tdf_abs = df_type[:]\n","\tdf_abs.suppression = abs(df_abs.suppression)\n","\treg_stats.append(plot_reg(df_abs, axf[bottom_mosaic[0][i]]))\n","\n","plt.savefig(dir_path+'figures/Fig5A_E.png', dpi=600)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"zd8VylDwyLWn"},"source":["# Fig 2H refer to l2_en_jaccard_hm.png\n","# =============================================================================\n","# \n","#  Summary: What is the dominant effect of each stress?\n","# \n","# =============================================================================\n","def sum_df(df):\n","\tdicto_sum = {}\n","\tfor ss in single:\n","\t\tori = [x for x in cross if ss in x]\n","\t\trelcross = [y for x in cross if ss in x for y in x if ss not in y]\n","\t\tsim = [ji_cal(df.loc[x[0]].gene, df.loc[x[1]].gene) for x in ori]\n","\t\tnov = [novel(df.loc[x[0]].gene, df.loc[x[1]].gene, df.loc[x].gene) for x in ori]\n","\t\tsup = [suppInX(df.loc[ss].gene, df.loc[x].gene, df.loc[ori[i]].gene) for i, x in enumerate(relcross)]\n","\t\tfor i, x in enumerate(relcross):\n","\t\t\tdicto_sum[ss+x] = [ss, sim[i], sup[i], nov[i]]\n","\t\t\t\n","\tdf_sum= pd.DataFrame.from_dict(dicto_sum, orient='index', columns=['stress', 'similarity', 'suppression', 'novel'])\n","\treturn df_sum\n","\n","cond=['similarity', 'suppression', 'novel']\n","df_list = [sum_df(df_U), sum_df(df_D)]\n","\n","# figs, axs = plt.subplots(3,2,\n","# \t\t\t\t\t\t sharex=True,\n","# \t\t\t\t\t\t sharey='row',\n","# \t\t\t\t\t\t constrained_layout=True,\n","# \t\t\t\t\t\t figsize=(7,6))\n","\n","# for i, x in enumerate(df_list):\n","# \tfor j, y in enumerate(cond):\n","# \t\tsns.violinplot(x='stress', y= y, data=x, ax= axs[j][i])\n","# for k in range(3):\n","# \taxs[k][1].set_ylabel('')\n","# for l in range(2):\n","# \tfor m in range(2):\n","# \t\taxs[l][m].set_xlabel('')\n","# axs[0][0].set_title('Upregulated', fontsize=14)\n","# axs[0][1].set_title('Downregulated', fontsize=14)\n","# plt.savefig(odir + 'venn_sum.png', dpi=600)\n","\n","# =============================================================================\n","# \n","# Enrichment\n","# \n","# =============================================================================\n","\n","from statsmodels.stats.multitest import multipletests\n","import math\n","import numpy as np\n","\n","ori_count = Counter([y for x in list(meranno.values()) for y in x])\n","mapbins = list(dicto.values())\n","\n","def sig_df(df, sigcol, merdict, mapbins):\n","\t\"\"\"\n","\tCalculates and correct mapman bin enrichment p-value for all stresses\n","\tReturns dataframe\n","\n","\tParameters\n","\t----------\n","\tdf : dataframe\n","\t\tdf containing genes and corresponding mapman bins of DEGs.\n","\tsigcol : str\n","\t\tcolumn name to use for enrichment\n","\tmerdict : dict\n","\t\tcorresponding dictionary of mapman annotation/ 2nd level bins to use\n","\tmapbins : list\n","\t\tlist of mapman annotation/bins to use\n","\n","\tReturns\n","\t-------\n","\tdf_sig : dataframe\n","\t\tdf summarising enrichment (corrected p-value) for each mapman bin (row)\n","\t\tand each stress (column).\n","\n","\t\"\"\"\n","\tsig_sum = {}\n","\tfor s in all_s:\n","\t\ts_count = Counter(df.loc[s][sigcol])\n","\t\tvalid_bins = list(s_count.keys()) # bins found in stress\n","\t\t# initialise count dicitonary\n","\t\tsig_count = {}\n","\t\tfor key in valid_bins:\n","\t\t\tsig_count[key] = 1\n","\t\t# random simulations\n","\t\tfor i in range(1000):\n","\t\t\tshuffle = list(merdict.values())\n","\t\t\trandom.shuffle(shuffle)\n","\t\t\tsub = shuffle[:len(df.loc[s].gene)]\n","\t\t\tsub_count = Counter([y for x in sub for y in x])\n","\t\t\tfor mapman in valid_bins:\n","\t\t\t\tif sub_count[mapman] >= s_count[mapman]:\n","\t\t\t\t\tsig_count[mapman] += 1\n","\t\t# p-value calculation\n","\t\tpval_coll = []\n","\t\tfor mapman in valid_bins:\n","\t\t\tpval = sig_count[mapman]/1000\n","\t\t\t# correction for pval > 1\n","\t\t\tif pval <= 1:\n","\t\t\t\tpval_coll.append(pval)\n","\t\t\telse:\n","\t\t\t\tpval_coll.append(float(round(pval)))\n","\t\t\n","\t\t# BH correction for multiple testing\n","\t\ty = multipletests(pvals=pval_coll, alpha=0.05, method=\"fdr_bh\")[1]\n","\t\tall_bins_corr_pval = []\n","\t\tfor mapman in mapbins:\n","\t\t\tif mapman in valid_bins:\n","\t\t\t\tall_bins_corr_pval.append(y[valid_bins.index(mapman)])\n","\t\t\telse:\n","\t\t\t\tall_bins_corr_pval.append(None)\n","\t\tsig_sum[s] = all_bins_corr_pval\n","\n","\tdf_sig = pd.DataFrame.from_dict(sig_sum, orient='index', columns=mapbins)\n","\treturn df_sig\n","\t\n","def chunk(uval, dval):\n","\tif math.isnan(uval) and math.isnan(dval):\n","\t\t# not differentially regulated\n","\t\tcat = 0\n","\telif  uval >= 0.05 and (dval >= 0.05 or math.isnan(dval)):\n","\t\t# not enriched\n","\t\tcat = 0\n","\telif  dval >= 0.05 and (uval >= 0.05 or math.isnan(uval)):\n","\t\t# not enriched\n","\t\tcat = 0\n","\telif  uval < 0.05 and dval < 0.05:\n","\t\t# differentially up and downregulated in bin\n","\t\tcat = 2\n","\telif dval < 0.05:\n","\t\t# differentially downregualted\n","\t\tcat = 1\n","\telif uval < 0.05:\n","\t\t# differentially upregulated\n","\t\tcat = 3\n","\treturn cat\n","\n","from matplotlib.colors import ListedColormap\n","cmap = ListedColormap([\"lightgray\", \"royalblue\", \"violet\", \"firebrick\"])\n","catno = 4\n","cbarticks = [(x/(catno*2))*(catno-1) for x in range(1,catno*2,2)]\n","\n","# =============================================================================\n","# \n","# Enrichment (Part 2: 2nd level Mapman)\n","# \n","# =============================================================================\n","mapbins2 = list(set([y for x in list(merbin.values()) for y in x]))\n","mapbins2.sort(key=lambda x: (int(x.split('.')[0]), int(x.split('.')[1])))\n","\n","df_sig_U2 = sig_df(df_U, 'mapbin2', merbin, mapbins2)\n","df_sig_D2 = sig_df(df_D, 'mapbin2', merbin, mapbins2)\n","df_sig_U2 = df_sig_U2.fillna(value=np.nan)\n","df_sig_D2 = df_sig_D2.fillna(value=np.nan)\n","\n","cat_dict2 = {}\n","for mapman in list(df_sig_U2.columns):\n","\tcat_col = []\n","\tfor stress in list(df_sig_U2.index):\n","\t\tuval, dval = df_sig_U2.loc[stress, mapman], df_sig_D2.loc[stress, mapman]\n","\t\tcat_col.append(chunk(uval,dval))\n","\tcat_dict2[mapman] = cat_col\n","\n","df_combined_sig2 = pd.DataFrame.from_dict(cat_dict2, orient='index', columns=list(df_sig_U2.index))\n","df_combined_sig2 = df_combined_sig2.loc[df_combined_sig2.max(axis=1) > 0,:]\n","df_combined_sig2 = df_combined_sig2.loc[(df_combined_sig2 > 0).sum(axis=1) >2,:]\n","df_combined_sig2.reset_index(inplace=True)\n","df_combined_sig2['index'] = df_combined_sig2['index'].apply(lambda x: map2anno[x])\n","df_combined_sig2.set_index('index', inplace=True)\n","\n","# =============================================================================\n","# \n","# Plotting 2nd level mapman enrichment (df_combined_sig2)\n","# \n","# =============================================================================\n","from scipy.cluster.hierarchy import dendrogram, linkage\n","from scipy.spatial.distance import squareform\n","from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n","\n","def jdistprep(df, axis):\n","\t'''\n","\tConvert df to sets (for calculation of JD of X axis)\n","\n","\tParameters\n","\t----------\n","\tdf : dataframe\n","\t\tdataframe of categorical variables to be converted to sets.\n","\taxis : int \n","\t\taxis to do sets on, 0 by column (default), 1 by row\n","\n","\tReturns\n","\t-------\n","\tdicto : dict\n","\t\tdictionary containing list of column values.\n","\n","\t'''\n","\tif axis == 1:\n","\t\tdf = df.T\n","\tdxkeys = df.columns.to_list()\n","\tdykeys = df.index.to_list()\n","\tdicto = {}\n","\tfor col in dxkeys:\n","\t\tdicto[col] = [dykeys[i] + '_' + str(x) for i, x in enumerate(df[col].to_list())]\n","\treturn [dicto, dxkeys]\n","\n","def jdist(df, axis=0):\n","\t'''\n","\tConstruct jaccard distance square matrix\n","\n","\tParameters\n","\t----------\n","\tdf : df\n","\t\tdataframe to be used for jiprep/ jdist calculation.\n","\taxis : int\n","\t\taxis to do sets on, 0 by column (default), 1 by row\n","\n","\tReturns\n","\t-------\n","\tlinkage_matrix : list\n","\t\tcondensed jaccard distance matrix.\n","\tjlist : list\n","\t\tlist of list (jaccard distance square matrix)\n","\tdicto : dict\n","\t\tdictionary of list\n","\n","\t'''\n","\tdicto, dxkeys = jdistprep(df, axis)\n","\tjlist = []\n","\tfor key in dxkeys:\n","\t\tcol = []\n","\t\tfor key2 in dxkeys:\n","\t\t\tset1, set2 = dicto[key], dicto[key2]\n","\t\t\tset1x = set([x for x in set1 if x.split('_')[1] != '0'])\n","\t\t\tset2x = set([x for x in set2 if x.split('_')[1] != '0'])\n","\t\t\tcol.append(1 - ji_cal(set1x, set2x))\n","\t\tjlist.append(col)\n","\tdists = squareform(jlist)\n","\tlinkage_matrix = linkage(dists, \"single\")\n","\treturn linkage_matrix, jlist, dicto\n","\n","def plot_dendro(linkage_matrix, ax, orient):\n","\t'''\n","\tPlots dendrogram into subplot\n","\n","\tParameters\n","\t----------\n","\tmat : list of lists\n","\t\tContains the square matrix of jaccard distances.\n","\tax : axes\n","\t\taxis of subplot to plot to.\n","\torient : str\n","\t\torientation of dendrogram to be plotted.\n","\n","\tReturns\n","\t-------\n","\tNone.\n","\n","\t'''\n","\t\n","\tdendrogram(linkage_matrix, no_labels=True, ax=ax, orientation=orient, color_threshold=0, above_threshold_color='#000000')\n","\n","xmat, xlist, xdict = jdist(df_combined_sig2)\n","ymat, ylist , ydict = jdist(df_combined_sig2, axis=1)\n","\n","yden = dendrogram(ymat, labels=df_combined_sig2.index.to_list(), orientation='left') #, color_threshold=0, above_threshold_color='#000000'\n","plt.show()\n","xden = dendrogram(xmat, labels=df_combined_sig2.columns.to_list(), orientation='top') #, color_threshold=0, above_threshold_color='#000000'\n","plt.show()\n","\n","yorder = yden['ivl']\n","xorder = xden['ivl']\n","df_sig2_reordered = df_combined_sig2[xorder]\n","df_sig2_reordered = df_sig2_reordered.reindex(yorder[::-1])\n","\n","fig, ax = plt.subplots(2,2,\n","\t\t\t\t\t   figsize=(7.5,8.5), # (width, height)\n","\t\t\t\t\t   constrained_layout=True,\n","\t\t\t\t\t   gridspec_kw={'width_ratios': [1.5, 5],'height_ratios': [1, 5]}) # constrained_layout=True,\n","ax0, ax1, ax2, ax3 = ax.flatten()\n","for i in [ax0, ax1, ax2]:\n","\ti.axis('off')\n","\n","plot_dendro(xmat, ax1, 'top')\n","plot_dendro(ymat, ax2, 'left')\n","# heatmap, tick and tick labels\n","hplot = ax3.imshow(df_sig2_reordered, cmap=cmap)\n","ax3.yaxis.tick_right()\n","ax3.set_ylabel(\"\")\n","ax3.set_xticks(np.arange(0, len(df_sig2_reordered.columns), 1))\n","ax3.set_yticks(np.arange(0, len(df_sig2_reordered), 1))\n","\n","\n","xcolour = ['k'] + ['firebrick']*4 + ['gray']*6 + ['mediumseagreen']*4 + ['k']*2 + ['darkorange']*3 + ['k']*2 + ['royalblue']*3\n","ax3.set_xticklabels(df_sig2_reordered.columns.to_list(), rotation=90)\n","\n","for i, tick_label in enumerate(ax3.get_xticklabels()):\n","\ttick_text = tick_label.get_text()\n","\ttick_label.set_color(xcolour[i])\n","\t\n","anno_long = ['annotated', 'cellulose', 'biosynthesis', 'hemicellulose', 'pectin', 'channels', 'degradation']\n","ax3.set_yticklabels([(lambda x: x.split('.')[1].lower() if x.split('.')[1] not in anno_long else x.lower())(x) for x in df_sig2_reordered.index.to_list()])\n","# cbar plotting and control\n","axins = inset_axes(ax0,\n","\t\t\t\t\twidth=\"40%\",  # width = 50% of parent_bbox width\n","\t\t\t\t\theight=\"90%\",  # height : 5%\n","\t\t\t\t\tloc = 'center')\n","cbar = fig.colorbar(hplot, cax=axins, ticks = cbarticks)\n","cbar.ax.set_yticklabels(['N', 'D', 'UD', 'U'])\n","\n","plt.savefig(dir_path+'figures/Fig5F.png', dpi=600, bbox_inches='tight') # no N"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Supp. Fig. 10: Expression of GRN TFs across experiments (clustered)"],"metadata":{"id":"bs2FjU8DQVlf"}},{"cell_type":"code","source":["import pandas as pd\n","import matplotlib.pyplot as plt\n","import seaborn as sns\n","from collections import Counter\n","\n","nw_anno = pd.read_csv(dir_path + 'prep_files/condensed_deg_anno_nwonly.txt', header=0, sep=\"\\t\")\n","tf_enrich = pd.read_csv(dir_path + 'prep_files/TF_group_enrich_anno.txt', header=0, sep=\"\\t\")\n","bin_enrich = pd.read_csv(dir_path + 'prep_files/group_TF_enrich_anno.txt', header=0, sep=\"\\t\")\n","\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t   'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN', 'N',\n","\t\t\t\t   'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","\n","stresses = [\"D\", \"H\", \"C\", \"L\", \"M\", \"S\", \"N\"]\n","clustered_stresses = [y for x in stresses for y in all_stresses if x in y]\n","\n","# plot tf\n","tfdf = nw_anno[-nw_anno.TF_anno.isna()]\n","#tfdf.sort_values(by=[\"TF_anno\"], inplace=True)\n","tfdf[\"longname\"] = tfdf.gene + \" (\" + tfdf.TF_anno + \")\"\n","\n","# clustering\n","tfplot = tfdf[[\"longname\"] + all_stresses]\n","tfplot.set_index(\"longname\", inplace=True)\n","\n","from scipy.cluster.hierarchy import dendrogram, linkage\n","from scipy.spatial.distance import squareform\n","from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n","import numpy as np\n","\n","def ji_cal(a, b):\n","\t# jaccard index calculation\n","\treturn len(a&b) / len(a|b)\n","\n","def jdistprep(df, axis):\n","\t'''\n","\tConvert df to sets (for calculation of JD of X axis)\n","\n","\tParameters\n","\t----------\n","\tdf : dataframe\n","\t\tdataframe of categorical variables to be converted to sets.\n","\taxis : int \n","\t\taxis to do sets on, 0 by column (default), 1 by row\n","\n","\tReturns\n","\t-------\n","\tdicto : dict\n","\t\tdictionary containing list of column values.\n","\n","\t'''\n","\tif axis == 1:\n","\t\tdf = df.T\n","\tdxkeys = df.columns.to_list()\n","\tdykeys = df.index.to_list()\n","\tdicto = {}\n","\tfor col in dxkeys:\n","\t\tdicto[col] = [dykeys[i] + '_' + str(x) for i, x in enumerate(df[col].to_list())]\n","\treturn [dicto, dxkeys]\n","\n","def jdist(df, axis=0):\n","\t'''\n","\tConstruct jaccard distance square matrix\n","\n","\tParameters\n","\t----------\n","\tdf : df\n","\t\tdataframe to be used for jiprep/ jdist calculation.\n","\taxis : int\n","\t\taxis to do sets on, 0 by column (default), 1 by row\n","\n","\tReturns\n","\t-------\n","\tlinkage_matrix : list\n","\t\tcondensed jaccard distance matrix.\n","\tjlist : list\n","\t\tlist of list (jaccard distance square matrix)\n","\tdicto : dict\n","\t\tdictionary of list\n","\n","\t'''\n","\tdicto, dxkeys = jdistprep(df, axis)\n","\tjlist = []\n","\tfor key in dxkeys:\n","\t\tcol = []\n","\t\tfor key2 in dxkeys:\n","\t\t\tset1, set2 = dicto[key], dicto[key2]\n","\t\t\tset1x = set([x for x in set1 if x.split('_')[1] != 'nan'])\n","\t\t\tset2x = set([x for x in set2 if x.split('_')[1] != 'nan'])\n","\t\t\tcol.append(1 - ji_cal(set1x, set2x))\n","\t\tjlist.append(col)\n","\tdists = squareform(jlist)\n","\tlinkage_matrix = linkage(dists, \"single\")\n","\treturn linkage_matrix, jlist, dicto\n","\n","def plot_dendro(linkage_matrix, ax, orient):\n","\t'''\n","\tPlots dendrogram into subplot\n","\n","\tParameters\n","\t----------\n","\tmat : list of lists\n","\t\tContains the square matrix of jaccard distances.\n","\tax : axes\n","\t\taxis of subplot to plot to.\n","\torient : str\n","\t\torientation of dendrogram to be plotted.\n","\n","\tReturns\n","\t-------\n","\tNone.\n","\n","\t'''\n","\t\n","\tdendrogram(linkage_matrix, no_labels=True, ax=ax, orientation=orient, color_threshold=0, above_threshold_color='#000000')\n","\n","#%% TFs with stress specific regulation\n","spec_path = elnet_dir + 'TF_Scond_spec.txt'\n","spec_df = pd.read_csv(spec_path, header=0, sep=\"\\t\")\n","spec_df[\"longname\"] = spec_df.gene + \" (\" + spec_df.TF_anno + \")\"\n","spec_TF = spec_df[spec_df.specificity > 0.7].longname.unique().tolist()\n","tfplot_subset = tfplot.loc[spec_TF]\n","\n","ymat, ylist, ydict = jdist(tfplot_subset, axis=1)\n","yden = dendrogram(ymat, labels=tfplot_subset.index.to_list(), orientation='left')\n","plt.show()\n","\n","yorder = yden['ivl']\n","tfplot_reordered = tfplot_subset.reindex(yorder[::-1])\n","\n","# duplicate stress columns\n","tfplot_reordered.replace({\"DOWN\": 1, \"UP\": 2}, inplace=True)\n","tfplot_reordered.fillna(0, inplace=True)\n","tfplot_all = tfplot_reordered[clustered_stresses]\n","\n","\n","#%%\n","# plotting\n","fig, ax = plt.subplots(2,2,\n","\t\t\t\t\t   figsize=(10,13), # (width, height)\n","\t\t\t\t\t   constrained_layout=True,\n","\t\t\t\t\t   gridspec_kw={'width_ratios': [1, 4],'height_ratios': [1, 7]}) # constrained_layout=True,\n","ax0, ax1, ax2, ax3 = ax.flatten()\n","for i in [ax0, ax1, ax2]:\n","\ti.axis('off')\n","\n","plot_dendro(ymat, ax2, 'left')\n","# heatmap, tick and tick labels\n","\n","hplot = ax3.imshow(tfplot_all, cmap=cmap)\n","# https://stackoverflow.com/questions/10354397/python-matplotlib-y-axis-ticks-on-right-side-of-plot\n","ax3.yaxis.tick_right()\n","ax3.set_ylabel(\"\")\n","ax3.set_xticks(np.arange(0, len(tfplot_all.columns), 1))\n","ax3.set_yticks(np.arange(0, len(tfplot_all), 1))\n","\n","ax3.set_xticklabels(tfplot_all.columns.to_list(), rotation=90)\n","ax3.set_yticklabels(tfplot_all.index.to_list())\n","\n","# xlabel colour\n","xcolour = ['royalblue']*6 + [\"darkkhaki\"]*7 + [\"mediumorchid\"]*7 + [\"firebrick\"]*5 + [\"k\"]*6 + [\"mediumseagreen\"]*5 + [\"darkorange\"]*7\n","for i, tick_label in enumerate(ax3.get_xticklabels()):\n","\ttick_text = tick_label.get_text()\n","\ttick_label.set_color(xcolour[i])\n","# ylabel colour\n","TF_fam_count = Counter(tfdf.TF_anno.tolist())\n","TF_fam_big = [k for k, v in TF_fam_count.items() if v > 3]\n","dict_col = {k : palette[i] for i, k in enumerate(TF_fam_big) if TF_fam_count[k] > 3}\n","ycolour = [dict_col[x.split(\"(\")[1].split(\")\")[0]] if x.split(\"(\")[1].split(\")\")[0] in dict_col else palette[-1] for  x in tfplot_all.index.to_list()]\n","\n","for i, tick_label in enumerate(ax3.get_yticklabels()):\n","\ttick_text = tick_label.get_text()\n","\ttick_label.set_color(ycolour[i])\n","\n","# cbar plotting and control\n","# https://matplotlib.org/stable/gallery/axes_grid1/demo_colorbar_with_inset_locator.html\n","axins = inset_axes(ax0,\n","\t\t\t\t\twidth=\"40%\",  # width = 50% of parent_bbox width\n","\t\t\t\t\theight=\"90%\",  # height : 5%\n","\t\t\t\t\tloc = 'center')\n","cbar = fig.colorbar(hplot, cax=axins, ticks = cbarticks)\n","cbar.ax.set_yticklabels([\"N\", 'D', 'U'])\n","plt.savefig(dir_path + 'figures/' + 'FigS10.png', dpi=600)"],"metadata":{"id":"ak1hXz5GQeEc"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 3B: Expression of GRN TFs across experiments"],"metadata":{"id":"XqxINfaLBErg"}},{"cell_type":"code","source":["# -*- coding: utf-8 -*-\n","\"\"\"\n","Created on Wed Apr 27 12:15:31 2022\n","\n","@author: Qiao Wen\n","Essentially, calculate if TF is significantly enriched for a bin\n","\"\"\"\n","\n","import pandas as pd\n","from ast import literal_eval\n","from collections import Counter\n","import random\n","\n","anno_path = dir_path + 'prep_files/condensed_deg_anno.txt'\n","anno_df = pd.read_csv(anno_path, header=0, sep=\"\\t\")\n","anno_df.merbin = anno_df.merbin.apply(literal_eval)\n","\n","nw_path = elnet_dir + 'union_0.8_ignoreAGRIS_topTF.txt'\n","nw_df = pd.read_csv(nw_path, header=0, sep=\"\\t\")\n","nw_genes = list(set(nw_df.predicted.to_list() + nw_df[\"Gene.ID\"].to_list()))\n","\n","nw_anno = anno_df[anno_df.gene.isin(nw_genes)]\n","bin_counts = Counter([x for x in nw_anno[nw_anno.TF_anno.isna()].merbinname.to_list()])\n","nw_anno[\"group\"] = [\" and \".join(literal_eval(x)) if bin_counts[x] > 2 else \"Other bin combinations\" for x in nw_anno.merbinname.to_list()]\n","nw_anno.to_csv(dir_path + 'prep_files/condensed_deg_anno_nwonly.txt', index=False, sep=\"\\t\")\n","\n","TF_list = nw_anno[-nw_anno.TF_anno.isna()].gene.to_list()\n","gene_groups = nw_anno[nw_anno.TF_anno.isna()].group.unique().tolist()\n","all_groups = nw_anno[nw_anno.TF_anno.isna()].group.to_list()\n","\n","# enrichment TF-bin\n","tf_enrich = pd.DataFrame(columns=[\"TF\", \"group\", \"pval\"])\n","for tf in TF_list:\n","\tsubset = nw_df[nw_df[\"Gene.ID\"] == tf]\n","\tannot = pd.merge(subset, nw_anno[[\"gene\", \"group\"]], left_on = ['predicted'], right_on = ['gene'])\n","\tactual_group = Counter(annot.group.to_list())\n","\tdict_count = {k:0 for k in actual_group.keys()}\n","\tfor i in range(1000):\n","\t\trandom_group = Counter(random.sample(all_groups, len(annot)))\n","\t\tfor g in actual_group.keys():\n","\t\t\tif random_group[g] >= actual_group[g]:\n","\t\t\t\tdict_count[g] += 1\n","\tfor k, v in dict_count.items():\n","\t\ttf_enrich.loc[len(tf_enrich)] = [tf, k, v/1000]\n","\n","merged_df = pd.merge(nw_df, nw_anno[[\"gene\", \"group\"]], left_on = ['predicted'], right_on = ['gene'])\n","count_df = merged_df.groupby([\"Gene.ID\", \"group\", \"TF_stat\"]).predicted.count()\n","\n","def get_details(tf, group):\n","\tstats = Counter(merged_df[(merged_df[\"Gene.ID\"] == tf) & (merged_df[\"group\"] == group)].TF_stat.to_list())\n","\tint_size = sum(stats.values())\n","\tratio = (stats[\"Activator\"] - stats[\"Repressor\"])/int_size\n","\tif ratio > 0.6:\n","\t\tratio_stat = \"VP\"\n","\telif ratio > 0.2:\n","\t\tratio_stat = \"P\"\n","\telif ratio > -0.2:\n","\t\tratio_stat = \"A\"\n","\telif ratio > -0.6:\n","\t\tratio_stat = \"N\"\n","\telse:\n","\t\tratio_stat = \"VN\"\n","\treturn([int_size, ratio, ratio_stat, tf + \" (meta) \" + group])\n","\n","tf_enrich[\"details\"] = tf_enrich.apply(lambda x: get_details(x.TF, x.group), axis=1)\n","tf_enrich[\"int_size\"] = tf_enrich.apply(lambda x: x.details[0], axis=1)\n","tf_enrich[\"ratio\"] = tf_enrich.apply(lambda x: x.details[1], axis=1)\n","tf_enrich[\"ratio_stat\"] = tf_enrich.apply(lambda x: x.details[2], axis=1)\n","tf_enrich[\"edge_name\"] = tf_enrich.apply(lambda x: x.details[3], axis=1)\n","\n","tf_enrich = tf_enrich[['edge_name', 'TF', 'group', 'pval', 'int_size', 'ratio', 'ratio_stat']]\n","tf_enrich.to_csv(dir_path + 'prep_files/TF_group_enrich_anno.txt', index=False, sep=\"\\t\")\n","\n","# enrichment bin-TF family\n","bin_enrich = pd.DataFrame(columns=[\"TF_fam\", \"group\", \"pval\"])\n","merged_df[\"TF_fam\"] = [nw_anno[nw_anno.gene == x].TF_anno.to_list()[0] for x in merged_df[\"Gene.ID\"].to_list()]\n","noTFs_as_predicted = merged_df[-merged_df.predicted.isin(TF_list)]\n","all_fams = noTFs_as_predicted.TF_fam.to_list()\n","for group in gene_groups:\n","\tsubset = noTFs_as_predicted[noTFs_as_predicted.group == group]\n","\tactual_fam = Counter(Counter(subset.TF_fam.to_list()))\n","\tdict_fam = {k:0 for k in actual_fam.keys()}\n","\tfor i in range(1000):\n","\t\trandom_fam = Counter(random.sample(all_fams, len(subset)))\n","\t\tfor g in actual_fam.keys():\n","\t\t\tif random_fam[g] >= actual_fam[g]:\n","\t\t\t\tdict_fam[g] += 1\n","\tfor k, v in dict_fam.items():\n","\t\tbin_enrich.loc[len(bin_enrich)] = [k, group, v/1000]\n","\n","group_size_dict = {k:len(noTFs_as_predicted[noTFs_as_predicted[\"group\"] == k]) for k in gene_groups}\n","def get_fam_details(fam, group):\n","\tstats = Counter(noTFs_as_predicted[(noTFs_as_predicted[\"TF_fam\"] == fam) & (noTFs_as_predicted[\"group\"] == group)].TF_stat.to_list())\n","\tint_size = sum(stats.values())\n","\tgroup_size = group_size_dict[group]\n","\trel_int_size = int_size/group_size\n","\tratio = (stats[\"Activator\"] - stats[\"Repressor\"])/int_size\n","\tif ratio > 0.6:\n","\t\tratio_stat = \"VP\"\n","\telif ratio > 0.2:\n","\t\tratio_stat = \"P\"\n","\telif ratio > -0.2:\n","\t\tratio_stat = \"A\"\n","\telif ratio > -0.6:\n","\t\tratio_stat = \"N\"\n","\telse:\n","\t\tratio_stat = \"VN\"\n","\treturn([int_size, rel_int_size, ratio, ratio_stat])\n","\n","bin_enrich[\"details\"] = bin_enrich.apply(lambda x: get_fam_details(x.TF_fam, x.group), axis=1)\n","bin_enrich[\"int_size\"] = bin_enrich.apply(lambda x: x.details[0], axis=1)\n","bin_enrich[\"rel_int_size\"] = bin_enrich.apply(lambda x: x.details[1], axis=1)\n","bin_enrich[\"ratio\"] = bin_enrich.apply(lambda x: x.details[2], axis=1)\n","bin_enrich[\"ratio_stat\"] = bin_enrich.apply(lambda x: x.details[3], axis=1)\n","\n","bin_enrich = bin_enrich[['TF_fam', 'group', 'pval', 'int_size', 'rel_int_size', 'ratio', 'ratio_stat']]\n","bin_enrich.to_csv(dir_path + 'prep_files/group_TF_enrich_anno.txt', index=False, sep=\"\\t\")"],"metadata":{"id":"uBsgiZgeD50e"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Selected TF representatives from 3C\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import numpy as np\n","\n","nw_anno = pd.read_csv(dir_path + 'prep_files/condensed_deg_anno_nwonly.txt', header=0, sep=\"\\t\")\n","tf_enrich = pd.read_csv(dir_path + 'prep_files/TF_group_enrich_anno.txt', header=0, sep=\"\\t\")\n","bin_enrich = pd.read_csv(dir_path + 'prep_files/group_TF_enrich_anno.txt', header=0, sep=\"\\t\")\n","\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t   'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN', 'N',\n","\t\t\t\t   'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","stresses = [\"D\", \"H\", \"C\", \"L\", \"M\", \"S\", \"N\"]\n","clustered_stresses = [y for x in stresses for y in all_stresses if x in y]\n","\n","# plot tf\n","tfdf = nw_anno[-nw_anno.TF_anno.isna()]\n","tfdf[\"longname\"] = tfdf.gene + \" (\" + tfdf.TF_anno + \")\"\n","\n","tf_list_ordered = [\"Mp3g21490.1\", \"Mp4g21220.1\", \"Mp6g02620.1\", \"Mp8g18310.1\",\n","\t\t\t\t   \"Mp4g00180.1\", \"Mp8g01770.1\", \"Mp7g00860.1\", \"Mp6g04650.1\",\n","\t\t\t\t   \"Mp5g18910.1\", \"Mp5g12480.1\", \"Mp3g06860.1\", \"Mp1g13740.1\",\n","\t\t\t\t   \"Mp8g14220.1\", \"Mp1g02860.1\", \"Mp8g04130.1\", \"Mp2g20960.1\",\n","\t\t\t\t   \"Mp4g12490.1\", \"Mp1g27060.1\", \"Mp5g14820.1\", \"Mp6g08290.1\",\n","\t\t\t\t   \"Mp3g10500.1\", \"Mp6g08290.1\", \"Mp3g10500.1\", \"Mp1g25720.1\",\n","\t\t\t\t   \"Mp1g13010.1\", \"Mp1g19730.1\", \"Mp3g19670.1\", \"Mp5g10280.1\",\n","\t\t\t\t   \"Mp1g04550.1\", \"Mp4g22280.1\", \"Mp1g25140.1\", \"Mp2g00890.1\"]\n","\n","tfdf_short = tfdf[tfdf.gene.isin(tf_list_ordered)]\n","tfdf_short.set_index(\"gene\", inplace=True)\n","tfdf_short_filtered = tfdf_short.reindex(tf_list_ordered)\n","tfdf_short_filtered = tfdf_short_filtered.iloc[:, 3:-2]\n","\n","# duplicate stress columns\n","tfdf_short_filtered.replace({\"DOWN\": 1, \"UP\": 2}, inplace=True)\n","tfdf_short_filtered.fillna(0, inplace=True)\n","tfplot_all = tfdf_short_filtered[clustered_stresses]\n","\n","# Plot heatmap\n","fig, ax = plt.subplots(figsize=(7,6))\n","\n","# heatmap, tick and tick labels\n","from matplotlib.colors import ListedColormap\n","cmap = ListedColormap([\"lightgray\", \"royalblue\", \"firebrick\"])\n","catno = 3\n","cbarticks = [(x/(catno*2))*(catno-1) for x in range(1,catno*2,2)]\n","\n","hplot = ax.imshow(tfplot_all, cmap=cmap)\n","ax.yaxis.tick_right()\n","ax.set_ylabel(\"\")\n","ax.set_xticks(np.arange(0, len(tfplot_all.columns), 1))\n","ax.set_yticks(np.arange(0, len(tfplot_all), 1))\n","\n","ax.set_xticklabels(tfplot_all.columns.to_list(), rotation=90)\n","ax.set_yticklabels(tfplot_all.index.to_list())\n","\n","xcolour = [\"k\"]*6 +  [\"firebrick\"]*5 + ['royalblue']*6 + [\"mediumseagreen\"]*5 + [\"mediumorchid\"]*7 + [\"darkkhaki\"]*7 + [\"darkorange\"]*7\n","for i, tick_label in enumerate(ax.get_xticklabels()):\n","\ttick_text = tick_label.get_text()\n","\ttick_label.set_color(xcolour[i])\n","\t\n","plt.savefig(dir_path + 'figures/' + 'Fig3B.png', dpi=600, bbox_inches=\"tight\")"],"metadata":{"id":"rU0ykuMFBIXH"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 3C: Specific expression in GRN TFs"],"metadata":{"id":"IkCNIWdIBtSL"}},{"cell_type":"code","source":["import pandas as pd\n","\n","nw_anno = pd.read_csv(dir_path + 'prep_files/condensed_deg_anno_nwonly.txt', header=0, sep=\"\\t\")\n","\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t   'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN', 'N',\n","\t\t\t\t   'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","full_stresses = [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen deficiency\"]\n","clustered_stresses = [[y for y in all_stresses if x in y] for x in stresses]\n","\n","# get TFs\n","tfdf = nw_anno[-nw_anno.TF_anno.isna()]\n","\n","def get_cluster_df(idx, cluster):\n","\tsubset = tfdf[[\"gene\", \"TF_anno\"] + cluster]\n","\tsubset.dropna(subset=cluster, how=\"all\", inplace=True)\n","\tcluster_len = len(cluster)\n","\tsubset[\"specificity\"] = subset.apply(lambda x: x.iloc[2:].value_counts().max()/cluster_len, axis=1)\n","\tsubset[\"top_direction\"] = subset.apply(lambda x: x.iloc[2:-1].value_counts().idxmax(), axis=1)\n","\tsubset[\"cluster\"] = full_stresses[idx]\n","\t\n","\treturn(subset[[\"cluster\", \"gene\", \"TF_anno\", \"specificity\", \"top_direction\"]])\n","\n","df_list = [get_cluster_df(i, cluster) for i, cluster in enumerate(clustered_stresses)]\n","cluster_df = pd.concat(df_list)\n","\n","cluster_df.to_csv(elnet_dir + 'TF_Scond_spec.txt', index=False, sep=\"\\t\")\n","\n","import numpy as np\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","spec_path = elnet_dir + 'TF_Scond_spec.txt'\n","spec_df = pd.read_csv(spec_path, header=0, sep=\"\\t\")\n","spec_df[\"longname\"] = spec_df.gene + \" (\" + spec_df.TF_anno + \")\"\n","spec_df_subset = spec_df[spec_df.specificity > 0.7]\n","spec_TF = spec_df_subset.longname.unique().tolist()\n","\n","def update_df(idx, col, val, df):\n","\tdf.loc[idx, col] = val\n","\t\n","# initialise dataframe\n","plot_spec_df = pd.DataFrame(np.nan, index = spec_TF, columns = stresses)\n","plot_dir_df = pd.DataFrame(np.nan, index = spec_TF, columns = stresses)\n","plot_spec_dir_df = pd.DataFrame(np.nan, index = spec_TF, columns = stresses)\n","# update dataframe\n","spec_df_subset.apply(lambda x: update_df(df=plot_spec_df, idx=x.longname, col=x.cluster[0], val=x.specificity), axis=1)\n","spec_df_subset.apply(lambda x: update_df(df=plot_dir_df, idx=x.longname, col=x.cluster[0], val=x.top_direction), axis=1)\n","spec_df_subset.apply(lambda x: update_df(df=plot_spec_dir_df, idx=x.longname, col=x.cluster[0], val=x.top_direction), axis=1)\n","# change values to suit plotting\n","plot_dir_df.replace({\"DOWN\": 1, \"UP\": 2}, inplace=True)\n","plot_dir_df.fillna(0, inplace=True)\n","\n","plot_spec_dir_df.replace({\"DOWN\": -1, \"UP\": 1}, inplace=True)\n","plot_dendro = plot_spec_dir_df * plot_spec_df\n","plot_dendro.fillna(0, inplace=True)\n","#%% colour palette\n","from collections import Counter\n","palette = sns.color_palette()\n","palette.append((166/255,163/255,162/255))\n","TF_list =  plot_spec_df.index.to_list()\n","TF_fam_count = Counter([x.split(\"(\")[1].split(\")\")[0] for x in TF_list])\n","TF_fam_big = [k for k, v in TF_fam_count.items() if v > 3]\n","dict_col = {k : palette[i] for i, k in enumerate(TF_fam_big) if TF_fam_count[k] > 3}\n","\n","# to get status of TF (ycolour)\n","elnet_path = elnet_dir + union_0.8_ignoreAGRIS_topTF.txt'\n","elnet_df = pd.read_csv(elnet_path, header=0, sep=\"\\t\")\n","\n","def get_TFstat(tf):\n","\tstat_dict = Counter(elnet_df[elnet_df[\"Gene.ID\"] == tf].TF_stat)\n","\tTF_stat_max = max(stat_dict, key=stat_dict.get)\n","\tif TF_stat_max == 'Activator':\n","\t\treturn(\"limegreen\")\n","\telif TF_stat_max == 'Ambiguous':\n","\t\treturn(\"yellow\")\n","\telse:\n","\t\treturn(\"firebrick\")\n","\n","# Count number of Activators, Ambiguous and Repressors in all 95 Mpo TFs.\n","tf_list = elnet_df[\"Gene.ID\"].unique().tolist()\n","\n","def get_TFstat_str(tf):\n","\tstat_dict = Counter(elnet_df[elnet_df[\"Gene.ID\"] == tf].TF_stat)\n","\tTF_stat_max = max(stat_dict, key=stat_dict.get)\n","\treturn TF_stat_max\n","\n","tf_stat_str_list = [get_TFstat_str(x) for x in tf_list]\n","stat_count = Counter(tf_stat_str_list) # Counter({'Activator': 75, 'Repressor': 19, 'Ambiguous': 1})\n","\t\n","ycolour = [get_TFstat(tf) for tf in [x.split(\" (\")[0] for x in TF_list]]\n","\n","#%% Plotting\n","g = sns.clustermap(plot_dendro, cmap=\"coolwarm\",\n","\t\t\t\t   yticklabels=True, row_colors = ycolour,\n","\t\t\t\t   figsize=(8, 14), cbar_pos=(0.02, 0.9, 0.05, 0.08),\n","\t\t\t\t   dendrogram_ratio=(0.2,0.1)) \n","\n","ax = g.ax_heatmap\n","\n","for tick_label in ax.get_yticklabels():\n","\ttick_text = tick_label.get_text().split(\"(\")[1].split(\")\")[0]\n","\tif tick_text in dict_col:\n","\t\ttick_label.set_color(dict_col[tick_text])\n","\telse:\n","\t\ttick_label.set_color(palette[-1])\n","plt.savefig(dir_path + 'figures/' + 'Fig3C.png', dpi=600)\n","\n","#%% Fill df with more info\n","plot_dendro[\"spec_count\"] = plot_dendro.apply(lambda x: sum(x != 0), axis=1)\n","plot_dendro.to_csv(elnet_dir + 'TF_Scond_specificity_07.txt', sep=\"\\t\")"],"metadata":{"id":"twO61lVoBv_i"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Supp. Fig 12: Influence of TFs"],"metadata":{"id":"RW2j4eR4QqUf"}},{"cell_type":"code","source":["\"\"\"\n","Extent of influence for each transcription factor\n","\"\"\"\n","\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import numpy as np\n","\n","tf_nw_p = dir_path + 'prep_files/network_cf22_spec.txt'\n","top_nw_p = elnet_dir + 'union_0.8_ignoreAGRIS_topTF.txt'\n","\n","tf_nw = pd.read_csv(tf_nw_p, header=0, sep=\"\\t\")\n","top_nw = pd.read_csv(top_nw_p, header=0, sep=\"\\t\") \n","\n","tf_list = top_nw[\"Gene.ID\"].unique().tolist()\n","\n","# get number of direct genes\n","tf_direct_genes = [top_nw[top_nw[\"Gene.ID\"] == tf].predicted.to_list() for tf in tf_list]\n","tf_direct_gcount = [len(x) for x in tf_direct_genes]\n","tf_direct_tf = [[gene for gene in top_nw[top_nw[\"Gene.ID\"] == tf].predicted.to_list() if gene in tf_list] for tf in tf_list]\n","tf_direct_tfcount = [len(x) for x in tf_direct_tf]\n","direct_genedf = pd.DataFrame({\"gene_list\" : tf_direct_genes,\n","\t\t\t\t\t\t \"gene_count\" : tf_direct_gcount,\n","\t\t\t\t\t\t \"tf_list\" : tf_direct_tf,\n","\t\t\t\t\t\t \"tf_count\" : tf_direct_tfcount},\n","\t\t\t\t\t\t index=tf_list)\n","# get number of direct TFs\n","tf_direct_tf = [[gene for gene in tf_nw[tf_nw[\"Gene.ID\"] == tf].predicted.to_list() if gene in tf_list] for tf in tf_list]\n","tf_direct_tfcount = [len(x) for x in tf_direct_tf]\n","tf_incoming_tf = [[gene for gene in tf_nw[tf_nw.predicted == tf][\"Gene.ID\"].to_list() if gene in tf_list] for tf in tf_list]\n","tf_incoming_tfcount = [len(x) for x in tf_incoming_tf]\n","direct_tfdf = pd.DataFrame({\"tf_list\" : tf_direct_tf,\n","\t\t\t\t\t\t \"tf_count\" : tf_direct_tfcount,\n","\t\t\t\t\t\t \"incoming_tf_list\": tf_incoming_tf,\n","\t\t\t\t\t\t \"incoming_tf_count\": tf_incoming_tfcount},\n","\t\t\t\t\t\t index=tf_list)\n","\n","# plot number of outgoing edges (TFs)\n","tf_count = direct_tfdf.tf_count.value_counts()\n","tf_count_index = list(tf_count.keys())\n","tf_count_index.sort()\n","tf_count_values = [tf_count[x] for x in tf_count_index]\n","plt.bar(tf_count_index, tf_count_values)\n","plt.title(\"Number of TFs in first neighbourhood\")\n","plt.xlabel(\"Count\")\n","plt.ylabel(\"Frequency\")\n","#plt.savefig('G:/My Drive/Projects/Marchantia_2019/grn/TF_hierarchy/nbh_tf.png', dpi=600)\n","plt.show()\n","\n","# get number of outgoing edges (genes and TFs)\n","\n","all_direct_genes = [list(set(direct_genedf.loc[tf].gene_list + direct_tfdf.loc[tf].tf_list + [z for y in [direct_genedf.loc[x].gene_list for x in direct_tfdf.loc[tf].tf_list] for z in y])) for tf in tf_list]\n","len_all_direct_genes = [len(x) for x in all_direct_genes]\n","\t\n","direct_tfdf[\"direct_genes\"] = all_direct_genes\n","direct_tfdf[\"gene_count\"] = len_all_direct_genes\n","\n","direct_tfdf.to_csv(dir_path + 'prep_files/influence_nw_anno.txt', sep=\"\\t\")\n","\n","#%% get stress specificity source (for labelling)\n","desc_path = elnet_dir + 'TF_Scond_specificity_07.txt'\n","desc_file = pd.read_csv(desc_path, header=0, sep=\"\\t\")\n","desc_file.columns = [\"desc\"] + desc_file.columns.to_list()[1:]\n","desc_file[\"gene\"] = desc_file.desc.apply(lambda x: x.split(\" \")[0])\n","\n","stress_index = desc_file.iloc[:,1:8].columns.to_list()\n","def stress_source(df_row):\n","\tstatus = [stress_index[i] for i, x in enumerate(list(df_row)) if x != 0]\n","\tdesc = \", \".join(status)\n","\treturn desc\n","desc_file[\"stress_source\"] = desc_file.iloc[:,1:8].apply(lambda x: stress_source(x), axis=1)\n","\n","#%% get stress score: estimated impact (and similarity?)\n","valid_genes = desc_file.gene.to_list()\n","direct_tfdf[\"stress_source\"] = [desc_file[desc_file.gene == tf].stress_source.values[0] if tf in valid_genes else \"X\" for tf in direct_tfdf.index.to_list()]\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","stress_tf_list = [direct_tfdf[[s in x for x in direct_tfdf.stress_source.to_list()]] for s in stresses]\n","stress_tf_count = [len(x) for x in stress_tf_list]\n","stress_tfnbh_list = [list(set([z for y in direct_tfdf[[s in x for x in direct_tfdf.stress_source.to_list()]].tf_list.to_list() for z in y + direct_tfdf[[s in x for x in direct_tfdf.stress_source.to_list()]].index.to_list()])) for s in stresses]\n","stress_tfnbh_count = [len(x) for x in stress_tfnbh_list]\n","# [40, 9, 16, 19, 34, 59, 56]\n","stress_nbh_list = [list(set([y for x in direct_genedf.loc[stress_tfnbh_list[i]].gene_list.to_list() + stress_tfnbh_list[i] for y in x])) for i, s in enumerate(stresses)]\n","stress_nbh_count = [len(x) for x in stress_nbh_list]\n","\n","stress_stats = pd.DataFrame({\"stress\" : stresses,\n","\t\t\t\t\t\t\t \"tf\" : stress_tf_list,\n","\t\t\t\t\t\t\t \"tf_count\": stress_tf_count,\n","\t\t\t\t\t\t\t \"tf_nbh\": stress_tfnbh_list,\n","\t\t\t\t\t\t\t \"tf_nbh_count\": stress_tfnbh_count,\n","\t\t\t\t\t\t\t \"gene_nbh\": stress_nbh_list,\n","\t\t\t\t\t\t\t \"gene_nbh_count\": stress_nbh_count\n","\t\t\t\t\t\t\t })\n","\n","stress_stats.plot.bar(\"stress\", logy=True)\n","plt.legend([\"TFs\", \"TFs in neighbourhood\", \"Genes in neighbourhood\"], bbox_to_anchor=(1.0, 1.02))\n","plt.ylabel(\"log(Count)\")\n","plt.xticks([i for i in range(len(stresses))], [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen\\ndeficiency\"], rotation=\"horizontal\")\n","plt.savefig(dir_path + 'figures/' + 'FigS12.png', dpi=600, bbox_inches=\"tight\")"],"metadata":{"id":"hevrKQcXQxfD"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Supp. Fig 13: Robustly responding second-level MapMan bins across the 7 abiotic stresses"],"metadata":{"id":"bhj_ByC1QrnY"}},{"cell_type":"code","source":["from ast import literal_eval\n","import random\n","from collections import defaultdict, Counter\n","from statsmodels.stats.multitest import multipletests\n","import math\n","import numpy as np\n","import pandas as pd\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","\n","dicto = literal_eval(open(dir_path + 'mercator/merdict.txt', \"r\").read())\n","\n","meranno = defaultdict(list)\n","merbin = defaultdict(list)\n","map2anno = {}\n","\n","mfile = dir_path + 'mercator/MpoProt.results.txt'\n","merfile = open(mfile, 'r')\n","merfile.readline()\n","for line in merfile:\n","\tlinecon = line.rstrip().replace(\"'\", \"\").split(\"\\t\")\n","\tif len(linecon) == 5:\n","\t\tbincode, name, identifier, desc, ptype = linecon\n","\t\tmeranno[identifier].append(dicto[int(bincode.split('.')[0])])\n","\t\tmerbin[identifier].append('.'.join(bincode.split('.')[:2]))\n","\tif len(linecon[0].split('.')) == 2:\n","\t\tmap2anno[linecon[0]] = linecon[1]\n","\n","#%% Initialise\n","wdir = mpo_path\n","data = pd.read_csv(wdir + 'resSig_compiled.txt', sep = '\\t')\n","\n","all_s = [x.split(\"_\")[1] for x in data.stress.unique()]\n","single = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 1]\n","cross = [x.split(\"_\")[1] for x in data.stress.unique() if len(x.split(\"_\")[1]) == 2]\n","\n","data.annotation = data.annotation.apply(literal_eval)\n","data[\"mername\"] = data.annotation.apply(lambda x: [dicto[int(y[0].split('.')[0])] for y in x])\n","\n","dict_A = defaultdict(list)\n","dict_U = defaultdict(list)\n","dict_D = defaultdict(list)\n","\n","def sum_to_dict(dicto, stress, reg):\n","\tif reg == \"ALL\":\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress)]\n","\telse:\n","\t\tsubset = data[(data.stress == \"Mpo_\" + stress) & (data.L2FC_D2 == reg)]\n","\tdicto[stress].append(set(subset.gene.to_list()))\n","\tdicto[stress].append([y for x in subset.mername.to_list() for y in x])\n","\tdicto[stress].append(['.'.join(y[0].split('.')[:2]) for x in subset.annotation.to_list() for y in x])\n","\n","def dict_to_df(dicto):\n","\tdf = pd.DataFrame.from_dict(dicto, orient='index', columns=[\"gene\", \"mername\", \"mapbin2\"])\n","\treturn df\n","\n","for s in all_s:\n","\tsum_to_dict(dict_A, s, \"ALL\")\n","\tsum_to_dict(dict_U, s, \"UP\")\n","\tsum_to_dict(dict_D, s, \"DOWN\")\n","\n","df_U = dict_to_df(dict_U)\n","df_D = dict_to_df(dict_D)\n","\n","#%% Definitions for enrichment\n","def sig_df(df, sigcol, merdict, mapbins):\n","\t\"\"\"\n","\tCalculates and correct mapman bin enrichment p-value for all stresses\n","\tReturns dataframe\n","\n","\tParameters\n","\t----------\n","\tdf : dataframe\n","\t\tdf containing genes and corresponding mapman bins of DEGs.\n","\tsigcol : str\n","\t\tcolumn name to use for enrichment\n","\tmerdict : dict\n","\t\tcorresponding dictionary of mapman annotation/ 2nd level bins to use\n","\tmapbins : list\n","\t\tlist of mapman annotation/bins to use\n","\n","\tReturns\n","\t-------\n","\tdf_sig : dataframe\n","\t\tdf summarising enrichment (corrected p-value) for each mapman bin (row)\n","\t\tand each stress (column).\n","\n","\t\"\"\"\n","\tsig_sum = {}\n","\tfor s in all_s:\n","\t\ts_count = Counter(df.loc[s][sigcol])\n","\t\tvalid_bins = list(s_count.keys()) # bins found in stress\n","\t\t# initialise count dicitonary\n","\t\tsig_count = {}\n","\t\tfor key in valid_bins:\n","\t\t\tsig_count[key] = 1\n","\t\t# random simulations\n","\t\tfor i in range(1000):\n","\t\t\tshuffle = list(merdict.values())\n","\t\t\trandom.shuffle(shuffle)\n","\t\t\tsub = shuffle[:len(df.loc[s].gene)]\n","\t\t\tsub_count = Counter([y for x in sub for y in x])\n","\t\t\tfor mapman in valid_bins:\n","\t\t\t\tif sub_count[mapman] >= s_count[mapman]:\n","\t\t\t\t\tsig_count[mapman] += 1\n","\t\t# p-value calculation\n","\t\tpval_coll = []\n","\t\tfor mapman in valid_bins:\n","\t\t\tpval = sig_count[mapman]/1000\n","\t\t\t# correction for pval > 1\n","\t\t\tif pval <= 1:\n","\t\t\t\tpval_coll.append(pval)\n","\t\t\telse:\n","\t\t\t\tpval_coll.append(float(round(pval)))\n","\t\t\n","\t\t# BH correction for multiple testing\n","\t\ty = multipletests(pvals=pval_coll, alpha=0.05, method=\"fdr_bh\")[1]\n","\t\tall_bins_corr_pval = []\n","\t\tfor mapman in mapbins:\n","\t\t\tif mapman in valid_bins:\n","\t\t\t\tall_bins_corr_pval.append(y[valid_bins.index(mapman)])\n","\t\t\telse:\n","\t\t\t\tall_bins_corr_pval.append(None)\n","\t\tsig_sum[s] = all_bins_corr_pval\n","\n","\tdf_sig = pd.DataFrame.from_dict(sig_sum, orient='index', columns=mapbins)\n","\treturn df_sig\n","\t\n","def chunk(uval, dval):\n","\tif math.isnan(uval) and math.isnan(dval):\n","\t\t# not differentially regulated\n","\t\tcat = 0\n","\telif  uval >= 0.05 and (dval >= 0.05 or math.isnan(dval)):\n","\t\t# not enriched\n","\t\tcat = 0\n","\telif  dval >= 0.05 and (uval >= 0.05 or math.isnan(uval)):\n","\t\t# not enriched\n","\t\tcat = 0\n","\telif  uval < 0.05 and dval < 0.05:\n","\t\t# differentially up and downregulated in bin\n","\t\tcat = 2\n","\telif dval < 0.05:\n","\t\t# differentially downregualted\n","\t\tcat = 1\n","\telif uval < 0.05:\n","\t\t# differentially upregulated\n","\t\tcat = 3\n","\treturn cat\n","\n","\n","# =============================================================================\n","# \n","# Enrichment (Part 2: 2nd level Mapman)\n","# \n","# =============================================================================\n","mapbins2 = list(set([y for x in list(merbin.values()) for y in x]))\n","mapbins2.sort(key=lambda x: (int(x.split('.')[0]), int(x.split('.')[1])))\n","\n","df_sig_U2 = sig_df(df_U, 'mapbin2', merbin, mapbins2)\n","df_sig_D2 = sig_df(df_D, 'mapbin2', merbin, mapbins2)\n","df_sig_U2 = df_sig_U2.fillna(value=np.nan)\n","df_sig_D2 = df_sig_D2.fillna(value=np.nan)\n","\n","cat_dict2 = {}\n","for mapman in list(df_sig_U2.columns):\n","\tcat_col = []\n","\tfor stress in list(df_sig_U2.index):\n","\t\tuval, dval = df_sig_U2.loc[stress, mapman], df_sig_D2.loc[stress, mapman]\n","\t\tcat_col.append(chunk(uval,dval))\n","\tcat_dict2[mapman] = cat_col\n","\n","df_combined_sig2 = pd.DataFrame.from_dict(cat_dict2, orient='index', columns=list(df_sig_U2.index))\n","df_combined_sig2 = df_combined_sig2.loc[df_combined_sig2.max(axis=1) > 0,:]\n","df_combined_sig2 = df_combined_sig2.loc[(df_combined_sig2 > 0).sum(axis=1) >2,:]\n","\n","df_combined_sig2.reset_index(inplace=True)\n","df_combined_sig2['index'] = df_combined_sig2['index'].apply(lambda x: map2anno[x])\n","df_combined_sig2.set_index('index', inplace=True)\n","df_combined_sig2.replace({0: np.nan, 1: \"DOWN\", 2: \"UP_DOWN\", 3: \"UP\"}, inplace=True)\n","\n","df_combined_sig2.to_csv(mpo_path + 'map2_collapsed_df.txt', sep=\"\\t\")\n","\n","#%% collapse to cluster\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t   'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN', 'N',\n","\t\t\t\t   'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","full_stresses = [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen deficiency\"]\n","clustered_stresses = [[y for y in all_stresses if x in y] for x in stresses]\n","\n","def get_cluster_df(idx, cluster, df):\n","\tsubset = df[cluster]\n","\tsubset.dropna(subset=cluster, how=\"all\", inplace=True)\n","\tcluster_len = len(cluster)\n","\t# https://www.codegrepper.com/code-examples/python/find+max+value+index+in+value+count+pandas\n","\tsubset[\"specificity\"] = subset.apply(lambda x: x.value_counts().max()/cluster_len, axis=1)\n","\tsubset[\"top_direction\"] = subset.apply(lambda x: x.iloc[:-1].value_counts().idxmax(), axis=1)\n","\tsubset[\"cluster\"] = full_stresses[idx]\n","\tsubset.reset_index(inplace=True)\n","\treturn(subset[[\"index\", \"cluster\", \"specificity\", \"top_direction\"]])\n","\n","df_list = [get_cluster_df(i, cluster, df_combined_sig2) for i, cluster in enumerate(clustered_stresses)]\n","cluster_df = pd.concat(df_list)\n","cluster_df.to_csv(mpo_path + \"mapman_spec_df.txt\", sep=\"\\t\", index=None)\n","#%% preparation for plotting\n","def update_df(idx, col, val, df):\n","\tdf.loc[idx, col] = val\n","\n","mapbin2_list = [x.split(\".\")[1] for x in cluster_df[\"index\"].unique().tolist()]\n","\n","#%% For filtered heatmap; specificity > 0.7\n","\n","cluster_df = cluster_df[cluster_df.specificity > 0.7]\n","cluster_df.to_csv(mpo_path + \"mapman_spec_df_07.txt\", sep=\"\\t\", index=None)\n","mapbin2_list = [x.split(\".\")[1] for x in cluster_df[\"index\"].unique().tolist()]\n","\n","# initialise dataframe\n","plot_spec_df = pd.DataFrame(np.nan, index = mapbin2_list, columns = stresses)\n","plot_spec_dir_df = pd.DataFrame(np.nan, index = mapbin2_list, columns = stresses)\n","# update dataframe\n","cluster_df.apply(lambda x: update_df(df=plot_spec_df, idx=x[\"index\"].split(\".\")[1], col=x.cluster[0], val=x.specificity), axis=1)\n","cluster_df.apply(lambda x: update_df(df=plot_spec_dir_df, idx=x[\"index\"].split(\".\")[1], col=x.cluster[0], val=x.top_direction), axis=1)\n","# change values to suit plotting\n","plot_spec_dir_df.replace({\"DOWN\": -1, \"UP_DOWN\": 0, \"UP\": 1}, inplace=True)\n","plot_dendro = plot_spec_dir_df * plot_spec_df\n","plot_dendro.fillna(0.00000000001, inplace=True)\n","plot_dendro.dropna(how=\"all\", inplace=True)\n","\n","# Plot\n","sns.clustermap(plot_dendro, cmap=\"coolwarm\",\n","\t\t\t   yticklabels=True, mask = plot_dendro == 0.00000000001, figsize=(6, 8),\n","\t\t\t   cbar_pos=(0.02, 0.9, 0.05, 0.08), dendrogram_ratio=(0.2,0.1)\n","\t\t\t   )\n","plt.savefig(dir_path + 'figures/' + 'FigS13.png', dpi=600)"],"metadata":{"id":"dSoyfrJtubic"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 4 B,C,D; Supp. Fig. 14, 15 Bipartite networks for robustly expressed TFs and biological processes\n","(Fig 4C, D and Supp. Fig. 15 is processed here but visualised in cytoscape)"],"metadata":{"id":"51xf_0EgCXyH"}},{"cell_type":"code","source":["# Preparation for Fig 4C\n","import pandas as pd\n","from ast import literal_eval\n","from collections import defaultdict, Counter\n","import random\n","from statsmodels.stats.multitest import multipletests\n","import numpy as np\n","\n","mapman_spec07_path = mpo_path + \"mapman_spec_df_07.txt\"\n","mapman_spec07_df = pd.read_csv(mapman_spec07_path, header=0, sep=\"\\t\")\n","\n","spec_lvl1bin = [x.split(\".\")[0] for x in mapman_spec07_df[\"index\"].to_list()]\n","\n","#%% mercator stuff\n","merdict = literal_eval(open(dir_path + 'mercator/merdict.txt', \"r\").read())\n","merbin = defaultdict(list)\n","map2anno = {}\n","\n","mfile = dir_path + 'mercator/results/MpoProt.results.txt'\n","merfile = open(mfile, 'r')\n","merfile.readline()\n","for line in merfile:\n","\tlinecon = line.rstrip().replace(\"'\", \"\").split(\"\\t\")\n","\tif len(linecon) == 5:\n","\t\tbincode, name, identifier, desc, ptype = linecon\n","\t\tmerbin[identifier].append('.'.join(name.split('.')[:2]))\n","\tif len(linecon[0].split('.')) == 2:\n","\t\tmap2anno[linecon[0]] = linecon[1]\n","\n","#%% TF-network enrichment (2nd level mapman)\n","nw_path = elnet_dir + \"union_0.8_ignoreAGRIS_topTF.txt\"\n","nw_df = pd.read_csv(nw_path, header=0, sep=\"\\t\")\n","nw_df[\"predicted_map2\"] = nw_df.apply(lambda x: merbin[x.predicted.lower()], axis=1)\n","\n","TF_list = nw_df[\"Gene.ID\"].unique().tolist()\n","anno_list = list(merbin.values())\n","sig_df = pd.DataFrame(0, index=map2anno.values(), columns=TF_list)\n","\n","for tf in TF_list:\n","\tnbh = nw_df[nw_df[\"Gene.ID\"] == tf]\n","\ts_count = Counter([y for x in nbh.predicted_map2.to_list() for y in x])\n","\tvalid_bins = list(s_count.keys())\n","\tsig_df[tf] = [1 if x in valid_bins else np.nan for x in sig_df.index.to_list()]\n","\t\n","\tfor i in range(1000):\n","\t\tsub = random.sample(anno_list, len(nbh))\n","\t\tsub_count = Counter([y for x in sub for y in x])\n","\t\tfor mapman in valid_bins:\n","\t\t\tif sub_count[mapman] >= s_count[mapman]:\n","\t\t\t\tsig_df.loc[mapman, tf] += 1\n","\n","sig_df.replace(1001, 1000, inplace=True)\n","sig_df /= 1000\n","\n","#%% BH correction\n","sig_df_corr = sig_df.copy()\n","df_idx = sig_df_corr.index.to_list()\n","for tf in TF_list:\n","\tpval_series = sig_df[tf].dropna()\n","\tpval_index = pval_series.index.tolist()\n","\tpval_list = pval_series.to_list()\n","\n","\ty = multipletests(pvals=pval_list, alpha=0.05, method=\"fdr_bh\")[1]\n","\tsig_df_corr[tf] = [y[pval_index.index(x)] if x in pval_index else np.nan for x in df_idx]\n","\t\n","#%% Flatten\n","sig_bin_flat = pd.DataFrame(columns=[\"TF\", \"map2bin\", \"corr_pval\"])\n","all_bin_flat = pd.DataFrame(columns=[\"TF\", \"map2bin\", \"corr_pval\"])\n","for tf in TF_list:\n","\tfor mapbin in df_idx:\n","\t\tpval = sig_df_corr.loc[mapbin, tf]\n","\t\tif not np.isnan(pval):\n","\t\t\tall_bin_flat.loc[len(all_bin_flat)] = [tf, mapbin, pval]\n","\t\tif pval < 0.05:\n","\t\t\tsig_bin_flat.loc[len(sig_bin_flat)] = [tf, mapbin, pval]\n","\n","sig_bin_flat.to_csv(elnet_dir + 'TF_map2bin_enrich_sig.txt', sep=\"\\t\", index=False)\n","all_bin_flat.to_csv(elnet_dir + 'TF_map2bin_enrich_all.txt', sep=\"\\t\", index=False)"],"metadata":{"id":"lZAxuwbWwHlR"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Fig 4B\n","import pandas as pd\n","from ast import literal_eval\n","from collections import defaultdict, Counter\n","import matplotlib.pyplot as plt\n","\n","TF_spec_path = elnet_dir + 'TF_Scond_specificity_07.txt'\n","TF_spec = pd.read_csv(TF_spec_path, header=0, sep=\"\\t\")\n","\n","mapman_spec_path = mpo_path + 'mapman_spec_df_07.txt'\n","mapman_spec = pd.read_csv(mapman_spec_path, header=0, sep=\"\\t\")\n","mapman_spec[\"top_dir_val\"] = mapman_spec.top_direction.replace({\"DOWN\": -1, \"UP_DOWN\": 0, \"UP\": 1})\n","mapman_spec[\"corr_specificity\"] = mapman_spec.specificity * mapman_spec.top_dir_val\n","\n","enrich_path = elnet_dir + 'TF_map2bin_enrich_sig.txt'\n","enriched_df = pd.read_csv(enrich_path, header=0, sep=\"\\t\")\n","\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","\n","#%% TF-enriched map2 (For Supp. Fig. 14 and 15)\n","for s in stresses:\n","\tTF_subset = TF_spec[TF_spec[s] != 0]\n","\tmapman_subset = mapman_spec[[x.startswith(s) for x in mapman_spec.cluster.to_list()]]\n","\t\n","\tspec_TF_long = [x for x in TF_subset.iloc[:,0].to_list()]\n","\tspec_TF = [x.split(\" (\")[0] for x in spec_TF_long]\n","\tspec_mapman = mapman_subset[\"index\"].to_list()\n","\t\n","\tenriched_subset = enriched_df[(enriched_df.TF.isin(spec_TF)) & (enriched_df.map2bin.isin(spec_mapman))]\n","\tif len(enriched_subset) != 0:\n","\t\tenriched_subset[\"TF_spec\"] = enriched_subset.apply(lambda x: TF_subset[TF_subset.iloc[:,0] == spec_TF_long[spec_TF.index(x.TF)]][s].values.tolist()[0], axis=1)\n","\t\tenriched_subset[\"mapman_spec\"] = enriched_subset.apply(lambda x: mapman_subset[(mapman_subset['index'] == x.map2bin) & ([x.startswith(s) for x in mapman_subset.cluster.to_list()])].corr_specificity.values.tolist()[0], axis=1)\n","\t\tenriched_subset[\"TF_dir\"] = enriched_subset.apply(lambda x: 'UP' if x.TF_spec > 0 else \"DOWN\", axis=1)\n","\t\tenriched_subset[\"mapman_dir\"] = enriched_subset.apply(lambda x: 'UP' if x.mapman_spec > 0 else \"DOWN\", axis=1)\n","\t\tenriched_subset.to_csv(elnet_dir + 'TF_map2_spec_bipartite_' + s + '.txt', sep=\"\\t\", index=False)\n","\t\n","enriched_all = enriched_df[(enriched_df.TF.isin([x.split(\" (\")[0] for x in TF_spec.iloc[:,0].to_list()])) & (enriched_df.map2bin.isin(mapman_spec[\"index\"].to_list()))]\n","\n","#%% TF-map2 (ratio of map2/all genes in map2)\n","\n","# mapman stuff\n","merdict = literal_eval(open(dir_path + 'mercator/merdict.txt', \"r\").read())\n","merbin = defaultdict(list)\n","map2anno = {}\n","\n","mfile = dir_path + 'mercator/MpoProt.results.txt'\n","merfile = open(mfile, 'r')\n","merfile.readline()\n","for line in merfile:\n","\tlinecon = line.rstrip().replace(\"'\", \"\").split(\"\\t\")\n","\tif len(linecon) == 5:\n","\t\tbincode, name, identifier, desc, ptype = linecon\n","\t\tmerbin[identifier].append('.'.join(name.split('.')[:2]))\n","\tif len(linecon[0].split('.')) == 2:\n","\t\tmap2anno[linecon[0]] = linecon[1]\n","\n","bin_size = Counter([y for x in list(merbin.values()) for y in x])\n","\n","# nw stuff\n","nw_path = elnet_dir + \"union_0.8_ignoreAGRIS_topTF.txt\"\n","nw_df = pd.read_csv(nw_path, header=0, sep=\"\\t\")\n","nw_df[\"predicted_map2\"] = nw_df.apply(lambda x: merbin[x.predicted.lower()], axis=1)\n","\n","#%% tf_stat anno\n","\n","spec_path = elnet_dir + 'TF_Scond_spec.txt'\n","spec_df = pd.read_csv(spec_path, header=0, sep=\"\\t\")\n","tf_unique = spec_df.drop_duplicates(subset=['gene'])[[\"gene\", \"TF_anno\"]]\n","\n","elnet_path = elnet_dir + \"union_0.8_ignoreAGRIS_topTF.txt'\n","elnet_df = pd.read_csv(elnet_path, header=0, sep=\"\\t\")\n","\n","def get_TFstat(tf):\n","\tstat_dict = Counter(elnet_df[elnet_df[\"Gene.ID\"] == tf].TF_stat)\n","\tTF_stat_max = max(stat_dict, key=stat_dict.get)\n","\tif TF_stat_max == 'Activator':\n","\t\treturn(\"limegreen\")\n","\telif TF_stat_max == 'Ambiguous':\n","\t\treturn(\"yellow\")\n","\telse:\n","\t\treturn(\"firebrick\")\n","tf_unique[\"stat_colour\"] = tf_unique.apply(lambda x: get_TFstat(x.gene), axis=1)\n","tf_unique[\"top_stat\"] = tf_unique.apply(lambda x: \"Activator\" if x.stat_colour == \"limegreen\" else \"Ambiguous\" if x.stat_colour == \"yellow\" else \"Repressor\", axis=1)\n","\n","tf_unique.to_csv(elnet_dir + 'TF_top_stat_anno.txt', sep=\"\\t\", index=False)\n","\n","#%% calculate\n","df_collate = []\n","for s in stresses:\n","\tcounter = pd.DataFrame(columns=[\"TF\", \"Mapman\", \"TF_specificity\", \"Mapman_specificity\", \"Mapman_ratio\", \"TF_type\"])\n","\tTF_subset = TF_spec[TF_spec[s] != 0][[\"Unnamed: 0\", s]]\n","\tmapman_subset = mapman_spec[[x.startswith(s) for x in mapman_spec.cluster.to_list()]]\n","\tfor i, tf in enumerate([x.split(\" (\")[0] for x in TF_subset.iloc[:,0].to_list()]):\n","\t\ttfsub = nw_df[nw_df[\"Gene.ID\"] == tf]\n","\t\ttf_spec = TF_subset.iloc[i, 1]\n","\t\ttf_type = tf_unique[tf_unique.gene == tf].top_stat.to_list()[0]\n","\t\tbin_count = Counter([y for x in tfsub.predicted_map2.to_list() for y in x])\n","\t\tbin_norm = {k: v/bin_size[k] for k, v in bin_count.items()}\n","\t\tfor j, m in enumerate(mapman_subset[\"index\"].to_list()):\n","\t\t\tif m in bin_norm:\n","\t\t\t\tmap2_spec = mapman_subset.iloc[j, 5]\n","\t\t\t\tmap2_ratio = bin_norm[m]\n","\t\t\t\tcounter.loc[len(counter)] = [tf, m, tf_spec, map2_spec, map2_ratio, tf_type]\n","\t# bipartite graphs\n","    counter.to_csv(elnet_dir + 'TF_map2_spec_bipartite_ratio_' + s + '.txt', sep=\"\\t\", index=False)\n","\tcounter['stress'] = s\n","\tdf_collate.append(counter)\n","\tTF_unique = counter.drop_duplicates(subset=['TF'])\n","\tmapman_unique = counter.drop_duplicates(subset=['Mapman'])\n","\tnode = TF_unique.TF.to_list() + mapman_unique.Mapman.to_list()\n","\tspecificity = TF_unique.TF_specificity.to_list() + mapman_unique.Mapman_specificity.to_list()\n","\tcombined = [\"\\t\".join([node[i], str(specificity[i])]) + \"\\n\" for i in range(len(node))]\n","\t# annotation for graphs\n","    with open(elnet_dir + 'TF_map2_spec_bipartite_ratio_anno' + s + '.txt', \"w+\") as annof:\n","\t\tfor line in combined:\n","\t\t\tannof.write(line)\n","\n","df_all = pd.concat(df_collate)\n","df_cut05 = df_all[df_all.Mapman_ratio >= 0.05]\n","\n","df_all.Mapman_ratio.plot.hist()\n","plt.xlabel(\"Number of genes regulated by TFx in binY/ Total number of genes in binY\")\n","plt.savefig(dir_path + 'figures/' + 'FigS14.png', dpi=600)\n","\n","def expected_outcome(tf_spec, mapman_spec, tf_type):\n","\tif tf_spec > 0:\n","\t\tif mapman_spec > 0:\n","\t\t\tif tf_type == \"Activator\":\n","\t\t\t\treturn \"Y\"\n","\t\t\telse:\n","\t\t\t\treturn \"N\"\n","\t\telif mapman_spec < 0:\n","\t\t\tif tf_type == \"Repressor\":\n","\t\t\t\treturn \"Y\"\n","\t\t\telse:\n","\t\t\t\treturn \"N\"\n","\t\telse:\n","\t\t\treturn(\"U\")\n","\telif tf_spec < 0:\n","\t\tif mapman_spec < 0:\n","\t\t\tif tf_type == \"Activator\":\n","\t\t\t\treturn \"Y\"\n","\t\t\telse:\n","\t\t\t\treturn \"N\"\n","\t\telif mapman_spec > 0:\n","\t\t\tif tf_type == \"Repressor\":\n","\t\t\t\treturn \"Y\"\n","\t\t\telse:\n","\t\t\t\treturn \"N\"\n","\t\telse:\n","\t\t\treturn(\"U\")\n","\n","df_cut05[\"expected\"] = df_cut05.apply(lambda x: expected_outcome(x.TF_specificity, x.Mapman_specificity, x.TF_type) , axis=1)\n","count = Counter(df_cut05.expected)\n","plt.pie(count.values(), labels=[\"Expected\" if x == \"Y\" else \"Not expected\" if x ==  \"N\" else \"Ambiguous\" for x in count.keys()])\n","count_vals = list(count.values())[::-1]\n","count_keys = [\"Expected\" if x == \"Y\" else \"Not expected\" if x ==  \"N\" else \"Ambiguous\" for x in count.keys()][::-1]\n","bottom_list = [0, count_vals[0], count_vals[1] + count_vals[0]]\n","# Fig 4B\n","fig, ax = plt.subplots(figsize=(1,6))\n","for i, x in enumerate(count_vals):\n","\tax.bar(\"relationship\", x, width = 0.3, label=count_keys[i], bottom=bottom_list[i])\n","plt.ylabel(\"Count\")\n","ax.legend(bbox_to_anchor=(1.3, -0.06))\n","\n","plt.savefig(dir_path + 'figures/' + 'Fig4B.png', bbox_inches='tight', dpi=600)"],"metadata":{"id":"V6kL2bmVCegm"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 5: Annotation of Ath orthologs with evidence from literature"],"metadata":{"id":"AtxfxK3uKxuG"}},{"cell_type":"code","source":["import pandas as pd\n","import numpy as np\n","import matplotlib.pyplot as plt\n","from matplotlib.colors import ListedColormap\n","import seaborn as sns\n","from matplotlib import colors\n","\t\t\n","hm_path = dir_path + 'prep_files/anno_hm_compat.txt'\n","hm_df = pd.read_csv(hm_path, header=[0,1], index_col=[0,1], sep=\"\\t\")\n","hm_df.replace(np.nan, 0, inplace=True)\n","hm_df.replace(\"L\", 2, inplace=True)\n","hm_df.replace(\"GO\", 2, inplace=True)\n","hm_df.replace(\"T\", 1, inplace=True)\n","\n","xlabels = [x[1].replace(\"\\r\",\"\") for x in hm_df.columns.to_list()[1:]]\n","xlabels = [x.replace(\"Nitrogen\", \"Nitrogen\\ndefeciency\") for x in xlabels]\n","\n","ylabels = [x[0] if type(x[0]) != float else \"No specific response\" for x in hm_df.index.to_list()]\n","ylabel_unique = []\n","ylabel_pos = []\n","\n","for lab in list(set(ylabels)):\n","\tylabel_unique.append(lab)\n","\tpos = [i for i, x in enumerate(ylabels) if x == lab]\n","\tif pos == 1:\n","\t\tylabel_pos.append(pos)\n","\telse:\n","\t\tylabel_pos.append((pos[0] + pos[-1])/2)\n","ylabel_unique_mod = [x.replace(\"Light/Nitrogen\", \"           Light/Nitrogen\") for x in ylabel_unique]\n","ylabel_unique_mod = [x.replace(\"Cold/Nitrogen\", \"           Cold/Nitrogen\") for x in ylabel_unique_mod]\n","\n","yunqiue_ordered = []\n","for l in ylabels:\n","\tif l not in yunqiue_ordered:\n","\t\tyunqiue_ordered.append(l)\n","\n","#%% plotting\n","\n","row_col_list = [\"orange\", \"lemonchiffon\", \"rosybrown\",\n","\t\t\t\t\"peachpuff\", \"orange\", \"lemonchiffon\",\n","\t\t\t\t\"rosybrown\", \"peachpuff\", \"orange\",\n","\t\t\t\t\"lemonchiffon\", \"rosybrown\", \"peachpuff\",\n","\t\t\t\t\"orange\", \"lemonchiffon\", \"gray\"]\n","\n","row_col_list = [colors.CSS4_COLORS[x] for x in row_col_list]\n","yrowcol = [row_col_list[yunqiue_ordered.index(y)] for y in ylabels]\n","\n","colour_list = ListedColormap([\"white\", \"royalblue\", \"limegreen\"])\n","g = sns.clustermap(hm_df.iloc[:,1:].values, cmap=colour_list,\n","\t\t\t\t   col_cluster=False, row_cluster=False,\n","\t\t\t\t   cbar_kws={\"ticks\": [1/3, 3/3, 5/3], },\n","\t\t\t\t   figsize=(6,10), row_colors=yrowcol)\n","\n","g.ax_heatmap.axes.set_yticks(ylabel_pos, ylabel_unique_mod)\n","g.ax_heatmap.axes.set_xticks([i+0.5 for i in range(len(xlabels))], xlabels, rotation=90)\n","g.cax.set_yticklabels(['NA', 'Observed', 'Literature'])\n","plt.savefig(dir_path + 'figures/' + 'Fig5A.png', dpi=600, bbox_inches=\"tight\")\n","\n","#%% stacked bar by row\n","stresses_lit = hm_df.columns.to_list()[1:8]\n","stresses_obs = hm_df.columns.to_list()[8:]\n","stress_lit = [\"Heat\", \"Cold\", \"Salt\", \"Mannitol\", \"Light\", \"Dark\", \"Nitrogen\", \"nan\"]\n","stress_obs = [\"Heat\", \"Cold\", \"Salt\", \"Mannitol\", \"nan\"]\n","\n","def get_counts(stresses, stress):\n","\tcoll_df = pd.concat([hm_df[[s in str(x[0]) for x in hm_df.index.to_list()]][stresses].sum() for s in stress], axis=1)\n","\tcoll_df.columns = stress[:-1] + [\"No specific response\"]\n","\tcoll_df.index = [x[1].replace(\"\\r\\n\", \"\") for x in coll_df.index.to_list()]\n","\tcoll_df_per = pd.DataFrame(index=coll_df.index.to_list())\n","\tlit_sum = coll_df.sum()\n","\tfor col in coll_df.columns.to_list():\n","\t\tcoll_df_per[col] = coll_df[col].apply(lambda x: x/lit_sum[col] if lit_sum[col] > 0 else 0)\n","\treturn(coll_df_per)\n","\n","# Plotting\n","lit_col = colors.ListedColormap([\"steelblue\", \"darkorange\", \"mediumseagreen\", \"crimson\", \"mediumpurple\",\"saddlebrown\", \"hotpink\"])\n","obs_col = colors.ListedColormap([\"steelblue\", \"darkorange\", \"mediumpurple\",\"saddlebrown\"])\n","\n","col_lit = get_counts(stresses_lit, stress_lit)\n","col_obs = get_counts(stresses_obs, stress_obs)\n","col_lit.T.plot.bar(stacked=True, cmap=lit_col)\n","plt.legend(bbox_to_anchor = (1., 1.))\n","plt.ylabel(\"Ratio of Arabidopsis responses\")\n","plt.xlabel(\"Response in Marchantia\")\n","plt.savefig(dir_path + 'figures/' + 'Fig5B.png', dpi=600, bbox_inches=\"tight\")\n","\n","col_obs.T.plot.bar(stacked=True, cmap=obs_col, figsize=(4.9,5.2), fontsize=13)\n","plt.legend(bbox_to_anchor = (1., 1.))\n","plt.ylabel(\"Ratio of Arabidopsis responses\", fontsize=13)\n","plt.xlabel(\"Response in Marchantia\", fontsize=13)\n","plt.savefig(dir_path + 'figures/' + 'Fig5C.png', dpi=600, bbox_inches=\"tight\")\n","\n","#%% for supp table S8\n","for stress in [stress_lit, stress_obs]:\n","\tcoll_df = pd.concat([hm_df[[s in str(x[0]) for x in hm_df.index.to_list()]][stresses].sum() for s in stress], axis=1)\n","\tcoll_df.columns = stress[:-1] + [\"No specific response\"]\n","\tcoll_df.index = [x[1].replace(\"\\r\\n\", \"\") for x in coll_df.index.to_list()]\n","\tformatted = coll_df.T\n","\tformatted[\"Number of corresponding orthologs\"] = [len(hm_df[[s in str(x[0]) for x in hm_df.index.to_list()]]) for s in stress]\n","\tif len(stress) == 8:\n","\t\tfname = \"lit_count.txt\"\n","\telse:\n","\t\tfname = \"obs_count.txt\"\n","\tformatted.to_csv(dir_path + '/prep_files/' + fname, sep=\"\\t\")\n"],"metadata":{"id":"GLMMOHUoK5UM"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 6 A, B: Effects of combined stress in terms of significant L2FC"],"metadata":{"id":"zvufjLWgMPYJ"}},{"cell_type":"code","source":["# Figure 6A\n","import pandas as pd\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","from matplotlib import colors\n","from matplotlib import cm\n","from matplotlib.colors import ListedColormap\n","\n","#%% intialisation\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","stresses_long = [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen deficiency\"]\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN',\n","\t\t\t\t'N', 'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","combined_stress = [x for x in all_stresses if len(x) == 2]\n","\n","#%% get averaged L2FC files\n","deseq_dir = mpo_path\n","def deg_stat(H2, D2):\n","\tif H2 > 0 and D2 > 0:\n","\t\treturn(\"UP\")\n","\telif H2 < 0 and D2 < 0:\n","\t\treturn(\"DOWN\")\n","\telse:\n","\t\treturn(\"MIXED\")\n","\t\n","def get_resSig(stress):\n","\tH2df = pd.read_csv(deseq_dir + stress + \"controlH2_resSig.tsv\", header=0, sep=\"\\t\")\n","\tD2df = pd.read_csv(deseq_dir + stress + \"controlD2_resSig.tsv\", header=0, sep=\"\\t\")\n","\tH2df.columns = ['gene', 'baseMean', 'L2FC_H2', 'lfcSE', 'stat', 'pvalue', 'padj']\n","\tD2df.columns = ['gene', 'baseMean', 'L2FC_D2', 'lfcSE', 'stat', 'pvalue', 'padj']\n","\tmerged = pd.merge(H2df[['gene', 'L2FC_H2']], D2df[['gene', 'L2FC_D2']], on='gene', how='outer')\n","\tmerged.dropna(inplace=True)\n","\tmerged['deg_stat'] = merged.apply(lambda x: deg_stat(x.L2FC_H2, x.L2FC_D2) , axis=1)\n","\tmerged_filtered = merged[merged.deg_stat != \"MIXED\"]\n","\tmerged_filtered[\"L2FC_avg\"] = merged_filtered.apply(lambda x: sum([x.L2FC_H2, x.L2FC_D2])/2, axis=1)\n","\tmerged_filtered[\"stress\"] = stress\n","\treturn merged_filtered\n","\n","resSig_list = [get_resSig(x) for x in all_stresses]\n","resSig_df = pd.concat(resSig_list)\n","resSig_df.to_csv(deseq_dir + \"resSig_L2FC_compiled.txt\", sep=\"\\t\", index = False)\n","\n","#%% get stats for each gene\n","cds_path = dir_path + 'prep_files/Mpo.cds.fasta'\n","all_genes = [x[1:-1] for x in open(cds_path, \"r\").readlines() if \">\" in x]\n","\n","deg_path = mpo_path + 'resSig_L2FC_compiled.txt'\n","deg_df = pd.read_csv(deg_path, header=0, sep=\"\\t\")\n","\n","int_class_df = pd.DataFrame(columns = [\"gene\", \"Xname\", \"Yname\", \"XYname\", \"Xl2fc\", \"Yl2fc\", \"XYl2fc\"])\n","\n","for sxy in combined_stress:\n","\tprint(\"Calculating: \" + sxy)\n","\tsx = sxy[0]\n","\tsy = sxy[1]\n","\tfor g in all_genes:\n","\t\tstat_series = [deg_df[(deg_df.stress == x) & (deg_df.gene == g)].L2FC_avg for x in [sx, sy, sxy]]\n","\t\tstat_list = [\"NC\" if len(x) == 0 else x.to_list()[0] for x in stat_series]\n","\t\tint_class_df.loc[len(int_class_df)] = [g, sx, sy, sxy] + stat_list\n","\n","int_class_df.to_csv(mpo_path + 'stress_int_class_l2fc.txt', sep=\"\\t\", index=None)\n","\n","# get reverse of XY i.e. YX for the sake of counting\n","int_class_rev_df = int_class_df[[\"gene\", \"Yname\", \"Xname\", \"XYname\", \"Yl2fc\", \"Xl2fc\", \"XYl2fc\"]]\n","int_class_rev_df.columns = [\"gene\", \"Xname\", \"Yname\", \"XYname\", \"Xl2fc\", \"Yl2fc\", \"XYl2fc\"]\n","\n","int_class_all_df = pd.concat([int_class_df, int_class_rev_df])\n","int_class_all_df[\"Xstat\"] = [x if x == \"NC\" else \"UP\" if float(x) > 0 else \"DOWN\" for x in int_class_all_df.Xl2fc.to_list()]\n","int_class_all_df[\"Ystat\"] = [x if x == \"NC\" else \"UP\" if float(x) > 0 else \"DOWN\" for x in int_class_all_df.Yl2fc.to_list()]\n","int_class_all_df[\"XYstat\"] = [x if x == \"NC\" else \"UP\" if float(x) > 0 else \"DOWN\" for x in int_class_all_df.XYl2fc.to_list()]\n","int_class_all_df[\"Xstat_Ystat\"] = int_class_all_df.Xstat + \"_\" + int_class_all_df.Ystat\n","int_class_all_df[\"all_stat\"] = int_class_all_df[\"Xstat_Ystat\"] + \"_\" + int_class_all_df[\"XYstat\"]\n","\n","\n","for col in [\"Xl2fc\", \"Yl2fc\", \"XYl2fc\"]:\n","\tint_class_all_df[col] = int_class_all_df[col].replace(\"NC\", 0)\n","\tint_class_all_df[col] = int_class_all_df[col].apply(float)\n","\n","#%% all_stat\n","all_reshaped = int_class_all_df.groupby([\"Xname\", \"all_stat\"])[[\"Xl2fc\", \"Yl2fc\", \"XYl2fc\"]].mean()\n","s_reshaped = [all_reshaped.xs(x, level=\"Xname\") for x in stresses]\n","merged_flat = pd.concat(s_reshaped, axis=1)\n","\n","col1 = [x for x in stresses for i in range(3)]\n","col2 = merged_flat.columns.to_list()\n","index = pd.MultiIndex.from_tuples(zip(col1, col2), names=[\"Xname\", \"l2fc\"])\n","\n","merged_flat.columns = index\n","\n","plt.figure(figsize=(6, 7))\n","sns.heatmap(merged_flat, cmap=\"coolwarm\", yticklabels=True)\n","\n","#%% Xstat_Ystat\n","\n","two_reshaped = int_class_all_df.groupby([\"Xname\", \"Xstat_Ystat\"])[[\"Xl2fc\", \"Yl2fc\", \"XYl2fc\"]].mean()\n","s_reshaped = [two_reshaped.xs(x, level=\"Xname\") for x in stresses]\n","merged_flat = pd.concat(s_reshaped, axis=1)\n","\n","col1 = [x for x in stresses for i in range(3)]\n","col2 = merged_flat.columns.to_list()\n","index = pd.MultiIndex.from_tuples(zip(col1, col2), names=[\"Xname\", \"l2fc\"])\n","\n","merged_flat.columns = index\n","\n","merged_flat.to_csv(mpo_path + 'stress_int_class_l2fc_avg.txt', sep=\"\\t\")\n","\n","plt.figure(figsize=(11, 2.3))\n","sns.heatmap(merged_flat, cmap=\"coolwarm\", yticklabels=True, annot=merged_flat, fmt=\",.1f\")\n","plt.savefig(dir_path + 'figures/' + 'Fig6A.png', dpi=600)"],"metadata":{"id":"zVajCOu5MtAE"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#Figure 6B: Scatter plot\n","\n","int_class_all_df = pd.read_csv(mpo_path + 'stress_int_class_l2fc_extended.txt', sep=\"\\t\", header=0)\n","# calculate difference between Sxy and Sx or Sy (depending if it is higher than both or lower than both)\n","int_class_all_df[\"XY_val\"] = int_class_all_df.apply(lambda x: x.XYl2fc - max(x.Xl2fc, x.Yl2fc) if x.XYl2fc > x.Xl2fc and x.XYl2fc > x.Yl2fc else x.XYl2fc - min(x.Xl2fc, x.Yl2fc) if x.XYl2fc < x.Xl2fc and x.XYl2fc < x.Yl2fc else 0, axis=1)\n","int_class_all_df.to_csv(mpo_path + 'stress_int_class_l2fc_extended.txt', sep=\"\\t\", index=None)\n","\n","# Set limits for scatter plot for general overview\n","int_class_all_df_sub = int_class_all_df[(abs(int_class_all_df.Xl2fc) < 20) & (abs(int_class_all_df.Yl2fc) < 20)]\n","int_class_all_df_sub[\"XY_relstat\"] = int_class_all_df_sub.apply(lambda x: \"firebrick\" if x.XY_val > 0 else \"royalblue\" if x.XY_val < 0 else \"lightgray\", axis=1)\n","\n","#%% plot\n","int_class_all_df_sub.plot.scatter(x=\"Xl2fc\", y=\"Yl2fc\", marker='.',\n","\t\t\t\t\t\t\t\t  c=\"XY_relstat\", linewidths=0, alpha=0.2,\n","\t\t\t\t\t\t\t\t  xlim=(-12,12), ylim=(-12,12))\n","plt.savefig(dir_path + 'figures/' + 'Fig6B.png', dpi=600)"],"metadata":{"id":"Pj2zY_uFNd8X"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 6C: Classification of stress interactions"],"metadata":{"id":"UFx24ozgPE-T"}},{"cell_type":"code","source":["\"\"\"\n","Classification of interactions\n","\"\"\"\n","\n","import pandas as pd\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","stresses_long = [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen deficiency\"]\n","all_stresses = ['C', 'CD', 'CL', 'CM', 'CN', 'CS', 'D', 'H', 'HD',\n","\t\t\t\t'HM', 'HN', 'HS', 'L', 'LS', 'M', 'MD', 'ML','MN',\n","\t\t\t\t'N', 'ND', 'NL', 'S', 'SD', 'SM', 'SN']\n","combined_stress = [x for x in all_stresses if len(x) == 2]\n","\n","cds_path = dir_path + 'prep_files/Mpo.cds.fasta'\n","all_genes = [x[1:-1] for x in open(cds_path, \"r\").readlines() if \">\" in x]\n","deg_path = mpo_path + 'resSig_compiled.txt'\n","deg_df = pd.read_csv(deg_path, header=0, sep=\"\\t\")\n","int_class_df = pd.DataFrame(columns = [\"gene\", \"Xname\", \"Yname\", \"XYname\", \"Xstat\", \"Ystat\", \"XYstat\"])\n","\n","# get status (up-, down-regulated and no change) for each gene of stress X, Y and XY \n","for sxy in combined_stress:\n","\tprint(\"Calculating: \" + sxy)\n","\tsx = sxy[0]\n","\tsy = sxy[1]\n","\tfor g in all_genes:\n","\t\tstat_series = [deg_df[(deg_df.stress == \"Mpo_\" + x) & (deg_df.gene == g)].L2FC_D2 for x in [sx, sy, sxy]]\n","\t\tstat_list = [\"NC\" if len(x) == 0 else x.to_list()[0] for x in stat_series]\n","\t\tint_class_df.loc[len(int_class_df)] = [g, sx, sy, sxy] + stat_list\n","\n","int_class_df.to_csv(mpo_path + 'stress_int_classification.txt', sep=\"\\t\", index=None)\n","\n","# get reverse of XY i.e. YX for the sake of counting\n","int_class_rev_df = int_class_df[[\"gene\", \"Yname\", \"Xname\", \"XYname\", \"Ystat\", \"Xstat\", \"XYstat\"]]\n","int_class_rev_df.columns = [\"gene\", \"Xname\", \"Yname\", \"XYname\", \"Xstat\", \"Ystat\", \"XYstat\"]\n","\n","int_class_all_df = pd.concat([int_class_df, int_class_rev_df])\n","int_class_all_df[\"Xstat_Ystat\"] = int_class_all_df.Xstat + \"_\" + int_class_all_df.Ystat\n","\n","int_class_all_df.to_csv(mpo_path + 'stress_int_classification_extended.txt', sep=\"\\t\", index=None)\n","\n","int_class_count = int_class_all_df.groupby([\"Xname\", \"Xstat_Ystat\", \"XYstat\"]).count()[\"Yname\"]\n","int_class_count.to_csv(mpo_path + 'stress_int_classification_count.txt', sep=\"\\t\")\n","\n","#%% For collated plot, separate norm\n","int_class_all_df[\"all_stat\"] = int_class_all_df.Xstat + \"_\" + int_class_all_df.Ystat + \"_\" + int_class_all_df.XYstat\n","collapsed_count = int_class_all_df.groupby([\"Xname\", \"all_stat\"]).count()[\"Yname\"]\n","all_unstack = collapsed_count.unstack()\n","\n","#%% Plotting\n","fig, ax = plt.subplots(9, 1, sharex=True, constrained_layout=True, figsize=(8,10))\n","axes = ax.flatten()\n","\n","for i in range(0,27,3):\n","\tsubset = all_unstack.iloc[:,i:i+3]\n","\ttotal_int_stress = subset.sum(axis=1)\n","\tall_unstack_norm = subset.apply(lambda x: x/total_int_stress)\n","\tsubset = subset.T\n","\tall_unstack_norm = all_unstack_norm.T\n","\tsns.heatmap(all_unstack_norm, cmap=\"Blues\", annot=subset, fmt=\",.0f\", ax=axes[int(i/3)])\n","for axis in range(8):\n","\taxes[axis].set_xlabel('')\n","for axis in range(9):\n","\taxes[axis].set_ylabel('')\n","plt.savefig(dir_path + 'figures/' + 'Fig6C.png', dpi=600)"],"metadata":{"id":"WME4DLfAPJtd"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 6D: Linear regression of all experiments\n","Note: package not compatible to colab, please run locally"],"metadata":{"id":"jKN2LWQ5NmHW"}},{"cell_type":"code","source":["# Figure 6D\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import sklearn.linear_model\n","import pandas as pd\n","from sklearn.utils import shuffle\n","from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error\n","import seaborn as sns\n","\n","# directory paths\n","deseq_dir = mpo_path\n","linreg_dir = deseq_dir + \"linreg/\"\n","linreg_path_safe = linreg_dir.replace(' ', '\\ ')\n","!mkdir $linreg_path_safe\n","\n","# data paths\n","avg_path = deseq_dir + \"stress_int_class_l2fc_avg.txt\"\n","norm_path = deseq_dir + \"stress_int_class_l2fc_allres.txt\"\n","extended_path = deseq_dir + \"stress_int_class_l2fc_extended_allres.txt\"\n","\n","#%% for averaged linreg\n","avg_df = pd.read_csv(avg_path, sep=\"\\t\", header=[0,1], index_col = 0)\n","stress_col = avg_df.columns.get_level_values(0).unique().tolist()\n","avg_flat = pd.concat([avg_df[x] for x in stress_col])\n","avg_flat_shuffled = shuffle(avg_flat)\n","\n","X_train = avg_flat_shuffled[[\"Xl2fc\", \"Yl2fc\"]].values\n","Y_train = avg_flat_shuffled[\"XYl2fc\"].values\n","\n","# Training\n","model = sklearn.linear_model.LinearRegression()\n","model.fit(X_train, Y_train)\n","\n","#print('Score R2:', model.score(xtrain, ytrain))\n","coefs = model.coef_\n","intercept = model.intercept_\n","\n","Y_pred = model.predict(X_train)\n","MAE = round(mean_absolute_error(Y_train, Y_pred), 2)\n","RMSE = round(mean_squared_error(Y_train, Y_pred, squared=False), 2)\n","r2 = round(r2_score(Y_train, Y_pred), 2)\n","\n","minx, maxx = avg_flat_shuffled.Xl2fc.min()-1, avg_flat_shuffled.Xl2fc.max()+1\n","miny, maxy = avg_flat_shuffled.Yl2fc.min()-1, avg_flat_shuffled.Yl2fc.max()+1\n","\n","# for mesh\n","xs = np.tile(np.arange(minx, maxx), (len(np.arange(miny, maxy)),1))\n","ys = np.tile(np.arange(miny, maxy), (len(np.arange(minx, maxx)),1)).T\n","zs = xs*coefs[0]+ys*coefs[1]+intercept\n","\n","# Plot\n","fig = plt.figure(constrained_layout=True)\n","ax = fig.add_subplot(111, projection='3d')\n","ax.set_xlabel(\"Sx\")\n","ax.set_ylabel(\"Sy\")\n","ax.set_zlabel(\"Sxy\")\n","\n","ax.plot_surface(xs,ys,zs, alpha=0.5)\n","ax.stem(X_train[:,0], X_train[:,1], Y_pred, bottom=-6,\n","\t\tbasefmt=\" \", markerfmt=\"C1o\", linefmt=\"C1-\")\n","eqn_str = \"Sxy = {:.2f} + {:.2f} Sx + {:.2f} Sy\".format(intercept, coefs[0], coefs[1])\n","err_str = \"MAE: {}    RMSE: {}   R\".format(MAE, RMSE) + u\"\\u00b2\" + \": {}\".format(r2)\n","\n","z_bound = ax.get_zbound()\n","x_lower_bound = ax.get_xbound()[0]\n","y_bound = ax.get_ybound()\n","z_bound_dist = z_bound[1] - z_bound[0]\n","ax.text(x=x_lower_bound - x_lower_bound/4, y=(y_bound[1] - y_bound[0])/2, z = z_bound[1] + z_bound_dist/6, s=eqn_str)\n","ax.text(x=x_lower_bound - x_lower_bound/4, y=(y_bound[1] - y_bound[0])/2, z = z_bound[1] + 0.7, s=err_str)\n","plt.savefig(dir_path + 'figures/' + 'Fig6D.png', dpi=600, bbox_inches=\"tight\")\n","plt.show()"],"metadata":{"id":"dUrZIMyTN2Li"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### Figure 7 A-H: Linear regression by stress"],"metadata":{"id":"JlaKpugiOhT4"}},{"cell_type":"code","source":["#%% Plotting\n","def plot_3d(X_test, Y_test, title, mae, rmse, r2, xs, ys, zs, coefs, intercept):\n","\tfig = plt.figure(constrained_layout=True)\n","\tax = fig.add_subplot(111, projection='3d')\n","\tax.set_xlabel(\"Sx\")\n","\tax.set_ylabel(\"Sy\")\n","\tax.set_zlabel(\"Sxy\")\n","\t#plt.ylim((-12,12))\n","\t#plt.xlim((-12,12))\n","\t\n","\tax.plot_surface(xs,ys,zs, alpha=0.5)\n","\tplane_zmin = min([y for x in zs.tolist() for y in x])\n","\tif plane_zmin < min(Y_test):\n","\t\tstem_bottom = plane_zmin\n","\telse:\n","\t\tstem_bottom = min(Y_test)\n","\t\n","\txvals = X_test[:,0]\n","\tyvals = X_test[:,1]\n","\tzvals = Y_test\n","\tfor i in range(len(xvals)):\n","\t\tax.plot([xvals[i], xvals[i]], [yvals[i], yvals[i]], [stem_bottom, zvals[i]], \n","\t\t  '-', linewidth=1, color='darkgray', alpha=.1)\n","\n","\t# plotting a circle on the top of each stem\n","\tax.plot(xvals, yvals, zvals, 'o', markersize=4, color='orange',label='ib', alpha=.1)\n","\teqn_str = \"Sxy = {:.2f} + {:.2f} Sx + {:.2f} Sy\".format(intercept, coefs[0], coefs[1])\n","\terr_str = \"MAE: {}    RMSE: {}   R\".format(mae, rmse) + u\"\\u00b2\" + \": {}\".format(r2)\n","\t\n","\tz_bound = ax.get_zbound()\n","\tx_lower_bound = ax.get_xbound()[0]\n","\ty_bound = ax.get_ybound()\n","\tz_bound_dist = z_bound[1] - z_bound[0]\n","\t# https://glowingpython.blogspot.com/2012/12/3d-stem-plot.html\n","\tax.text(x=x_lower_bound - x_lower_bound/4, y=(y_bound[1] - y_bound[0])/3, z = z_bound[1] + z_bound_dist/8, s=eqn_str)\n","\tax.text(x=x_lower_bound - x_lower_bound/4, y=(y_bound[1] - y_bound[0])/3, z = z_bound[1], s=err_str)\n","\tplt.savefig(dir_path + 'figures/' + \"Fig7_\" + title + \".png\", dpi=600, bbox_inches=\"tight\")\n","\tplt.show()\n","\n","#%% Process and plot\n","stresses = [\"C\", \"S\", \"M\", \"H\", \"D\", \"L\", \"N\"]\n","stresses_long = [\"Cold\", \"Salt\", \"Mannitol\", \"Heat\", \"Dark\", \"Light\", \"Nitrogen\\ndeficiency\"]\n","\n","data_df = pd.read_csv(extended_path, header=0, sep=\"\\t\")\n","sub_len = [len(data_df[data_df.Xname == x]) for x in stresses]\n","def process_data(data_type, idx, permu, j=None):\n","\ttitle = stresses_long[idx]\n","\tif permu == False:\n","\t\tshuffled_df = shuffle(data_df[data_df.Xname == data_type])\n","\telse:\n","\t\tif  j % 1000 == 0:\n","\t\t\tprint(\"Doing \" + data_type)\n","\t\telif j+1 % 100 == 0:\n","\t\t\tprint(\"\\t{}/1000\".format(j+1))\n","\t\tshuffled_df = shuffle(data_df).iloc[:sub_len[idx]]\n","\t\n","\t# Data preparation\n","\txtrain = shuffled_df[[\"Xl2fc\", \"Yl2fc\"]].values\n","\tytrain = shuffled_df[\"XYl2fc\"].values\n","\t\n","\t# Training\n","\tmodel = sklearn.linear_model.LinearRegression()\n","\tmodel.fit(xtrain, ytrain)\n","\n","\tcoefs = model.coef_\n","\tintercept = model.intercept_\n","\t\n","\tY_pred = model.predict(xtrain)\n","\tMAE = round(mean_absolute_error(ytrain, Y_pred), 2)\n","\tRMSE = round(mean_squared_error(ytrain, Y_pred, squared=False), 2)\n","\tr2 = round(r2_score(ytrain, Y_pred), 2)\n","\t\n","\tif permu == False:\n","\t\tminx, maxx = -12, 12\n","\t\tminy, maxy = -12, 12\n","\t\t\n","        # for mesh\n","\t\txs = np.tile(np.arange(minx, maxx), (len(np.arange(miny, maxy)),1))\n","\t\tys = np.tile(np.arange(miny, maxy), (len(np.arange(minx, maxx)),1)).T\n","\t\tzs = xs*coefs[0]+ys*coefs[1]+intercept\n","\t\t\n","\t\tplot_sub = shuffled_df[((shuffled_df.Xl2fc < 12) & (shuffled_df.Xl2fc > -12)) & ((shuffled_df.Yl2fc < 12) & (shuffled_df.Yl2fc > -12))]\n","\t\tplot_x = plot_sub[[\"Xl2fc\", \"Yl2fc\"]].values\n","\t\tplot_y = plot_sub[\"XYl2fc\"].values\n","\t\tplot_3d(plot_x, plot_y, title, MAE, RMSE, r2, xs, ys, zs, coefs, intercept)\n","\t\n","\treturn [coefs[0], coefs[1], intercept, MAE, RMSE, r2]\n","\n","#%% Actual plotting\n","ori_out = []\n","s_start = 0\n","for i, item in enumerate(stresses[s_start:], s_start):\n","\tprint(\"Doing \" + item)\n","\tori_out.append(process_data(item, i, False))\n","\n","ori_df = pd.DataFrame.from_dict({\"yval\" : [x[0] for x in ori_out] + [x[1] for x in ori_out] + [x[-1] for x in ori_out], \"Legend\": [\"Sx\"]*len(ori_out) + [\"Sy\"]*len(ori_out) + [\"R\" +  u\"\\u00b2\"]*len(ori_out), \"stress\": stresses_long + stresses_long + stresses_long})\n","#%%\n","colors = [\"royalblue\", \"crimson\", \"lightsteelblue\"]\n","col_pal = sns.color_palette(colors)\n","sns.set(font_scale=1.25)\n","sns.catplot(data=ori_df, kind=\"bar\",  x=\"stress\", y=\"yval\", hue=\"Legend\", aspect=1.7, palette=col_pal)\n","plt.ylabel(\"Magnitude\")\n","plt.xlabel(\"Sx\")\n","plt.savefig(dir_path + 'figures/' + \"Fig7H.png\", dpi=600)"],"metadata":{"id":"M3TjPPppOYpw"},"execution_count":null,"outputs":[]}]}