############################# # Aspen Coyle, afcoyle@uw.edu # Roberts lab, UW-SAFS # 2021-02-01 ############################# # This file runs GO-MWU. Commands copied from the GO-MWU.R file # in the GO-MWU Github repository, available at https://github.com/z0on/GO_MWU # We will run GO-MWU twice - once for each comparison library(ape) library(tidyverse) # Need to be in same directory as all other GO-MWU files - # both data files and analysis files setwd("1_6_running_GO-MWU") #### GO-MWU Run 1: Elevated Day 2 vs. Ambient Days 0+2, Individual Libraries Only -------------------------- # Edit these to match your data file names: input="cbai2.0_elev2_vs_amb02_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_elev2_vs_amb02_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 94 GO terms at 10% FDR # We'll now manually save this file in the graphs/GOMWU_results folder grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.1, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) #plot(results[[2]],cex=0.6) #abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 2: Elevated Day 2 vs. Ambient Day 0+2+17 + Elevated Day 0 + Lowered Day 0 -------------------------- # Edit these to match your data file names: input="cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 52 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 3: Elevated Day 0 vs. Elevated Day 2, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_elev0_vs_elev2_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_elev0_vs_elev2_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 74 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 4: Ambient Day 0 vs. Ambient Day 2, Individual Libraries Only -------------------------- # Edit these to match your data file names: input="cbai2.0_amb0_vs_amb2_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_amb0_vs_amb2_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 16 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 5: Ambient Day 0 vs. Ambient Day 17, Individual Libraries Only -------------------------- # Edit these to match your data file names: input="cbai2.0_amb0_vs_amb17_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_amb0_vs_amb17_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 152 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.01, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.005, # FDR cutoff to print in regular (not italic) font. level3=0.0001, # FDR cutoff to print in large bold font. txtsize=0.9, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 6: Ambient Day 2 vs. Ambient Day 17, Individual Libraries Only -------------------------- # Edit these to match your data file names: input="cbai2.0_amb2_vs_amb17_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_amb2_vs_amb17_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 177 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.01, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.005, # FDR cutoff to print in regular (not italic) font. level3=0.001, # FDR cutoff to print in large bold font. txtsize=1.0, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 7: Ambient Day 2 vs. Elevated Day 2, Individual Libraries Only -------------------------- # Edit these to match your data file names: input="cbai2.0_amb2_vs_elev2_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_amb2_vs_elev2_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 71 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 8: Decreased Day 0 vs. Decreased Day 2, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_vs_low2_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_vs_low2_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 163 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.01, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.005, # FDR cutoff to print in regular (not italic) font. level3=0.001, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 9: Decreased Day 0 vs. Decreased Day 2 + 17, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_vs_low217_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_vs_low217_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 60 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 10: Decreased Day 0 + Ambient Day 0 + 2 + 17 vs. Decreased Day 2 + 17, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_amb0217_vs_low217_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_amb0217_vs_low217_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 58 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 11: Decreased Day 2 + 17 vs. Ambient Day 2 + 17, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low217_vs_amb217_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low217_vs_amb217_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 75 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 9: Decreased Day 0 vs. Ambient Day 0, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_vs_amb0_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_vs_amb0_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 122 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.01, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.005, # FDR cutoff to print in regular (not italic) font. level3=0.001, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 10: Decreased Day 0 vs. Decreased Day 17, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_vs_low17_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_vs_low17_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 0 GO terms at 10% FDR # Ending analysis here # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 11: Decreased Day 0 + Ambient Day 0 vs. Decreased Day 17 + Ambient Day 17, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_amb0_vs_low17_amb17_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_amb0_vs_low17_amb17_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 94 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.1, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.05, # FDR cutoff to print in regular (not italic) font. level3=0.01, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_")) #### GO-MWU Run 12: Decreased Day 0 + Elevated Day 0 vs. Decreased Day 2 + Elevated Day 2, Indiv. Libraries Only --------------------------------------- # Edit these to match your data file names: input="cbai2.0_low0_elev0_vs_low2_elev2_indiv_l2FC.csv" # two columns of comma-separated values: gene id, continuous measure of significance. To perform standard GO enrichment analysis based on Fisher's exact test, use binary measure (0 or 1, i.e., either sgnificant or not). goAnnotations="cbai2.0_low0_elev0_vs_low2_elev2_indiv_GOIDs_norepeats.txt" # two-column, tab-delimited, one line per gene, multiple GO terms separated by semicolon. If you have multiple lines per gene, use nrify_GOtable.pl prior to running this script. goDatabase="go.obo" # download from http://www.geneontology.org/GO.downloads.ontology.shtml goDivision="BP" # either MF, or BP, or CC source("gomwu.functions.R") # ------------- Calculating stats # It might take a few minutes for MF and BP. Do not rerun it if you just want to replot the data with different cutoffs, go straight to gomwuPlot. If you change any of the numeric values below, delete the files that were generated in previos runs first. gomwuStats(input, goDatabase, goAnnotations, goDivision, perlPath="C:/Users/acoyl/Documents/GradSchool/RobertsLab/Tools/perl/bin/perl.exe", # replace with full path to perl executable if it is not in your system's PATH already largest=0.1, # a GO category will not be considered if it contains more than this fraction of the total number of genes smallest=5, # a GO category should contain at least this many genes to be considered clusterCutHeight=0.25, # threshold for merging similar (gene-sharing) terms. See README for details. # Alternative="g" # by default the MWU test is two-tailed; specify "g" or "l" of you want to test for "greater" or "less" instead. # Module=TRUE,Alternative="g" # un-remark this if you are analyzing a SIGNED WGCNA module (values: 0 for not in module genes, kME for in-module genes). In the call to gomwuPlot below, specify absValue=0.001 (count number of "good genes" that fall into the module) # Module=TRUE # un-remark this if you are analyzing an UNSIGNED WGCNA module ) # --------------- Results # 142 GO terms at 10% FDR grDevices::windows() results=gomwuPlot(input,goAnnotations,goDivision, #absValue=0.05, # genes with the measure value exceeding this will be counted as "good genes". This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes"). absValue=1, # un-remark this if you are using log2-fold changes level1=0.01, # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue. level2=0.005, # FDR cutoff to print in regular (not italic) font. level3=0.001, # FDR cutoff to print in large bold font. txtsize=1.2, # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text) for better "word cloud" effect treeHeight=0.5, # height of the hierarchical clustering tree # colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") # these are default colors, un-remar and change if needed ) # manually rescale the plot so the tree matches the text # if there are too many categories displayed, try make it more stringent with level1=0.05,level2=0.01,level3=0.001. # text representation of results, with actual adjusted p-values results[[1]] # ------- extracting representative GOs # this module chooses GO terms that best represent *independent* groups of significant GO terms pcut=1e-2 # adjusted pvalue cutoff for representative GO hcut=0.9 # height at which cut the GO terms tree to get "independent groups". # plotting the GO tree with the cut level (un-remark the next two lines to plot) # plot(results[[2]],cex=0.6) # abline(h=hcut,col="red") # cutting ct=cutree(results[[2]],h=hcut) annots=c();ci=1 for (ci in unique(ct)) { message(ci) rn=names(ct)[ct==ci] obs=grep("obsolete",rn) if(length(obs)>0) { rn=rn[-obs] } if (length(rn)==0) {next} rr=results[[1]][rn,] bestrr=rr[which(rr$pval==min(rr$pval)),] best=1 if(nrow(bestrr)>1) { nns=sub(" .+","",row.names(bestrr)) fr=c() for (i in 1:length(nns)) { fr=c(fr,eval(parse(text=nns[i]))) } best=which(fr==max(fr)) } if (bestrr$pval[best]<=pcut) { annots=c(annots,sub("\\d+\\/\\d+ ","",row.names(bestrr)[best]))} } mwus=read.table(paste("MWU",goDivision,input,sep="_"),header=T) bestGOs=mwus[mwus$name %in% annots,] bestGOs # Move the 3 files we created to a permanent folder, since GO-MWU automatically puts them in # the same folder you run the script in file_loc <- input %>% str_remove("cbai2.0_") %>% str_remove("_l2FC.csv") filepath <- paste0("../../output/GO-MWU_output/cbai_transcriptomev2.0/", file_loc, "/") # We're encountering some issues with a double-named filename (something like "dissim_BP_cbai2.0_amb0217_elev0_low0_vs_elev2_l2FC.csv_cbai2.0_amb0217_elev0_low0_vs_elev2_GOIDs_norepeats.txt") # It's a problem with the function, but I don't want to touch the prebuilt GO-MWU stuff, # so I'm just going to remove part of the name files <- list.files(getwd(), pattern = "BP_") newfiles <- gsub("\\.csv_cbai.*", ".txt", files) file.rename(files, newfiles) file.copy(list.files(getwd(), pattern = "BP_"), filepath) file.remove(list.files(getwd(), pattern = "BP_"))