bensard_figures_2019

MPC EXPRESSION ANALYSIS NOTEBOOK

Paper: Regulation of Tumor Initiation by the Mitochondrial Pyruvate Carrier
Paper Authors: Claire L. Bensard and Dona R. Wisidigama, et. al.

Analysis performed within the scope of this notebook by Jordan A. Berg and Alex J. Bott
Description: The following notebook contains the code required to replicate specified plots and supplements

Import Dependencies

import os
import re
import pandas as pd 
import xpressplot as xp
from scipy import stats
from math import sqrt
from statistics import mean, stdev
from sklearn import preprocessing
from statsmodels.stats.multitest import multipletests
import GEOparse
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

#Some microarray data was also processed using Alt Analyze (http://www.altanalyze.org/)

#Set path to this directory for accessing and saving files
__path__  = os.getcwd() + '/'

Import Human Microarray Dataset GSE8671

The GEO-accessible data is not normalized. We therefore used Alt Analyze (http://www.altanalyze.org/) to RMA normalize probe signal.

#Get data
df_GSE8671 = xp.get_df(__path__ + 'data/GSE8671_rma_normalized.csv', delimiter=',') #RMA normalized with Alt Analyze
info_GSE8671 = xp.get_info(__path__ + "data/sample_info_gse8671.csv", delimiter=',')
df_GSE8671_c = xp.keep_labels(df_GSE8671, info_GSE8671, label_list=['Normal','Adenoma'])
df_GSE8671_clean = xp.clean_df(df_GSE8671_c)

#Collapse multi-mapping probes
df_GSE8671_collapsed = xp.probe_collapse(df_GSE8671_clean, __path__ + "data/GPL570.txt")
df_GSE8671_collapsed.to_csv(__path__ + "data/collapsed_GSE8671.csv",sep=',')

#Scale dataset
df_GSE8671_scaled, df_GSE8671_labeled = xp.prep_data(df_GSE8671_collapsed, info_GSE8671)

#sort by sample labels
info_GSE8671_sorted = info_GSE8671.copy()
info_GSE8671_sorted = info_GSE8671_sorted.loc[info_GSE8671_sorted[1].isin(['Adenoma', 'Normal'])]
info_GSE8671_sorted = info_GSE8671_sorted.sort_values([1], ascending=False)
info_GSE8671_sorted_list = info_GSE8671_sorted[0].tolist()
df_GSE8671_scaled_sorted = df_GSE8671_scaled[info_GSE8671_sorted_list]

#set palette
gse8671_colors = {'Adenoma': (0.8705882352941177, 0.5607843137254902, 0.0196078431372549),
        'Normal': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275)}

xp.check_samples(df_GSE8671_clean)

png

Import Human Microarray Dataset GSE20916

def get_geo(
    geo_id,
    output_info=False, output_path="./"):

    # Get data
    gse = GEOparse.get_GEO(geo=str(geo_id).upper(), destdir=output_path) # Import GSE dataset

    data = gse.pivot_samples('VALUE')
    data = xp.clean_df(data)

    # Get metadata
    # Write data to output file
    if output_info != False:
        with open(str(geo_id).upper() + '.txt', 'w+') as f: # Save all information as text file for reference
            for gsm_name, gsm in gse.gsms.items():
                f.write(gsm_name + '\n')
                for key, value in gsm.metadata.items():
                    f.write(" - %s : %s" % (key, ", ".join(value)) + '\n')

    # Populate metadata with sample ids and names
    metadata = pd.DataFrame(columns=['gsm', 'title']) # Create dataframe
    gsm_list, title_list, data_processing_list = [], [], []
    for gsm_name, gsm in gse.gsms.items():
        for key, value in gsm.metadata.items():
            if key == 'title':
                title_list.append(''.join(value))
            if key == 'geo_accession':
                gsm_list.append(''.join(value))
            if key == 'data_processing':
                data_processing_list.append(''.join(value))

    metadata['gsm'], metadata['title'] = gsm_list, title_list
    metadata.columns = range(metadata.shape[1])

    # Output processing style
    print('Data processing summary:\n' + str(set(data_processing_list))) # To determine if all samples have undergone the sample data processing

    # Clean data
    del data.columns.name
    del data.index.name

    # Clean metadata
    metadata[1] = metadata[1].apply(lambda x: x[0:(re.search("\d", x).start()) - 1])

    return data, metadata

