Notebook RCM Part 1 ccRCC Figure 4
Plots for figure 4 Early vs Late¶
- Correlation between RNA and DNA methylation
- Correlation between RNA and protein
- ORA overlap between protein DNA methylation and RNA
We also want to calculate the DNA methylation Beta value as opposed to just the M value change since we want to mainly use this!
In [20]:
# Need to add in entrez gene ID and also make labels for genes for ORA
# Imports
from scircm import * # Note if you have a mac M1 use from sircle import * and you won't be able to do 7,8
import seaborn as sns
from sciutil import *
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
import os
# Setup file locations and label of the cancer
u = SciUtil()
cancer = 'ClearCellRenalCellCarcinoma'
input_dir = 'Input_RCM'
output_dir = 'Output_Data'
supp_dir = 'Required_Refs'
fig_dir = 'Output_Figures'
regLabel = 'RG2_Changes_filtered'
# Use same files as in RCM
DNA_methylation_file = os.path.join('Input_Methylation', f'{cancer}_filtered_DCpG_EarlyLate.csv')
RNA_file = os.path.join('Input_RNAseq', f'{cancer}_filtered_DE_RNA_EarlyLate.csv')
protein_file = os.path.join('Input_Protein', f'{cancer}_filtered_DA_Protein_EarlyLate.csv')
# Names of columns in each of the above files
logFC_RNA_column = "logFC_rna"
padj_RNA_column = "padj_rna"
betadiff_DNA_methylation_column = "beta_diff"
padj_DNA_methylation_column = 'adj.P.Val'
logFC_protein_column = "logFC_protein"
padj_protein_column = "padj_protein"
# NOTE: all of the above files MUST have this column i.e. they must all have the samely named gene ID column
gene_id_column = "gene_name"
# Use same cutoffs as in RCM
rna_padj_cutoff = 0.05
prot_padj_cutoff = 0.05
meth_padj_cutoff = 0.05
rna_logfc_cutoff = 0.25
prot_logfc_cutoff = 0.1
meth_diff_cutoff = 0.05
# Join the DCpG file with the epic manifest and then filter the file
# Annotate the gene names to entrez gene IDs using annotation file from HG38
annot = pd.read_csv(os.path.join(supp_dir, 'hsapiens_gene_ensembl-GRCh38.p13.csv'))
annot = annot.dropna(subset=['external_gene_name', 'entrezgene_id'])
annot = annot.drop_duplicates(subset='external_gene_name')
name_to_entrez = dict(zip(annot.external_gene_name, annot.entrezgene_id))
/var/folders/gq/6ljhmvm1713fykdjqbl188pm0000gn/T/ipykernel_56074/33642123.py:49: DtypeWarning: Columns (2) have mixed types. Specify dtype option on import or set low_memory=False. annot = pd.read_csv(os.path.join(supp_dir, 'hsapiens_gene_ensembl-GRCh38.p13.csv'))
Figure 1A¶
Plot and calculate the DNA methylation and protein correlation
In [21]:
# RNA and Methylation correlation
from sciviso import *
from scipy.stats import pearsonr, spearmanr
from sciutil import *
import matplotlib.pyplot as plt
u = SciUtil()
m_df = pd.read_csv(DNA_methylation_file)
df = pd.read_csv(RNA_file)
m_df.set_index(gene_id_column, inplace=True)
df.set_index(gene_id_column, inplace=True)
# Join on gene name
joined_df = m_df.join(df, how='inner')
# Correlation
u.dp(['Pearsons R:', pearsonr(joined_df[betadiff_DNA_methylation_column].values, joined_df[logFC_RNA_column].values)])
u.dp(['Spearmans rho:', spearmanr(joined_df[betadiff_DNA_methylation_column].values, joined_df[logFC_RNA_column].values)])
# plot and print
sc = Scatterplot(joined_df, betadiff_DNA_methylation_column, logFC_RNA_column,
'DNA & RNA', 'DNA M. diff', 'RNA logFC',
config={'figsize': (1.5, 1.5), 'opacity': 0.2},
colour='black')
sc.plot()
plt.savefig(os.path.join(fig_dir, f'Figure4A_MethylationRNA_EarlyLate.svg'))
plt.show()
# Do the same with RNA and protein
p_df = pd.read_csv(protein_file)
df = pd.read_csv(RNA_file)
p_df.set_index(gene_id_column, inplace=True)
df.set_index(gene_id_column, inplace=True)
# Join on gene name
joined_df = p_df.join(df, how='inner')
# Correlation
u.dp(['Pearsons R:', pearsonr(joined_df[logFC_protein_column].values, joined_df[logFC_RNA_column].values)])
u.dp(['Spearmans rho:', spearmanr(joined_df[logFC_protein_column].values, joined_df[logFC_RNA_column].values)])
# plot and print
sc = Scatterplot(joined_df, logFC_protein_column, logFC_RNA_column, 'Protein & RNA',
'Protein logFC', 'RNA logFC',
config={'figsize': (1.5, 1.5), 'opacity': 0.2},
colour='black')
sc.plot()
plt.savefig(os.path.join(fig_dir, f'Figure4A_ProteinRNA_EarlyLate.svg'))
-------------------------------------------------------------------------------- Pearsons R: PearsonRResult(statistic=-0.08615947225078259, pvalue=0.038379169773448005) -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- Spearmans rho: SignificanceResult(statistic=-0.10924771284138445, pvalue=0.00857171446474042) --------------------------------------------------------------------------------
/Users/ariane/opt/miniconda3/envs/roundround/lib/python3.10/site-packages/sciviso/scatterplot.py:114: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored scatter = ax.scatter(vis_df[x].values, vis_df[y].values, c=self.colour, alpha=self.opacity,
-------------------------------------------------------------------------------- Pearsons R: PearsonRResult(statistic=0.3077873177806897, pvalue=8.364488521279044e-186) -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- Spearmans rho: SignificanceResult(statistic=0.375517201841577, pvalue=1.1726957157934435e-282) --------------------------------------------------------------------------------
/Users/ariane/opt/miniconda3/envs/roundround/lib/python3.10/site-packages/sciviso/scatterplot.py:114: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored scatter = ax.scatter(vis_df[x].values, vis_df[y].values, c=self.colour, alpha=self.opacity,
Figure 1B¶
Plot and visualise the overlap of ORA for the top terms, first we need to make "groups" for the different ORAs and top genes.
We do this for each of the above files (i.e. DNA methylation, gene expression and protein)
In [23]:
# -------------------------------------------------------------------
# Make ORA groups for DNA methylation
# -------------------------------------------------------------------
# Add in entrez gene ID and a label for each gene for DNA methylation
m_df = pd.read_csv(DNA_methylation_file)
entrez, labels = [], []
padjs, logfcs = m_df[padj_DNA_methylation_column].values, m_df[betadiff_DNA_methylation_column].values
for i, gene in enumerate(m_df[gene_id_column].values):
label = 'NA'
entrez.append(name_to_entrez.get(gene))
if padjs[i] < meth_padj_cutoff:
if logfcs[i] < (-1 * meth_diff_cutoff):
label = 'Hypo' # Hypo methylation
elif logfcs[i] > meth_diff_cutoff:
label = 'Hyper' # Hyper methylation i.e. "maybe" driving down reg.
labels.append(label)
m_df['ORA_label'] = labels
m_df['entrezgene_id'] = entrez
m_df = m_df[[gene_id_column, 'entrezgene_id', 'ORA_label', padj_DNA_methylation_column,
betadiff_DNA_methylation_column]]
m_df = m_df.dropna()
m_df.to_csv(os.path.join(output_dir, f'DNAMethylation_{cancer}_ORA_EarlyLate.csv'), index=False)
# -------------------------------------------------------------------
# Make ORA groups for Gene expression
# -------------------------------------------------------------------
df = pd.read_csv(RNA_file)
entrez, labels = [], []
padjs, logfcs = df[padj_RNA_column].values, df[logFC_RNA_column].values
for i, gene in enumerate(df[gene_id_column].values):
label = 'NA'
entrez.append(name_to_entrez.get(gene))
if padjs[i] < rna_padj_cutoff:
if logfcs[i] < (-1 * rna_logfc_cutoff):
label = 'DOWN' # Repressed
elif logfcs[i] > rna_logfc_cutoff:
label = 'UP' # Increased Expression
labels.append(label)
df['ORA_label'] = labels
df['entrezgene_id'] = entrez
df = df[[gene_id_column, 'entrezgene_id', 'ORA_label', padj_RNA_column, logFC_RNA_column]]
df = df.dropna()
df.to_csv(os.path.join(output_dir, f'RNA_{cancer}_ORA_EarlyLate.csv'), index=False)
# -------------------------------------------------------------------
# Make ORA groups for Protein
# -------------------------------------------------------------------
df = pd.read_csv(protein_file)
entrez, labels = [], []
padjs, logfcs = df[padj_protein_column].values, df[logFC_protein_column].values
for i, gene in enumerate(df[gene_id_column].values):
label = 'NA'
entrez.append(name_to_entrez.get(gene))
if padjs[i] < prot_padj_cutoff:
if logfcs[i] < (-1 * prot_logfc_cutoff):
label = 'DOWN' # Repressed
elif logfcs[i] > prot_logfc_cutoff:
label = 'UP' # Increased expression
labels.append(label)
df['ORA_label'] = labels
df['entrezgene_id'] = entrez
df = df[[gene_id_column, 'entrezgene_id', 'ORA_label', padj_protein_column, logFC_protein_column]]
df = df.dropna()
df.to_csv(os.path.join(output_dir, f'Protein_{cancer}_ORA_EarlyLate.csv'), index=False)
In [24]:
from matplotlib_venn import venn3
ora_dir = os.path.join(output_dir, 'ORA_EarlyLate')
# Get the overlapping GO terms that were significant
p_go = pd.read_csv(os.path.join(output_dir, f'Protein_{cancer}_ORA_EarlyLate.csv'))
p_go = p_go[p_go['ORA_label'] == 'UP']['gene_name'].values # Get sig terms
m_go = pd.read_csv(os.path.join(output_dir, f'DNAMethylation_{cancer}_ORA_EarlyLate.csv'))
m_go = m_go[m_go['ORA_label'] == 'Hypo']['gene_name'].values # Get sig terms
r_go = pd.read_csv(os.path.join(output_dir, f'RNA_{cancer}_ORA_EarlyLate.csv'))
r_go = r_go[r_go['ORA_label'] == 'UP']['gene_name'].values # Get sig terms
venn3([set(p_go), set(r_go), set(m_go)], set_labels=['Proteomics', 'RNAseq', 'DNA Methylation'],
set_colors=('#0101d7', '#e68e25', '#6aaf44'),)
plt.rcParams['svg.fonttype'] = 'none'
plt.title('Up')
plt.savefig(os.path.join(fig_dir, f'Figure4B_UpGenes_EarlyLate.svg'))
plt.show()
# Get the overlapping GO terms that were significant
p_go = pd.read_csv(os.path.join(output_dir, f'Protein_{cancer}_ORA_EarlyLate.csv'))
p_go = p_go[p_go['ORA_label'] == 'DOWN']['gene_name'].values # Get sig terms
m_go = pd.read_csv(os.path.join(output_dir, f'DNAMethylation_{cancer}_ORA_EarlyLate.csv'))
m_go = m_go[m_go['ORA_label'] == 'Hyper']['gene_name'].values # Get sig terms
r_go = pd.read_csv(os.path.join(output_dir, f'RNA_{cancer}_ORA_EarlyLate.csv'))
r_go = r_go[r_go['ORA_label'] == 'DOWN']['gene_name'].values # Get sig terms
venn3([set(p_go), set(r_go), set(m_go)], set_labels=['Proteomics', 'RNAseq', 'DNA Methylation'],
set_colors=('#0101d7', '#e68e25', '#6aaf44'),)
plt.rcParams['svg.fonttype'] = 'none'
plt.title('Down')
plt.savefig(os.path.join(fig_dir, 'Figure4B_DownGenes_EarlyLate.svg'))
Look at overlap now (after going to R and running SiRCle ORA)¶
Code to run ORA on each of the above files:
{r}
library(EnhancedVolcano)
library(org.Hs.eg.db)#For Human load:
library(clusterProfiler)#To run ORA
library(enrichplot)#For emapplot, dotplot,...
library(ggplot2)#For saving and legends
library(tidyverse) # used for data manipulation
#Establish a couple of functions needed to perform the GSEA (made by Aurelien)
library(fgsea)
library(GSEABase)
# ----------------------------------------------------------------------------------------------
# File and naming setup
# ----------------------------------------------------------------------------------------------
GSEA_file_dir <- 'Required_Refs/GSEA/' # Will need to change the path if you're on windows
figure_dir <- 'Output_Figures'
input_data_dir <- 'Output_Data'
output_data_dir <- 'Output_Data/ORA/'
cancer <- 'Renal_cell_carcinoma__NOS'
sircleFileName <- file.path(input_data_dir, paste0('sircle_PorMandR_', cancer, '.csv'))
# Run the ORA
#Import the pathways:
KEGG <- read.csv(file.path(GSEA_file_dir, "c2.cp.kegg.v6.2.symbols.csv"))
Reactome <- read.csv(file.path(GSEA_file_dir, "c2.cp.reactome.v6.2.symbols.csv"))
Biocarta <- read.csv(file.path(GSEA_file_dir, "c2.cp.biocarta.v6.2.symbols.csv"))
Hallmarks <- read.csv(file.path(GSEA_file_dir, "h.all.v6.2.symbols.csv"))
GO_BP <- read.csv(file.path(GSEA_file_dir, "c5.go.bp.v7.2.symbols.csv"))
GO_CC <- read.csv(file.path(GSEA_file_dir, "c5.go.cc.v7.2.symbols.csv"))
GO_MF <- read.csv(file.path(GSEA_file_dir, "c5.go.mf.v7.2.symbols.csv"))
Metabolic <- read.csv(file.path(GSEA_file_dir, "41467_2016_BFncomms13041_MOESM340_ESM.csv"))
Correction <- read.csv(file.path(GSEA_file_dir, "41467_2016_BFncomms13041_MOESM340_ESM.csv"))
#Run the GSEA analysis
##1."KEGG", "Reactome", "Biocarta", "Hallmarks"
pathways <- rbind(KEGG, Reactome, Biocarta, Hallmarks)
pathway_list <- list()
for (pathway in unique(pathways$term)) {
pathway_list[[pathway]] <- as.character(pathways[pathways$term == pathway, 1])
}
sircleORAHuman <- function(filename, entrezId, title, regLabels="RegulatoryLabels", emptyRegLabel="", fileType="pdf",
minGSSize=10, qvalueCutoff=0.2, pvalueCutoff=0.05, showCatagory=30, outputFolder='') {
## ------------ Setup and installs ----------- ##
packages <- c("org.Hs.eg.db", "clusterProfiler", "svglite", "enrichplot")
install.packages(setdiff(packages, rownames(installed.packages())))
library(org.Hs.eg.db)
library(clusterProfiler)
library(svglite)
library(enrichplot)
## ------------ Run ----------- ##
# open the data
df <- read.csv(filename)
allGenes <- as.character(df[[entrezId]]) #
clusterGenes <- subset(df, ! df[[regLabels]] == emptyRegLabel)
clusterGenes <- subset(clusterGenes, ! clusterGenes[[regLabels]] == "Not-Background")
grps_labels <- unlist(unique(clusterGenes[regLabels]))
for(g in grps_labels) {
grpGenes <- subset(df, df[[regLabels]] == g)
print(g)
print(dim(grpGenes))
clusterGo <- enrichGO(gene = as.character(grpGenes[[entrezId]]),
universe = allGenes,
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
qvalueCutoff = 1.0,
pvalueCutoff = 1.0,
minGSSize = minGSSize,
readable = TRUE)
# We have a cutoff of all of them, and then only visualise the ones that the user wants...
clusterGoSummary <- data.frame(clusterGo)
write.csv(clusterGoSummary, paste(outputFolder, 'ClusterGoSummary_', g, title, '.csv', sep=""))#Export the ORA results as .csv
if (!(dim(clusterGoSummary)[1] == 0)) {#exclude df's that have no observations
Dotplot <- dotplot(clusterGo, showCategory=showCatagory) +
ggtitle(paste("Dotplot ", g, sep=""))
ggsave(file=paste(outputFolder, "SiRCle-ORA_Dotplot_Human_", g, title, ".", fileType, sep=""), plot=Dotplot, width=10, height=8)
x2 <- pairwise_termsim(clusterGo)
Emapplot <- emapplot(x2, pie_scale=1.5, layout = "nicely")+
ggtitle(paste("Emapplot ", g, sep=""))
ggsave(file=paste(outputFolder, "SiRCle-ORA_Emapplot_Human_", g, title, ".", fileType, sep="" ), plot=Emapplot, width=10, height=8)
Heatplot <- heatplot(clusterGo,showCategory=showCatagory) +
theme(axis.text.x =element_text(size=5), axis.text.y =element_text(size=8,face="bold"), axis.title=element_text(size=12,face="bold"))+
ggtitle(paste("Heatplot ", g, sep=""))
ggsave(file=paste(outputFolder, "SiRCle-ORA_Heatplot_Human_", g,title, ".", fileType, sep="" ), plot=Heatplot, width=10, height=8)
}
}
}
sircleFileName <- file.path(input_data_dir, paste0('Protein_', cancer, '_ORA.csv'))
test_title <- '_protein'
sircleORAHuman(sircleFileName, "entrezgene_id", test_title, 'ORA_label', emptyRegLabel="None", minGSSize=10, qvalueCutoff=0.1, outputFolder=output_data_dir)
sircleFileName <- file.path(input_data_dir, paste0('RNA_', cancer, '_ORA.csv'))
test_title <- '_RNA'
sircleORAHuman(sircleFileName, "entrezgene_id", test_title, 'ORA_label', emptyRegLabel="None", minGSSize=10, qvalueCutoff=0.1, outputFolder=output_data_dir)
sircleFileName <- file.path(input_data_dir, paste0('DNAMethylation_', cancer, '_ORA.csv'))
test_title <- '_Methylation'
sircleORAHuman(sircleFileName, "entrezgene_id", test_title, 'ORA_label', emptyRegLabel="None", minGSSize=10, qvalueCutoff=0.1, outputFolder=output_data_dir)
In [25]:
from matplotlib_venn import venn3
ora_dir = os.path.join(output_dir, 'ORA_EarlyLate')
# Get the overlapping GO terms that were significant
p_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_DOWN_protein.csv'))
p_go = p_go[p_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
m_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_Hyper_Methylation.csv'))
m_go = m_go[m_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
r_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_DOWN_RNA.csv'))
r_go = r_go[r_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
venn3([set(p_go), set(r_go), set(m_go)], set_labels=['Proteomics', 'RNAseq', 'DNA Methylation'],
set_colors=('#0101d7', '#e68e25', '#6aaf44'),)
plt.rcParams['svg.fonttype'] = 'none'
plt.title('Down')
plt.savefig(os.path.join(fig_dir, f'Figure4B_DownGO_EarlyLate.svg'))
plt.show()
# Get the overlapping GO terms that were significant
p_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_UP_protein.csv'))
p_go = p_go[p_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
m_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_Hypo_Methylation.csv'))
m_go = m_go[m_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
r_go = pd.read_csv(os.path.join(ora_dir, f'ClusterGoSummary_UP_RNA.csv'))
r_go = r_go[r_go['p.adjust'] < 0.05]['Description'].values # Get sig terms
venn3([set(p_go), set(r_go), set(m_go)], set_labels=['Proteomics', 'RNAseq', 'DNA Methylation'],
set_colors=('#0101d7', '#e68e25', '#6aaf44'),)
plt.rcParams['svg.fonttype'] = 'none'
plt.title('UP')
plt.savefig(os.path.join(fig_dir, 'Figure4B_UPGO_EarlyLate.svg'))
In [26]:
import seaborn as sns
import matplotlib as mpl
import matplotlib.cm as cm
import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (3,10)
sns.set_theme(style="whitegrid")
cmaps = {'Down regulated': {'CpG': 'darkorange', 'Protein': 'royalblue', 'RNA': 'teal'},
'Up regulated': {'CpG': 'limegreen', 'Protein': 'darkred', 'RNA': 'lightcoral'},
}
for c in comparisons:
combined_df = pd.DataFrame()
rna_df = pd.read_csv(os.path.join(ora_dir, f'{comparisons[c][0]}'))
rna_df['GeneRatio'] = [int(g.split('/')[0])/int(g.split('/')[1]) for g in rna_df['GeneRatio'].values]
rna_df['Data'] = 'RNA'
rna_df = rna_df[rna_df['p.adjust'] < 0.05]
rna_df.sort_values('p.adjust', inplace=True)
protein_df = pd.read_csv(os.path.join(ora_dir, f'{comparisons[c][1]}'))
protein_df['GeneRatio'] = [int(g.split('/')[0])/int(g.split('/')[1]) for g in protein_df['GeneRatio'].values]
protein_df['Data'] = 'Protein'
protein_df = protein_df[protein_df['p.adjust'] < 0.05]
protein_df.sort_values('p.adjust', inplace=True)
cpg_df = pd.read_csv(os.path.join(ora_dir, f'{comparisons[c][2]}'))
cpg_df['GeneRatio'] = [int(g.split('/')[0])/int(g.split('/')[1]) for g in cpg_df['GeneRatio'].values]
cpg_df['Data'] = 'CpG'
cpg_df = cpg_df[cpg_df['p.adjust'] < 0.05]
cpg_df.sort_values('p.adjust', inplace=True)
combined_df = pd.concat([combined_df, rna_df])
combined_df = pd.concat([combined_df, cpg_df])
combined_df = pd.concat([combined_df, protein_df])
# Concat all together
num = 10
rna_df = pd.concat([rna_df.head(num), cpg_df.head(num), protein_df.head(num)])
rna_df.sort_values('p.adjust', inplace=True)
max_c = max(rna_df['Count'].values)
min_c = min(rna_df['Count'].values)
size = [int(1 + ((int(c) - min_c)/(max_c - min_c))*300) for c in rna_df['Count'].values]
norm = mpl.colors.Normalize(vmin=np.min(rna_df['p.adjust'].values), vmax=np.max(rna_df['p.adjust'].values))
cmap = cm.RdBu
m = cm.ScalarMappable(norm=norm, cmap=cmap)
c_map = cmaps[c]
colours = [c_map.get(v) for v in rna_df['Data'].values]
plt.scatter(rna_df['GeneRatio'].values, rna_df['Description'].values, s=size, c=colours)#c=m.to_rgba(rna_df['p.adjust'].values))
plt.title(c)
combined_df.to_csv(f'{output_dir}{c}_Late-Early-Tumour.csv', index=False)
gmin = plt.scatter([], [], s=10 + int(300*(min_c/(max_c - min_c))), marker='o', color='#222')
#gmid = plt.scatter([], [], s=int(10 + 300*(((max_c - min_c)/2)/max_c)), marker='o', color='#222')
gmax = plt.scatter([], [], s=10 + int(300*(1)), marker='o', color='#222')
legend = plt.legend((gmin, gmax),
(str(min_c), str(max_c)),
scatterpoints=1,
loc='lower left',
ncol=1,
fontsize=8, bbox_to_anchor=(0, -0.1))
legend.set_title("No. Genes")
plt.gca().add_artist(legend)
plt.savefig(f'{output_dir}DotPlot_GO_{c}_EarlyLate.svg')
plt.show()
#sns.stripplot(data=rna_df.head(30), y='Description', x='GeneRatio', size=size) #, color='p.adjust', size='Count')
In [ ]: