Figure. CNV eQTL Effect Sizes



In [1]:

    
import copy
import cPickle
import os
import subprocess

import cdpybio as cpb
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
pd.options.mode.chained_assignment = None  # default='warn'
import pybedtools as pbt
import scipy.stats as stats
import seaborn as sns

import ciepy
import cardipspy as cpy

%matplotlib inline
%load_ext rpy2.ipython

dy_name = 'figure_cnv_eqtl_effect_sizes'
    
outdir = os.path.join(ciepy.root, 'output', dy_name)
cpy.makedir(outdir)

private_outdir = os.path.join(ciepy.root, 'private_output', dy_name)
cpy.makedir(private_outdir)

import socket
if socket.gethostname() == 'fl-hn1' or socket.gethostname() == 'fl-hn2':
    dy = os.path.join(ciepy.root, 'sandbox', 'tmp', dy_name)
    cpy.makedir(dy)
    pbt.set_tempdir(dy)

Each figure should be able to fit on a single 8.5 x 11 inch page. Please do not send figure panels as individual files. We use three standard widths for figures: 1 column, 85 mm; 1.5 column, 114 mm; and 2 column, 174 mm (the full width of the page). Although your figure size may be reduced in the print journal, please keep these widths in mind. For Previews and other three-column formats, these widths are also applicable, though the width of a single column will be 55 mm.



In [2]:

    
fn = os.path.join(ciepy.root, 'output/cnv_analysis/cnv_gene_variants.pickle')
cnv_gv = cPickle.load(open(fn))
fn = os.path.join(ciepy.root, 'output/cnv_analysis/combined_info.pickle')
combined_info = cPickle.load(open(fn))

sig_cnvs = set(cnv_gv.cnv_id)
not_sig_cnvs = set(combined_info.index) - sig_cnvs



In [3]:

    
sns.set_style('whitegrid')



In [6]:

    
fig = plt.figure(figsize=(6.85, 4.5), dpi=300)

gs = gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0])
ax.text(0, 0, 'Figure S5',
        size=16, va='bottom')
ciepy.clean_axis(ax)
ax.set_xticks([])
ax.set_yticks([])
gs.tight_layout(fig, rect=[0, 0.90, 0.5, 1])

gs = gridspec.GridSpec(2, 2)

tdf = cnv_gv.sort_values(by=['cnv_overlaps_gene', 'pvalue'], 
                         ascending=[False, True]).drop_duplicates(subset=['gene_id'])

a = tdf[tdf.cnv_overlaps_gene_cons]
b = tdf[tdf.cnv_overlaps_gene_cons == False]

# Genic, lead CNV
ax = fig.add_subplot(gs[0, 0])
ax.set_ylabel('Number of genes', fontsize=8)
ax.set_xlabel('$\\beta$', fontsize=8)
a.beta.hist(bins=np.arange(-2.75, 3, 0.25), ax=ax)
ax.grid(axis='x')
for t in ax.get_xticklabels() + ax.get_yticklabels():
    t.set_fontsize(8)
print('{:,} lead genic CNVs.'.format(a.shape[0]))
p = stats.binom_test((a.beta > 0).value_counts())
print('Effect sizes for genic lead CNVs are biased '
      '(p={:.3e}, binomial test).'.format(p))

# Intergenic, lead CNV
ax = fig.add_subplot(gs[0, 1])
b.beta.hist(bins=np.arange(-2.75, 3, 0.25), ax=ax)
ax.set_ylabel('Number of genes', fontsize=8)
ax.set_xlabel('$\\beta$', fontsize=8)
ax.grid(axis='x')
for t in ax.get_xticklabels() + ax.get_yticklabels():
    t.set_fontsize(8)
p = stats.binom_test((b.beta > 0).value_counts())
print('{:,} lead intergenic CNVs.'.format(b.shape[0]))
print('Effect sizes for intergenic lead CNVs are biased '
      '(p={:.3e}, binomial test).'.format(p))
    
a = cnv_gv[cnv_gv.cnv_overlaps_gene_cons]
b = cnv_gv[cnv_gv.cnv_overlaps_gene_cons == False]
    
# Genic, all CNV associations
ax = fig.add_subplot(gs[1, 0])
ax.set_ylabel('Number of CNVs', fontsize=8)
ax.set_xlabel('$\\beta$', fontsize=8)
a.beta.hist(bins=np.arange(-2.75, 3, 0.25), ax=ax)
ax.grid(axis='x')
for t in ax.get_xticklabels() + ax.get_yticklabels():
    t.set_fontsize(8)
p = stats.binom_test((a.beta > 0).value_counts())
print('{:,} genic CNVs.'.format(a.shape[0]))
print('Effect sizes for all genic CNV eQTLs are biased '
      '(p={:.3e}, binomial test).'.format(p))

# Intergenic, all CNV associations
ax = fig.add_subplot(gs[1, 1])
b.beta.hist(bins=np.arange(-2.75, 3, 0.25), ax=ax)
ax.set_ylabel('Number of CNVs', fontsize=8)
ax.set_xlabel('$\\beta$', fontsize=8)
ax.grid(axis='x')
for t in ax.get_xticklabels() + ax.get_yticklabels():
    t.set_fontsize(8)
p = stats.binom_test((b.beta > 0).value_counts())
print('{:,} intergenic CNVs.'.format(b.shape[0]))
print('Effect sizes for all intergenic CNV eQTLs are biased '
      '(p={:.3e}, binomial test).'.format(p))

gs.tight_layout(fig, rect=[0, 0, 1, 0.9])

t = fig.text(0.005, 0.86, 'A', weight='bold', 
             size=12)
t = fig.text(0.5, 0.86, 'B', weight='bold', 
             size=12)
t = fig.text(0.005, 0.42, 'C', weight='bold', 
             size=12)
t = fig.text(0.5, 0.42, 'D', weight='bold', 
             size=12)

fig.savefig(os.path.join(outdir, 'cnv_eqtl_effect_sizes.pdf'))
fig.savefig(os.path.join(outdir, 'cnv_eqtl_effect_sizes.png'), dpi=300)









    



102 lead genic CNVs.
Effect sizes for genic lead CNVs are biased (p=3.885e-05, binomial test).
145 lead intergenic CNVs.
Effect sizes for intergenic lead CNVs are biased (p=1.974e-02, binomial test).
370 genic CNVs.
Effect sizes for all genic CNV eQTLs are biased (p=5.330e-11, binomial test).
408 intergenic CNVs.
Effect sizes for all intergenic CNV eQTLs are biased (p=8.432e-13, binomial test).



In [5]:

    
sum(cnv_gv.gene_id.value_counts() > 1)









    Out[5]:





126