Upload data and script

README.md

 # RNAseq_ARA_TAZ
+RNA-seq Differential Expression, GSEA analysis and plots
 
+This repository contains R code for differential expression analysis using DESeq2 and Gene Set Enrichment Analysis (GSEA).
+
+The "data/" directory contains:
+
+- featureCounts_results_replica1.txt – Gene count matrix for replicate 1.
+- featureCounts_results_replica2.txt – Gene count matrix for replicate 2.
+- Decoder_replica1.txt – Sample annotation (metadata) for replicate 1.
+- Decoder_replica2.txt – Sample annotation (metadata) for replicate 2.
+- Villani_et_al_GeneSets_DCs.xlsx – Villani et al. dendritic cell gene signatures used for GSEA.
+- SenescenceSignatures.txt – Senescence-related gene signatures used for GSEA.
+
+# Citation:
+
+If you use this repository or data, please cite:
+
+Gilioli and Fusco et al., Therapy-Induced Senescence Promotes Immunogenicity in Acute Myeloid Leukemia Blasts through loss of PRC2-Mediated Gene Regulation
+
+Thank you.

data/Decoder_replica1.txt

+Sample	Group	Senescent	Patient
+C2	CTRL	Low	UPN02
+C3	CTRL	Low	UPN03
+C5	CTRL	Low	UPN05
+C10	CTRL	High	UPN10
+C11	CTRL	High	UPN11
+C13	CTRL	High	UPN13
+A2	TRT	Low	UPN02
+A3	TRT	Low	UPN03
+A5	TRT	Low	UPN05
+A10	TRT	High	UPN10
+A11	TRT	High	UPN11
+A13	TRT	High	UPN13
\ No newline at end of file

data/Decoder_replica2.txt

+Sample	Group	Senescent	Patient
+UPN05.ARA	ARA	Low	UPN05
+UPN05.TAZ	TAZ	Low	UPN05	
+UPN05.CTRL	CTRL	Low	UPN05
+X0826.ARA	ARA	Low	UPN03
+X0826.TAZ	TAZ	Low	UPN03
+X0826.CTRL	CTRL	Low	UPN03
+X0784.ARA	ARA	High	UPN22
+X0784.TAZ	TAZ	High	UPN22
+X0784.CTRL	CTRL	High	UPN22
+X0327.ARA	ARA	High	UPN10
+X0327.TAZ	TAZ	High	UPN10
+X0327.CTRL	CTRL	High	UPN10
+X0175.ARA	ARA	Low	UPN04
+X0175.TAZ	TAZ	Low	UPN04
+X0175.CTRL	CTRL	Low	UPN04

data/SenescenceSignatures.txt