# Get data from GEO database 
#df_GSE20916, info_GSE20916 = get_geo('GSE20916')

#df_GSE20916.to_csv(__path__ + "data/GSE20916_normalized.csv", sep=',')
#info_GSE20916.to_csv(__path__ + "data/sample_info_gse20916.csv", sep=',', index=False)

#Get data
df_GSE20916 = pd.read_csv(__path__ + "data/GSE20916_normalized.csv", index_col=0)
info_GSE20916 = pd.read_csv(__path__ + "data/sample_info_gse20916.csv")
info_GSE20916.columns = [0,1]
info_GSE20916[1] = info_GSE20916[1].str.capitalize() #Make sample types look nice
info_GSE20916 = info_GSE20916.replace('Normal_colon', 'Normal')
df_GSE20916_c = xp.keep_labels(df_GSE20916, info_GSE20916, label_list=['Normal','Adenoma','Adenocarcinoma'])
df_GSE20916_clean = xp.clean_df(df_GSE20916_c)

#Collapse multi-mapping probes
df_GSE20916_collapsed = xp.probe_collapse(df_GSE20916_clean, __path__ + "data/GPL570.txt")
df_GSE20916_collapsed.to_csv(__path__ + "data/collapsed_GSE20916.txt",sep='\t')

#Scale sorted dataset
df_GSE20916_scaled, df_GSE20916_labeled = xp.prep_data(df_GSE20916_collapsed, info_GSE20916)

#sort by sample labels
info_GSE20916_sorted = info_GSE20916.copy()
info_GSE20916_sorted = info_GSE20916_sorted.loc[info_GSE20916_sorted[1].isin(['Adenoma', 'Adenocarcinoma','Normal'])]
info_GSE20916_sorted = info_GSE20916_sorted.sort_values([1], ascending=False)
info_GSE20916_sorted_list = info_GSE20916_sorted[0].tolist()
df_GSE20916_scaled_sorted = df_GSE20916_scaled[info_GSE20916_sorted_list]

#Drop Adenocarcinomas
df_GSE20916_collapsed_noac = xp.drop_label(df_GSE20916_collapsed, info_GSE20916, 'Adenocarcinoma')

gse20916_colors = {'Adenocarcinoma': (0.5725490196078431, 0.5843137254901961, 0.5686274509803921),
        'Adenoma': (0.8705882352941177, 0.5607843137254902, 0.0196078431372549),
        'Normal': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275)}

xp.check_samples(df_GSE20916_clean)

png

Import Mouse APC Data

apc_df = pd.read_csv(__path__ + "data/custom_apc_normalized.csv", sep=",", index_col=0, low_memory=False)
apc_info = xp.get_info(__path__ + "data/apc_info.csv", delimiter=',')
apc_clean = xp.clean_df(apc_df)
apc_scaled, apc_labeled = xp.prep_data(apc_clean, apc_info)

apc_colors = {'APC-KO': (0.00392156862745098, 0.45098039215686275, 0.6980392156862745),
         'APC-WT': (0.8705882352941177, 0.5607843137254902, 0.0196078431372549),
         'Normal-KO': (0.8352941176470589, 0.3686274509803922, 0.0),
         'Normal-WT': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275)}

Import Mouse AOMDSS Data

aomdss_df = pd.read_csv(__path__ + "data/custom_aomdss_normalized.csv", sep=",", index_col=0, low_memory=False)
aomdss_info = xp.get_info(__path__ + "data/aomdss_info.csv", delimiter=',')
aomdss_clean = xp.clean_df(aomdss_df)
aomdss_scaled, aomdss_labeled = xp.prep_data(aomdss_clean, aomdss_info)

aomdss_colors = {'AOMDSS-KO': (0.00392156862745098, 0.45098039215686275, 0.6980392156862745),
         'AOMDSS-WT': (0.8705882352941177, 0.5607843137254902, 0.0196078431372549),
         'Normal-KO': (0.8352941176470589, 0.3686274509803922, 0.0),
         'Normal-WT': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275)}

Establish Gene Sets to Survey

with open(__path__ + "data/custom22.csv", 'r') as f:
    for line in f:
        custom22_mouse = line.split(",")
custom22_human = [x.upper() for x in custom22_mouse]

stem_genes = ['Ctnnb1','Notch1','Ascl2','Myc','Hopx','Sox9','Lgr5','Lef1','Mmp7','Axin2','Cd44','Ccnd1','Bmi1','Tert']
diff_genes = ['Fabp2','Atoh1','Muc2','Krt20','Chga','Vil1','MPC1','MPC2']