+gs_name	gene_symbol
+Hernandez	ACADVL
+Hernandez	ADPGK
+Hernandez	ARHGAP35
+Hernandez	ARID2
+Hernandez	ASCC1
+Hernandez	B4GALT7
+Hernandez	BCL2L2
+Hernandez	C2CD5
+Hernandez	CCND1
+Hernandez	CHMP5
+Hernandez	CNTLN
+Hernandez	CREBBP
+Hernandez	DDA1
+Hernandez	DGKA
+Hernandez	DYNLT3
+Hernandez	EFNB3
+Hernandez	FAM214B
+Hernandez	GBE1
+Hernandez	GDNF
+Hernandez	GSTM4
+Hernandez	ICE1
+Hernandez	KCTD3
+Hernandez	KLC1
+Hernandez	MEIS1
+Hernandez	MT-CYB
+Hernandez	NFIA
+Hernandez	NOL3
+Hernandez	P4HA2
+Hernandez	PATZ1
+Hernandez	PCIF1
+Hernandez	PDLIM4
+Hernandez	PDS5B
+Hernandez	PLK3
+Hernandez	PLXNA3
+Hernandez	POFUT2
+Hernandez	RAI14
+Hernandez	RHNO1
+Hernandez	SCOC
+Hernandez	SLC10A3
+Hernandez	SLC16A3
+Hernandez	SMO
+Hernandez	SPATA6
+Hernandez	SPIN4
+Hernandez	STAG1
+Hernandez	SUSD6
+Hernandez	TAF13
+Hernandez	TMEM87B
+Hernandez	TOLLIP
+Hernandez	TRDMT1
+Hernandez	TSPAN13
+Hernandez	UFM1
+Hernandez	USP6NL
+Hernandez	ZBTB7A
+Hernandez	ZC3H4
+Hernandez	ZNHIT1
+FRIDMAN_SENESCENCE_UP	ALDH1A3
+FRIDMAN_SENESCENCE_UP	AOPEP
+FRIDMAN_SENESCENCE_UP	CCN2
+FRIDMAN_SENESCENCE_UP	CCND1
+FRIDMAN_SENESCENCE_UP	CD44
+FRIDMAN_SENESCENCE_UP	CDKN1A
+FRIDMAN_SENESCENCE_UP	CDKN1C
+FRIDMAN_SENESCENCE_UP	CDKN2A
+FRIDMAN_SENESCENCE_UP	CDKN2B
+FRIDMAN_SENESCENCE_UP	CDKN2D
+FRIDMAN_SENESCENCE_UP	CITED2
+FRIDMAN_SENESCENCE_UP	CLTB
+FRIDMAN_SENESCENCE_UP	COL1A2
+FRIDMAN_SENESCENCE_UP	CREG1
+FRIDMAN_SENESCENCE_UP	CRYAB
+FRIDMAN_SENESCENCE_UP	CXCL14
+FRIDMAN_SENESCENCE_UP	CYP1B1
+FRIDMAN_SENESCENCE_UP	EIF2S2
+FRIDMAN_SENESCENCE_UP	ESM1
+FRIDMAN_SENESCENCE_UP	F3
+FRIDMAN_SENESCENCE_UP	FILIP1L
+FRIDMAN_SENESCENCE_UP	FN1
+FRIDMAN_SENESCENCE_UP	GSN
+FRIDMAN_SENESCENCE_UP	GUK1
+FRIDMAN_SENESCENCE_UP	HBS1L
+FRIDMAN_SENESCENCE_UP	HPS5
+FRIDMAN_SENESCENCE_UP	HSPA2
+FRIDMAN_SENESCENCE_UP	HTATIP2
+FRIDMAN_SENESCENCE_UP	IFI16
+FRIDMAN_SENESCENCE_UP	IFNG
+FRIDMAN_SENESCENCE_UP	IGFBP1
+FRIDMAN_SENESCENCE_UP	IGFBP2
+FRIDMAN_SENESCENCE_UP	IGFBP3
+FRIDMAN_SENESCENCE_UP	IGFBP4
+FRIDMAN_SENESCENCE_UP	IGFBP5
+FRIDMAN_SENESCENCE_UP	IGFBP6
+FRIDMAN_SENESCENCE_UP	IGFBP7
+FRIDMAN_SENESCENCE_UP	IGSF3
+FRIDMAN_SENESCENCE_UP	ING1
+FRIDMAN_SENESCENCE_UP	IRF5
+FRIDMAN_SENESCENCE_UP	IRF7
+FRIDMAN_SENESCENCE_UP	ISG15
+FRIDMAN_SENESCENCE_UP	MAP1LC3B
+FRIDMAN_SENESCENCE_UP	MAP2K3
+FRIDMAN_SENESCENCE_UP	MDM2
+FRIDMAN_SENESCENCE_UP	MMP1
+FRIDMAN_SENESCENCE_UP	NDN
+FRIDMAN_SENESCENCE_UP	NME2
+FRIDMAN_SENESCENCE_UP	NRG1
+FRIDMAN_SENESCENCE_UP	OPTN
+FRIDMAN_SENESCENCE_UP	PEA15
+FRIDMAN_SENESCENCE_UP	RAB13
+FRIDMAN_SENESCENCE_UP	RAB31
+FRIDMAN_SENESCENCE_UP	RAB5B
+FRIDMAN_SENESCENCE_UP	RABGGTA
+FRIDMAN_SENESCENCE_UP	RAC1
+FRIDMAN_SENESCENCE_UP	RBL2
+FRIDMAN_SENESCENCE_UP	RGL2
+FRIDMAN_SENESCENCE_UP	RHOB
+FRIDMAN_SENESCENCE_UP	RRAS
+FRIDMAN_SENESCENCE_UP	S100A11
+FRIDMAN_SENESCENCE_UP	SERPINB2
+FRIDMAN_SENESCENCE_UP	SERPINE1
+FRIDMAN_SENESCENCE_UP	SMPD1
+FRIDMAN_SENESCENCE_UP	SMURF2
+FRIDMAN_SENESCENCE_UP	SOD1
+FRIDMAN_SENESCENCE_UP	SPARC
+FRIDMAN_SENESCENCE_UP	STAT1
+FRIDMAN_SENESCENCE_UP	TES
+FRIDMAN_SENESCENCE_UP	TFAP2A
+FRIDMAN_SENESCENCE_UP	TGFB1I1
+FRIDMAN_SENESCENCE_UP	THBS1
+FRIDMAN_SENESCENCE_UP	TNFAIP2
+FRIDMAN_SENESCENCE_UP	TNFAIP3
+FRIDMAN_SENESCENCE_UP	TP53
+FRIDMAN_SENESCENCE_UP	TSPYL5
+FRIDMAN_SENESCENCE_UP	VIM
+SASP_SCHLEICH	ANG
+SASP_SCHLEICH	AREG
+SASP_SCHLEICH	ATF4
+SASP_SCHLEICH	BHLHE40
+SASP_SCHLEICH	CCL1
+SASP_SCHLEICH	CCL11
+SASP_SCHLEICH	CCL12
+SASP_SCHLEICH	CCL2
+SASP_SCHLEICH	CCL20
+SASP_SCHLEICH	CCL25
+SASP_SCHLEICH	CD55
+SASP_SCHLEICH	CD55B
+SASP_SCHLEICH	CD9
+SASP_SCHLEICH	CPE
+SASP_SCHLEICH	CSF2
+SASP_SCHLEICH	CSF2RB
+SASP_SCHLEICH	CSF2RB2
+SASP_SCHLEICH	CSF3
+SASP_SCHLEICH	CTSB
+SASP_SCHLEICH	CXCL1
+SASP_SCHLEICH	CXCL12
+SASP_SCHLEICH	CXCL13
+SASP_SCHLEICH	CXCL2
+SASP_SCHLEICH	CXCL5
+SASP_SCHLEICH	CXCR2
+SASP_SCHLEICH	EGF
+SASP_SCHLEICH	EGFR
+SASP_SCHLEICH	EREG
+SASP_SCHLEICH	ETS1
+SASP_SCHLEICH	ETS2
+SASP_SCHLEICH	FAS
+SASP_SCHLEICH	FGF2
+SASP_SCHLEICH	FGF7
+SASP_SCHLEICH	FN1
+SASP_SCHLEICH	GEM
+SASP_SCHLEICH	GMFG
+SASP_SCHLEICH	HGF
+SASP_SCHLEICH	ICAM1
+SASP_SCHLEICH	ID1
+SASP_SCHLEICH	IFNG
+SASP_SCHLEICH	IGF1
+SASP_SCHLEICH	IGF2
+SASP_SCHLEICH	IGF2R
+SASP_SCHLEICH	IGFBP1
+SASP_SCHLEICH	IGFBP2
+SASP_SCHLEICH	IGFBP3
+SASP_SCHLEICH	IGFBP4
+SASP_SCHLEICH	IGFBP5
+SASP_SCHLEICH	IGFBP6
+SASP_SCHLEICH	IGFBP7
+SASP_SCHLEICH	IL13
+SASP_SCHLEICH	IL15
+SASP_SCHLEICH	IL1A
+SASP_SCHLEICH	IL1B
+SASP_SCHLEICH	IL6
+SASP_SCHLEICH	IL6ST
+SASP_SCHLEICH	IL7
+SASP_SCHLEICH	ITGA2
+SASP_SCHLEICH	ITPKA
+SASP_SCHLEICH	KITLG
+SASP_SCHLEICH	MIF
+SASP_SCHLEICH	MMP10
+SASP_SCHLEICH	MMP12
+SASP_SCHLEICH	MMP14
+SASP_SCHLEICH	MMP1A
+SASP_SCHLEICH	MMP3
+SASP_SCHLEICH	NGF
+SASP_SCHLEICH	PECAM1
+SASP_SCHLEICH	PIGF
+SASP_SCHLEICH	PLAT
+SASP_SCHLEICH	PLAU
+SASP_SCHLEICH	PLAUR
+SASP_SCHLEICH	PTGES
+SASP_SCHLEICH	SERPINB2
+SASP_SCHLEICH	SERPINE1
+SASP_SCHLEICH	TGFB1
+SASP_SCHLEICH	TIMP1
+SASP_SCHLEICH	TIMP2
+SASP_SCHLEICH	TNFRSF11B
+SASP_SCHLEICH	TNFRSF1A
+SASP_SCHLEICH	TNFRSF1B
+SASP_SCHLEICH	VEGFA
+CLASSICAL_SASP	BGN
+CLASSICAL_SASP	CCL2
+CLASSICAL_SASP	CCL20
+CLASSICAL_SASP	COL1A1
+CLASSICAL_SASP	CXCL12
+CLASSICAL_SASP	DCN
+CLASSICAL_SASP	EGF
+CLASSICAL_SASP	EFEMP1
+CLASSICAL_SASP	FGF2
+CLASSICAL_SASP	FGF7
+CLASSICAL_SASP	FN1
+CLASSICAL_SASP	CSF1
+CLASSICAL_SASP	CTGF
+CLASSICAL_SASP	CXCL1
+CLASSICAL_SASP	GDNF
+CLASSICAL_SASP	IGFBP2
+CLASSICAL_SASP	IGFBP6
+CLASSICAL_SASP	IFNA1
+CLASSICAL_SASP	IFNB1
+CLASSICAL_SASP	IFNG
+CLASSICAL_SASP	IL1A
+CLASSICAL_SASP	IL1B
+CLASSICAL_SASP	IL6
+CLASSICAL_SASP	IL7
+CLASSICAL_SASP	CXCL8
+CLASSICAL_SASP	IL13
+CLASSICAL_SASP	IL15
+CLASSICAL_SASP	KITLG
+CLASSICAL_SASP	CCL2
+CLASSICAL_SASP	CCL9
+CLASSICAL_SASP	MMP2
+CLASSICAL_SASP	MMP3
+CLASSICAL_SASP	MMP9
+CLASSICAL_SASP	SERPINE1
+CLASSICAL_SASP	TGFB1
+CLASSICAL_SASP	THBS1
+CLASSICAL_SASP	TIMP2
+CLASSICAL_SASP	TNF
+CLASSICAL_SASP	VEGFA

script/RNAseq_Analysis.R

+#Import libraries
+library(DESeq2)
+library(ggplot2)
+library(ggpubr)
+library(reshape2)
+library(ggrepel)
+library(ggh4x)
+library(clusterProfiler)
+library(msigdbr)
+library(readxl)
+library(enrichplot)
+library(enrichR)
+
+#Set Working directory
+setwd(".")
+
+#PRC2 Subunits
+PRC2 = c("ASXL1", "EZH2", "EED", "SUZ12", "RBBP4", "RBBP7")
+
+#Read Villani et al. Signature
+villani_genes = read.xlsx("data/Villani_et_al_GeneSets_DCs.xlsx")
+
+#Read Senescence signatures
+senescence_signatures = read.table("data/SenescenceSignatures.txt")
+
+### Define Functions
+
+#Barplot Gene expression lfc
+barplot_gene_expression_lfc = function(res, genes, title){
+  p = ggplot(res[rownames(res) %in% genes, ], aes(x = Gene, y = log2FoldChange)) +
+    geom_bar(stat = "identity", color = "black") +
+    theme_bw(base_size = 15) +
+    xlab("")+
+    ylab("log2(FoldChange)") +
+    ggtitle(title) +
+    theme(text=element_text(family="Arial", face="plain", size=15))
+  return(p)
+}
+
+
+sort_by_lfc = function(frame, sorttype = "standard"){
+  frame = frame[!is.na(frame$padj),]
+  if (sorttype == "standard"){
+    original_gene_list = frame$log2FoldChange
+  }
+  if (sorttype == "composite"){
+    original_gene_list = sign(frame$log2FoldChange) * -log10(frame$pvalue)
+  }
+  names(original_gene_list) = rownames(frame)
+  gene_list <- na.omit(original_gene_list)
+  gene_list = gene_list[!is.infinite(gene_list)]
+  gene_list = sort(gene_list, decreasing = TRUE)
+  return(gene_list)
+  
+}
+
+#GSEA Analysis
+run_GSEA = function(frame.to.test, signatures.to.test, padj = 0.25, nPermSimple = 1000, 
+                    sorttype = "standard", width = 10, height = 10){
+  #Iterate over all comparisons
+  for (name in names(frame.to.test)){ 
+    
+    input = sort_by_lfc(frame.to.test[[name]], sorttype = sorttype) 
+    
+    #Iterate over all signatures
+    for (signature in names(signatures.to.test)){
+      print(paste("Testing...", name, "on", signature, sep = " " ))
+      
+      gse = GSEA(geneList = input,
+                 TERM2GENE = signatures.to.test[[signature]],
+                 by = "fgsea",
+                 maxGSSize= 10000000,
+                 minGSSize = 1,
+                 nPermSimple = nPermSimple, 
+                 pvalueCutoff = padj)
+      
+      gse.res = gse@result
+      gse.res$ID = factor(gse.res$ID, levels = gse.res$ID[order(gse.res$NES)])
+      gse.res = gse.res[order(gse.res$NES),]
+      
+      
+      p = ggplot(gse.res, aes(y = NES, x = ID, fill = p.adjust )) +
+        geom_bar(stat = "identity") +
+        coord_flip() +
+        theme_bw() +
+        scale_fill_gradient(high = "skyblue", low = "tomato") +
+        ggtitle(paste(signature, name, sep = " | "))
+      
+      p1 = gseaplot2(gse, geneSetID = 1:nrow(gse), pvalue_table = TRUE, title = paste(signature, name, sep = " | "), base_size = 20)
+      print(p)
+      print(p1)
+      
+    }
+    
+  }
+} 
+
+###
+
+
+### Replica 1
+
+#Read gene and decoder tables
+gene = read.table(file = "data/featureCounts_results_replica1.txt")
+decoder = read.delim(file = "data/Decoder_replica1.txt", header = T)
+rownames(decoder) = decoder$Sample
+  
+#Set same order in gene and decoder
+gene = gene[, rownames(decoder)]
+  
+#Add variables to decoder file
+decoder$Senescent = relevel(factor(decoder$Senescent), ref = "Low")
+decoder$Group = relevel(factor(decoder$Group), ref = "CTRL")
+decoder$Group2 = paste(decoder$Senescent, decoder$Group, sep = "_")
+decoder$Group2 = factor(decoder$Group2, levels = c("Low_CTRL", "Low_TRT", "High_CTRL", "High_TRT"))
+
+#Build DESeq2 object
+dds <- DESeqDataSetFromMatrix(countData = gene,
+                              colData = decoder,
+                              design= ~ Senescent + Group + Senescent:Group )
+#Run DESeq2
+dds <- DESeq(dds)
+
+dds1 <- DESeqDataSetFromMatrix(countData = gene,
+                               colData = decoder,
+                               design= ~ Group2 )
+#Run DESeq2
+dds1 <- DESeq(dds1)
+
+#Extract comparisons
+res.HighTRT.LowTRT = as.data.frame(DESeq2::results(dds, name =  "SenescentHigh.GroupTRT"))
+res.HighTRT.LowTRT$Gene = rownames(res.HighTRT.LowTRT)
+res.LowTRT.HighTRT = as.data.frame(DESeq2::results(dds, contrast =  c(0,0,0,-1)))
+res.LowTRT.HighTRT$Gene = rownames(res.LowTRT.HighTRT)
+res.TRT.CTRL = as.data.frame(DESeq2::results(dds, contrast = c("Group", "TRT", "CTRL")))
+res.TRT.CTRL$Gene = rownames(res.TRT.CTRL)
+res.High.Low = as.data.frame(DESeq2::results(dds, contrast = c("Senescent", "High", "Low")))
+res.High.Low$Gene = rownames(res.High.Low)
+res.LowTRT.LowCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "Low_TRT", "Low_CTRL")))
+res.LowTRT.LowCTRL$Gene = rownames(res.LowTRT.LowCTRL)
+res.HighTRT.HighCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "High_TRT", "High_CTRL")))
+res.HighTRT.HighCTRL$Gene = rownames(res.HighTRT.HighCTRL)
+res.HighCTRL.LowCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "High_CTRL", "Low_CTRL")))
+res.HighCTRL.LowCTRL$Gene = rownames(res.HighCTRL.LowCTRL)
+
+
+### Run GSEA analysis
+
+#Create list of dataframes
+frame.to.test = list("HighTRT.LowTRT" = res.HighTRT.LowTRT,
+                     "LowTRT.HighTRT" = res.LowTRT.HighTRT,
+                     "High.Low" = res.High.Low, 
+                     "TRT.CTRL" = res.TRT.CTRL, 
+                     "LowTRT.LowCTRL" = res.LowTRT.LowCTRL,
+                     "HighTRT.HighCTRL" = res.HighTRT.HighCTRL, 
+                     "HighCTRL.LowCTRL" = res.HighCTRL.LowCTRL)
+
+#Create list of signatures
+signatures.to.test = list("Villani" = villani_genes, 
+                          "Senescence" = senescence_signatures)
+
+#GSEA (GSEA Plot + GSEA Barplot)
+run_GSEA(frame.to.test, signatures.to.test)
+  
+###
+
+###
+
+
+  
+### Replica 2
+#Read  gene x counts table
+gene = read.table(file = "data/featureCounts_results_replica2.txt", header = T)
+gene = gene[,-c(2,3,4,5,6)]
+rownames(gene) = gene$Geneid
+gene$Geneid = NULL
+colnames(gene) = gsub("BAM.", "", colnames(gene))
+colnames(gene) = unlist(lapply(colnames(gene), function(x) Reduce(function(x,y) paste(x,y, sep = "."), 
+                                                                  unlist(strsplit(x, "\\."))[c(1,2)])
+))
+
+
+#Read decoder file
+decoder = read.table(file = "data/Decoder_replica2.txt", header = T)
+decoder$Sample = gsub("X", "", decoder$Sample)
+rownames(decoder) = decoder$Sample
+
+#Set same order in gene and decoder
+gene = gene[, rownames(decoder)]
+
+#Add variables to decoder file
+decoder$Senescent = relevel(factor(decoder$Senescent), ref = "Low")
+decoder$Group = relevel(factor(decoder$Group), ref = "CTRL")
+decoder$Group2 = paste(decoder$Senescent, decoder$Group, sep = "_")
+decoder$Group2 = factor(decoder$Group2, levels =  c("Low_CTRL", "Low_ARA", "Low_TAZ", "High_CTRL", "High_ARA", "High_TAZ"))
+
+
+#Build DESeq2 object
+dds <- DESeqDataSetFromMatrix(countData = gene,
+                              colData = decoder,
+                              design= ~ Senescent + Group + Senescent:Group )
+#Run DESeq2
+dds <- DESeq(dds)
+  
+#Build DESeq2 object
+dds1 <- DESeqDataSetFromMatrix(countData = gene,
+                               colData = decoder,
+                               design= ~ Group2 )
+#Run DESeq2
+dds1 <- DESeq(dds1)
+  
+#Extract results
+res.HighARA.LowARA = as.data.frame(DESeq2::results(dds, name =  "SenescentHigh.GroupARA"))
+res.HighARA.LowARA$Gene = rownames(res.HighARA.LowARA)
+res.HighTAZ.LowTAZ = as.data.frame(DESeq2::results(dds, name =  "SenescentHigh.GroupTAZ"))
+res.HighTAZ.LowTAZ$Gene = rownames(res.HighTAZ.LowTAZ)
+res.LowARA.HighARA = as.data.frame(DESeq2::results(dds, contrast = c(0,0,0,0,-1,0)))
+res.LowARA.HighARA$Gene = rownames(res.LowARA.HighARA)
+res.LowTAZ.HighTAZ = as.data.frame(DESeq2::results(dds, contrast = c(0,0,0,0,0,-1)))
+res.LowTAZ.HighTAZ$Gene = rownames(res.LowTAZ.HighTAZ)
+res.ARA.CTRL = as.data.frame(DESeq2::results(dds, contrast = c("Group", "ARA", "CTRL")))
+res.ARA.CTRL$Gene = rownames(res.ARA.CTRL)
+res.TAZ.CTRL = as.data.frame(DESeq2::results(dds, contrast = c("Group", "TAZ", "CTRL")))
+res.TAZ.CTRL$Gene = rownames(res.TAZ.CTRL)
+res.High.Low = as.data.frame(DESeq2::results(dds, contrast = c("Senescent", "High", "Low")))
+res.High.Low$Gene = rownames(res.High.Low)
+res.LowARA.LowCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "Low_ARA", "Low_CTRL")))
+res.LowARA.LowCTRL$Gene = rownames(res.LowARA.LowCTRL)
+res.HighARA.HighCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "High_ARA", "High_CTRL")))
+res.HighARA.HighCTRL$Gene = rownames(res.HighARA.HighCTRL)
+res.LowTAZ.LowCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "Low_TAZ", "Low_CTRL")))
+res.LowTAZ.LowCTRL$Gene = rownames(res.LowTAZ.LowCTRL)
+res.HighTAZ.HighCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "High_TAZ", "High_CTRL")))
+res.HighTAZ.HighCTRL$Gene = rownames(res.HighTAZ.HighCTRL)
+res.HighCTRL.LowCTRL = as.data.frame(DESeq2::results(dds1, contrast = c("Group2", "High_CTRL", "Low_CTRL")))
+res.HighCTRL.LowCTRL$Gene = rownames(res.HighCTRL.LowCTRL)
+
+
+#Do Barplot of PRC2 Subunits (Log2(FoldChange))
+tmp = res.HighCTRL.LowCTRL
+#Low vs High, invert FoldChange
+tmp$log2FoldChange = -tmp$log2FoldChange
+#Barplot
+barplot_gene_expression_lfc(tmp, PRC2, "Low Untreated vs High Untreated")
+
+  
+### Run GSEA  
+
+#Create list of dataframes
+frame.to.test = list("HighARA.LowARA" = res.HighARA.LowARA,
+                     "HighTAZ.LowTAZ" = res.HighTAZ.LowTAZ,
+                     "LowARA.HighARA" = res.LowARA.HighARA,
+                     "LowTAZ.HighTAZ" = res.LowTAZ.HighTAZ,
+                     "High.Low" = res.High.Low, 
+                     "ARA.CTRL" = res.ARA.CTRL,
+                     "TAZ.CTRL" = res.TAZ.CTRL, 
+                     "LowARA.LowCTRL" = res.LowARA.LowCTRL,
+                     "LowTAZ.LowCTRL" = res.LowTAZ.LowCTRL,
+                     "HighARA.HighCTRL" = res.HighARA.HighCTRL,
+                     "HighTAZ.HighCTRL" = res.HighTAZ.HighCTRL,
+                     "HighCTRL.LowCTRL" = res.HighCTRL.LowCTRL)
+
+#Create list of signatures
+signatures.to.test = list("Villani" = villani_genes, 
+                          "Senescence" = senescence_signatures)
+  
+  
+#Run GSEA analysis
+run_GSEA(frame.to.test, signatures.to.test, padj = 1, nPermSimple = 10000,
+         sorttype = "composite", width = 10, height = 8)
+  
+###
+
+###
+  
\ No newline at end of file