pyru = pd.read_csv('data/custom_pyruvate_list.csv',header=None)
pyru_list = pyru[0].tolist()

# Display number of missing genes from the microarray datasets from the gene list

print('Total pyruvate-related genes:',len(pyru_list))
print('Total in GSE8671:', df_GSE8671_labeled.reindex(pyru_list).dropna().shape[0])
print('Total in GSE20916:', df_GSE20916_labeled.reindex(pyru_list).dropna().shape[0])

Total pyruvate-related genes: 80
Total in GSE8671: 75
Total in GSE20916: 60

GSE8671_pyru = df_GSE8671_labeled.reindex(pyru_list)
GSE8671_pyru.to_csv('./data/GSE8671_pyruvate_list.csv')

GSE20916_pyru = df_GSE20916_labeled.reindex(pyru_list)
GSE20916_pyru.to_csv('./data/GSE20916_pyruvate_list.csv')

LDH and MPC Expression with Stats

"""
LDH and MPC Relationship in Human Dataset GSE8671
"""
xp.gene_overview(df_GSE8671_labeled, info_GSE8671, gene_name='MPC1',palette=gse8671_colors, 
                  order=['Normal','Adenoma'], grid=True, whitegrid=True,
                 save_fig=__path__ + 'plots/MPC1_boxswarm_GSE8671.pdf')

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE8671_labeled, info_GSE8671, gene_name='MPC2',palette=gse8671_colors, 
                  order=['Normal','Adenoma'],
                 save_fig=__path__ + 'plots/MPC2_boxswarm_GSE8671.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE8671_labeled, info_GSE8671, gene_name='LDHA',palette=gse8671_colors, 
                  order=['Normal','Adenoma'],
                 save_fig=__path__ + 'plots/LDHA_boxswarm_GSE8671.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE8671_labeled, info_GSE8671, gene_name='LDHB',palette=gse8671_colors, 
                  order=['Normal','Adenoma'],
                 save_fig=__path__ + 'plots/LDHB_boxswarm_GSE8671.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

#Statistical tests
GSE8671_stats = df_GSE8671_labeled.copy()
GSE8671_norm = GSE8671_stats[info_GSE8671.loc[info_GSE8671[1] == 'Normal'][0].tolist()]
GSE8671_ad = GSE8671_stats[info_GSE8671.loc[info_GSE8671[1] == 'Adenoma'][0].tolist()]

p_vals = []
genes = GSE8671_stats.index.tolist()

total = len(genes) - 1
counter = 1

for x in genes[1:]:
    s, p = stats.ranksums(xp.get_array(GSE8671_norm, str(x)), xp.get_array(GSE8671_ad, str(x)))
    p_vals.append(p) 
    print(counter,'/',total)
    counter += 1

fdr = multipletests(p_vals, alpha=0.05, method='fdr_bh')

sig_data = pd.DataFrame()
sig_data[0] = p_vals 
sig_data[1] = fdr[1]
sig_data.index = genes[1:]

print('GSE8671 FDR for Normal vs Adenoma')

#MPC1
print('MPC1: FDR=' + str(sig_data.loc['MPC1'][1]))

#MPC2
print('MPC2: FDR=' + str(sig_data.loc['MPC2'][1]))

#LDHA
print('LDHA: FDR=' + str(sig_data.loc['LDHA'][1]))

#LDHB
print('LDHB: FDR=' + str(sig_data.loc['LDHB'][1]))

GSE8671 FDR for Normal vs Adenoma
MPC1: FDR=1.2877272621302151e-09
MPC2: FDR=9.321124459705085e-08
LDHA: FDR=1.0004823540206768e-07
LDHB: FDR=7.549800011433636e-08

"""
LDH and MPC Relationship in Human Dataset GSE20916
"""
xp.gene_overview(df_GSE20916_labeled, info_GSE20916, gene_name='MPC1',palette=gse20916_colors, 
                  order=['Normal','Adenoma','Adenocarcinoma'],
                 save_fig=__path__ + 'plots/MPC1_boxswarm_GSE20916.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE20916_labeled, info_GSE20916, gene_name='MPC2',palette=gse20916_colors, 
                  order=['Normal','Adenoma','Adenocarcinoma'],
                 save_fig=__path__ + 'plots/MPC2_boxswarm_GSE20916.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE20916_labeled, info_GSE20916, gene_name='LDHA',palette=gse20916_colors, 
                  order=['Normal','Adenoma','Adenocarcinoma'],
                 save_fig=__path__ + 'plots/LDHA_boxswarm_GSE20916.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

xp.gene_overview(df_GSE20916_labeled, info_GSE20916, gene_name='LDHB',palette=gse20916_colors, 
                  order=['Normal','Adenoma','Adenocarcinoma'],
                 save_fig=__path__ + 'plots/LDHB_boxswarm_GSE20916.pdf', grid=True, whitegrid=True)

<Figure size 432x288 with 0 Axes>

png

#Statistical tests
GSE20916_stats = df_GSE20916_labeled.copy()
GSE20916_norm = GSE20916_stats[info_GSE20916.loc[info_GSE20916[1] == 'Normal'][0].tolist()]
GSE20916_ad = GSE20916_stats[info_GSE20916.loc[info_GSE20916[1] == 'Adenoma'][0].tolist()]
GSE20916_ac = GSE20916_stats[info_GSE20916.loc[info_GSE20916[1] == 'Adenocarcinoma'][0].tolist()]

p_vals_ad = []
p_vals_ac = []
genes = GSE20916_stats.index.tolist()

total = len(genes) - 1
counter = 1

for x in genes[1:]:
    s_ad, p_ad = stats.ranksums(xp.get_array(GSE20916_norm, str(x)), xp.get_array(GSE20916_ad, str(x)))
    p_vals_ad.append(p_ad) 
    
    s_ad, p_ac = stats.ranksums(xp.get_array(GSE20916_norm, str(x)), xp.get_array(GSE20916_ac, str(x)))
    p_vals_ac.append(p_ac)
    
    print(counter,'/',total)
    counter += 1

fdr_ad = multipletests(p_vals_ad, alpha=0.05, method='fdr_bh')
fdr_ac = multipletests(p_vals_ac, alpha=0.05, method='fdr_bh')

sig_data = pd.DataFrame()
sig_data[0] = p_vals_ad 
sig_data[1] = fdr_ad[1]
sig_data[2] = p_vals_ac
sig_data[3] = fdr_ac[1]

sig_data.index = genes[1:]

print('GSE20916 FDR for Normal vs Adenoma')

#MPC1
print('MPC1: FDR=' + str(sig_data.loc['MPC1'][1]))

#MPC2
print('MPC2: FDR=' + str(sig_data.loc['MPC2'][1]))

#LDHA
print('LDHA: FDR=' + str(sig_data.loc['LDHA'][1]))

#LDHB
print('LDHB: FDR=' + str(sig_data.loc['LDHB'][1]))

print('')
print('GSE20916 FDR for Normal vs Adenocarcinoma')

#MPC1
print('MPC1: FDR=' + str(sig_data.loc['MPC1'][3]))

#MPC2
print('MPC2: FDR=' + str(sig_data.loc['MPC2'][3]))

#LDHA
print('LDHA: FDR=' + str(sig_data.loc['LDHA'][3]))

#LDHB
print('LDHB: FDR=' + str(sig_data.loc['LDHB'][3]))

GSE20916 FDR for Normal vs Adenoma
MPC1: FDR=9.791300308807664e-08
MPC2: FDR=1.7172807817415693e-05
LDHA: FDR=0.0002574743207763419
LDHB: FDR=2.8134953763073485e-07

GSE20916 FDR for Normal vs Adenocarcinoma
MPC1: FDR=4.746379251491133e-06
MPC2: FDR=0.06426095127053551
LDHA: FDR=1.474973692220452e-05
LDHB: FDR=1.896386031707828e-05

Pyruvate Metabolism in Human Microarrays and Mouse Nanostring Experiments

"""
Pyruvate-related Gene Expression Changes in Human Dataset GSE8671
"""
xp.heatmap(df_GSE8671_scaled_sorted, info_GSE8671, sample_palette=gse8671_colors, gene_list=pyru_list, 
            figsize=(14,8), save_fig=(__path__ + 'plots/GSE8671_customPyruvate_heatmap_colclustered.pdf'),
            row_cluster=True, col_cluster=True, yticklabels=False, xticklabels=False,
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

xp.heatmap(df_GSE8671_scaled_sorted, info_GSE8671, sample_palette=gse8671_colors, gene_list=pyru_list, 
            figsize=(14,8), save_fig=(__path__ + 'plots/GSE8671_customPyruvate_heatmap.pdf'),
            row_cluster=True, col_cluster=False, yticklabels=False, xticklabels=False,
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

"""
Pyruvate-related Gene Expression Changes in Human Dataset GSE20916
"""
xp.heatmap(df_GSE20916_scaled_sorted, info_GSE20916, sample_palette=gse20916_colors, gene_list=pyru_list, 
            figsize=(30,8), save_fig=(__path__ + 'plots/GSE20916_customPyruvate_heatmap_colclustered.pdf'),
            row_cluster=True, col_cluster=True, yticklabels=False, xticklabels=False,
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

xp.heatmap(df_GSE20916_scaled_sorted, info_GSE20916, sample_palette=gse20916_colors, gene_list=pyru_list, 
            figsize=(30,8), save_fig=(__path__ + 'plots/GSE20916_customPyruvate_heatmap.pdf'),
            row_cluster=True, col_cluster=False, yticklabels=False, xticklabels=False,
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

"""
Pyruvate-related Gene Expression Changes in Human Dataset GSE8671 -- PCA
"""
xp.pca(df_GSE8671_scaled, info_GSE8671, palette=gse8671_colors,
        gene_list=pyru_list,
       save_scree=__path__ + 'plots/GSE8671_pyru_list_PCA_scree.pdf',
       title='GSE8671 Pyruvate-related Genes PCA', 
        save_fig=__path__ + 'plots/GSE8671_pyru_list_2CI_PCA.pdf')

png

"""
Overall Gene Expression Changes in Human Dataset GSE8671 -- PCA
"""
xp.pca(df_GSE8671_scaled, info_GSE8671, palette=gse8671_colors,
       save_scree=__path__ + 'plots/GSE8671_all_genes_PCA_scree.pdf',
       title='GSE8671 All Genes PCA', 
        save_fig=__path__ + 'plots/GSE8671_all_genes_2CI_PCA.pdf')

png

"""
Pyruvate-related Gene Expression Changes in Human Dataset GSE20916 -- PCA
"""
xp.pca(df_GSE20916_scaled, info_GSE20916, palette=gse20916_colors, order_legend=[1,3,2],
        gene_list=pyru_list,
       save_scree=__path__ + 'plots/GSE20916_pyru_list_PCA_scree.pdf',
       title='GSE20916 Pyruvate-related Genes PCA', 
        save_fig=__path__ + 'plots/GSE20916_pyru_list_2CI_PCA.pdf')

png

"""
Overall Gene Expression Changes in Human Dataset GSE20916 -- PCA
"""
xp.pca(df_GSE20916_scaled, info_GSE20916, palette=gse20916_colors, order_legend=[1,3,2],
       save_scree=__path__ + 'plots/GSE20916_all_genes_PCA_scree.pdf',
       title='GSE20916 All Genes PCA', 
        save_fig=__path__ + 'plots/GSE20916_all_genes_2CI_PCA.pdf')

png

Mouse Nanostring and Human Microarray Heatmaps – Custom22 Probe Genes

"""
Mouse AOMDSS Heatmap Custom22 Probes
"""
xp.heatmap(aomdss_scaled, aomdss_info, sample_palette=aomdss_colors, gene_list=custom22_mouse, 
            figsize=(16,6.5), save_fig=(__path__ + 'plots/AOMDSS_custom22_heatmap.pdf'), 
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

"""
Mouse APC Heatmap Custom22 Probes
"""
xp.heatmap(apc_scaled, apc_info, sample_palette=apc_colors, gene_list=custom22_mouse, 
            figsize=(16,6.5), save_fig=(__path__ + 'plots/APC_custom22_heatmap.pdf'), 
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

"""
Human GSE8671 Heatmap Custom22 Probes
"""
xp.heatmap(df_GSE8671_scaled, info_GSE8671, sample_palette=gse8671_colors, gene_list=custom22_human, 
            figsize=(22,8), save_fig=(__path__ + 'plots/GSE8671_custom22_heatmap.pdf'), 
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

"""
Human GSE20916 Heatmap Custom22 Probes
"""
xp.heatmap(df_GSE20916_scaled, info_GSE20916, sample_palette=gse20916_colors, gene_list=custom22_human, 
            figsize=(40,8), save_fig=(__path__ + 'plots/GSE20916_custom22_heatmap.pdf'), 
            cbar_kws={'label': 'z-score'}, linewidths=0)

<Figure size 432x288 with 0 Axes>

png

Mouse Model PCA Plots

"""
Mouse AOMDSS Heatmap Custom22 Probes -- PCA
"""
xp.pca(aomdss_scaled, aomdss_info, palette=aomdss_colors, order_legend=[3,4,2,1], n_components=2,
        gene_list=custom22_human, whitegrid=True, grid=False,
       save_scree=__path__ + 'plots/AOMDSS_custom22_PCA_scree.pdf',
       title='AOMDSS Custom22 Probes PCA', 
        save_fig=__path__ + 'plots/AOMDSS_custom22_2CI_PCA.pdf')

png

"""
Mouse APC Heatmap Custom22 Probes -- PCA
"""
xp.pca(apc_scaled, apc_info, palette=apc_colors, order_legend=[2,3,1,4], n_components=2,
        gene_list=custom22_human, whitegrid=True, grid=False,
       save_scree=__path__ + 'plots/APC_custom22_PCA_scree.pdf',
       title='APC Custom22 Probes PCA', 
        save_fig=__path__ + 'plots/APC_custom22_2CI_PCA.pdf')

png

"""
Human GSE8671 Custom22 Probes -- PCA
"""
xp.pca(df_GSE8671_scaled, info_GSE8671, palette=gse8671_colors, n_components=2,
        gene_list=custom22_human, whitegrid=True, grid=False,
       save_scree=__path__ + 'plots/GSE8671_custom22_PCA_scree.pdf',
       title='GSE8671 Custom22 Probes PCA', 
        save_fig=__path__ + 'plots/GSE8671_custom22_2CI_PCA.pdf')

png

"""
Human GSE20916 Custom22 Probes -- PCA
"""
xp.pca(df_GSE20916_scaled, info_GSE20916, palette=gse20916_colors, order_legend=[1,3,2], n_components=2,
        gene_list=custom22_human, whitegrid=True, grid=False,
       save_scree=__path__ + 'plots/GSE20916_custom22_PCA_scree.pdf',
       title='GSE20916 Custom22 Probes PCA', 
        save_fig=__path__ + 'plots/GSE20916_custom22_2CI_PCA.pdf')

png

"""
Human GSE8671 Jointplots -- MPC1 vs AXIN2
"""
xp.jointplot(df_GSE8671_collapsed, info_GSE8671, 'MPC1', 'AXIN2', palette=gse8671_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE8671_MPC1_AXIN2_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

"""
Human GSE8671 Jointplots -- MPC1 vs SOX9
"""
xp.jointplot(df_GSE8671_collapsed, info_GSE8671, 'MPC1', 'SOX9', palette=gse8671_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE8671_MPC1_SOX9_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

"""
Human GSE20916 Jointplots -- MPC1 vs AXIN2, no adenocarcinomas
"""
xp.jointplot(df_GSE20916_collapsed_noac, info_GSE20916, 'MPC1', 'AXIN2', palette=gse20916_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE20916_MPC1_AXIN2_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

"""
Human GSE20916 Jointplots -- MPC1 vs SOX9, no adenocarcinomas
"""
xp.jointplot(df_GSE20916_collapsed_noac, info_GSE20916, 'MPC1', 'SOX9', palette=gse20916_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE20916_MPC1_SOX9_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

"""
Human GSE20916 Jointplots -- MPC1 vs AXIN2
"""
xp.jointplot(df_GSE20916_collapsed, info_GSE20916, 'MPC1', 'AXIN2', palette=gse20916_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE20916_MPC1_AXIN2_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

"""
Human GSE20916 Jointplots -- MPC1 vs SOX9
"""
xp.jointplot(df_GSE20916_collapsed, info_GSE20916, 'MPC1', 'SOX9', palette=gse20916_colors, 
              order=['Normal','Adenoma'],
              save_fig=__path__ + 'plots/GSE20916_MPC1_SOX9_jointplot.pdf', 
              dpi=600, bbox_inches='tight', whitegrid=True, grid=True)

<Figure size 432x288 with 0 Axes>

<Figure size 720x720 with 0 Axes>

png

<Figure size 432x288 with 0 Axes>

Microarray Correlations

def add_fdr(file):
    
    data = pd.read_csv(
        file,
        sep=',')
    p_vals = data['p_value'].values.tolist()
    fdr = multipletests(p_vals, alpha=0.05, method='fdr_bh')
    data['fdr'] = fdr[1]
    
    data.to_csv(file)

    
#GSE20916/MPC1/ALL
xp.linreg(df_GSE20916_collapsed, 'MPC1', __path__ + 'data/GSE20916_MPC1_ALL_correlations.csv')
add_fdr(__path__ + 'data/GSE20916_MPC1_ALL_correlations.csv')

#GSE20916/MPC2/ALL
xp.linreg(df_GSE20916_collapsed, 'MPC2', __path__ + 'data/GSE20916_MPC2_ALL_correlations.csv')
add_fdr(__path__ + 'data/GSE20916_MPC2_ALL_correlations.csv')

#GSE20916/MPC1/NoAdenocarcinoma
df_GSE20916_collapsed_noAC = xp.drop_label(df_GSE20916_collapsed, info_GSE20916, 'Adenocarcinoma')
xp.linreg(df_GSE20916_collapsed_noAC, 'MPC1', __path__ + 'data/GSE20916_MPC1_noAC_correlations.csv')
add_fdr(__path__ + 'data/GSE20916_MPC1_noAC_correlations.csv')

#GSE20916/MPC2/NoAdenocarcinoma
df_GSE20916_collapsed_noAC = xp.drop_label(df_GSE20916_collapsed, info_GSE20916, 'Adenocarcinoma')
xp.linreg(df_GSE20916_collapsed_noAC, 'MPC2', __path__ + 'data/GSE20916_MPC2_noAC_correlations.csv')
add_fdr(__path__ + 'data/GSE20916_MPC2_noAC_correlations.csv')

#GSE8671/MPC1/ALL
xp.linreg(df_GSE8671_collapsed, 'MPC1', __path__ + 'data/GSE8671_MPC1_ALL_correlations.csv')
add_fdr(__path__ + 'data/GSE8671_MPC1_ALL_correlations.csv')

#GSE8671/MPC2/ALL
xp.linreg(df_GSE8671_collapsed, 'MPC2', __path__ + 'data/GSE8671_MPC2_ALL_correlations.csv')
add_fdr(__path__ + 'data/GSE8671_MPC2_ALL_correlations.csv')

TCGA GTEx Survey of MPC and LDH Expression

See doi: 10.1038/sdata.2018.61
Samples with z-score for gene set > 5 were excluded from heatmap to avoid signal swamp-out

# Import TCGA/GTEx Unity Data 

data_dir = '/Users/jordan/Desktop/project_trix/unified_data/unified_normalized/'
file_list = ['colon-rsem-fpkm-gtex.txt',
            'coad-rsem-fpkm-tcga.txt',
            'coad-rsem-fpkm-tcga-t.txt']
file_list = [str(data_dir) + str(x) for x in file_list]

#Initialize
df_list = []
meta = pd.DataFrame(columns = [0,1])
previous = 0
label_number = 0
drop_extra='Entrez_Gene_Id'
label_list = [
    'GTEX_Normal',
    'TCGA_Normal',
    'TCGA_Tumor']

#Import
for x in range(len(file_list)):

    #Get and clean data
    data = pd.read_csv(
        str(file_list[x]),
        index_col = 0,
        sep = '\t')
    del data.index.name

    #Drop extra column if needed
    if drop_extra != None:
        data = data.drop(labels=['Entrez_Gene_Id'], axis=1)

    #Add data to data list
    df_list.append(data)

    print('Size of ' + str(file_list[x]) + ': ' + str(data.shape))

    #Get metadata for current data matrix
    data_L = data.columns.tolist()

    for x in range(len(data_L)):
        meta.loc[previous + x] = [data_L[x], label_list[label_number]]

    #Update iterators
    previous = previous + len(data_L)
    label_number += 1

#Join dataframes
data = pd.concat(df_list, axis=1)
print('Size of final dataframe: ' + str(data.shape))
print('Size of final metadata table: ' + str(meta.shape))

data_scale = data.copy()
data_scale[data_scale.columns] = preprocessing.scale(data_scale[data_scale.columns], axis=1)

Size of /Users/jordan/Desktop/project_trix/unified_data/unified_normalized/colon-rsem-fpkm-gtex.txt: (19096, 339)
Size of /Users/jordan/Desktop/project_trix/unified_data/unified_normalized/coad-rsem-fpkm-tcga.txt: (19096, 41)
Size of /Users/jordan/Desktop/project_trix/unified_data/unified_normalized/coad-rsem-fpkm-tcga-t.txt: (19096, 285)
Size of final dataframe: (19096, 665)
Size of final metadata table: (665, 2)

sample_colors = {'GTEX_Normal': '#005f42',
                 'TCGA_Normal': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275),
                 'TCGA_Tumor': (0.5725490196078431, 0.5843137254901961, 0.5686274509803921)
                }

ldh = ['LDHA',
        'LDHB',
        'MPC1',
        'MPC2']
data_scale_ldh = data_scale.copy()
data_scale_ldh = data_scale_ldh.loc[ldh]

### To remove the constraint on expression magnitude, 
### remove the line below and run the code block
data_scale_ldh[data_scale_ldh > 5] = 5 # <-- Remove this one
 
xp.heatmap(
    data_scale_ldh,
    meta,
    sample_colors,
    col_cluster = False,
    xticklabels = False,
    figsize=(15,1.5))

colors = {'GTEx Normal': '#005f42',
          'TCGA Normal': (0.00784313725490196, 0.6196078431372549, 0.45098039215686275),
          'TCGA Tumor': (0.5725490196078431, 0.5843137254901961, 0.5686274509803921)
         }

g = lambda m,c: plt.plot([],[],marker='o', color=c, ls="none")[0]
handles_g = [g("s", list(colors.values())[i]) for i in range(len(list(colors.values())))]
plt.legend(handles_g, list(colors.keys()), bbox_to_anchor=(0, -7), loc=3, borderaxespad=0., title='Samples')

plt.savefig(__path__ + 'plots/tcga_mpc_ldh_heatmap.pdf',dpi=1800,bbox_inches='tight')
plt.show()

<Figure size 432x288 with 0 Axes>

png

# TCGA effect sizes and p-values 
# https://toptipbio.com/cohens-d/
# https://stackoverflow.com/questions/21532471/how-to-calculate-cohens-d-in-python
def cohen_d(comp, base):

    cohens_d = (mean(comp) - mean(base)) / (sqrt((stdev(comp) ** 2 + stdev(base) ** 2) / 2))

    return cohens_d

#Statistical tests
data_stats = data.copy()

normal_g = meta.loc[(meta[1] == 'GTEX_Normal')][0].tolist()
normal_t = meta.loc[(meta[1] == 'TCGA_Normal')][0].tolist()
normal = normal_g + normal_t
data_norm = data_stats[normal]

data_ad = data_stats[meta.loc[meta[1] == 'TCGA_Tumor'][0].tolist()]

p_vals = []
d_vals = []
genes = data_stats.index.tolist()

total = len(genes) - 1
counter = 1

for x in genes[1:]:
    s, p = stats.ranksums(xp.get_array(data_norm, str(x)), xp.get_array(data_ad, str(x)))
    p_vals.append(p) 
    d = cohen_d(xp.get_array(data_ad, str(x)), xp.get_array(data_norm, str(x)))
    d_vals.append(d)
    print(counter,'/',total)
    counter += 1

fdr = multipletests(p_vals, alpha=0.05, method='fdr_bh')

sig_data = pd.DataFrame()
sig_data[0] = p_vals 
sig_data[1] = fdr[1]
sig_data[2] = d_vals
sig_data.index = genes[1:]

print('GTEx/TCGA FDR for Normal vs Adenoma')

#MPC1
print('MPC1: FDR=' + str(sig_data.loc['MPC1'][1]))
print('\tD=' + str(sig_data.loc['MPC1'][2]))

#MPC2
print('MPC2: FDR=' + str(sig_data.loc['MPC2'][1]))
print('\tD=' + str(sig_data.loc['MPC2'][2]))

#LDHA
print('LDHA: FDR=' + str(sig_data.loc['LDHA'][1]))
print('\tD=' + str(sig_data.loc['LDHA'][2]))

#LDHB
print('LDHB: FDR=' + str(sig_data.loc['LDHB'][1]))
print('\tD=' + str(sig_data.loc['LDHB'][2]))

GTEx/TCGA FDR for Normal vs Adenoma
MPC1: FDR=7.705583486288622e-95
	D=-2.103597139871266
MPC2: FDR=3.7429732880948714e-32
	D=-0.9468630678233025
LDHA: FDR=1.369631843350729e-17
	D=0.6547615336558776
LDHB: FDR=2.8647008844045153e-27
	D=0.9163199268827382

This site is open source. Improve this page.