In [1]:

    

%matplotlib inline


    

In [2]:

    

#
# Covariates are: 
#   - Number of unique CpGs per cell
#   - Median Average Read CpG per cell (or mean if normally distributed)
#   - BS rate per cell
#   - CLL or Normal status per cell
#
# For distances (i.e. PDR difference between pairs, methylation difference between pairs), 
# the covariates are the mean between the two pairs. 
#


    

In [3]:

    

import glob
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib
pd.set_option('display.max_columns', 50) # print all rows


import os
os.chdir('/Users/evanbiederstedt/Downloads/RRBS_data_files')

import statsmodels.api as sm


    

In [4]:

    

stats = pd.read_csv("all_RRBS_statistics_final.csv")


    

In [5]:

    

stats.shape


    

    
Out[5]:





(438, 15)

In [6]:

    

normal = stats[stats["type"]=="normal"]
CLL = stats[stats["type"]=="CLL"]


    

In [7]:

    

len(normal)


    

    
Out[7]:





336

In [8]:

    

len(CLL)


    

    
Out[8]:





102

In [9]:

    

mcell_cpg = pd.read_csv("mcell_avgCpGs.csv")
pcell_cpg = pd.read_csv("pcell_avgCpGs.csv")
CD19cell_cpg = pd.read_csv("CD19_avgCpGs.csv")
normalB_cell_cpg = pd.read_csv("normalB_avgCpGs.csv")
trito_cell_cpg = pd.read_csv("trito_avgCpGs.csv")
cw154_cell_cpg = pd.read_csv("cw154_cpgs.csv")


    

In [10]:

    

mcell_cpg = mcell_cpg.drop(["Unnamed: 0"], axis=1)
pcell_cpg = pcell_cpg.drop(["Unnamed: 0"], axis=1)
CD19cell_cpg = CD19cell_cpg.drop(["Unnamed: 0"], axis=1)
normalB_cell_cpg = normalB_cell_cpg.drop(["Unnamed: 0"], axis=1)
trito_cell_cpg = trito_cell_cpg.drop(["Unnamed: 0"], axis=1)
cw154_cell_cpg = cw154_cell_cpg.drop(["Unnamed: 0"], axis=1)


    

In [11]:

    

cpg_total = pd.concat([mcell_cpg, pcell_cpg, CD19cell_cpg, normalB_cell_cpg, trito_cell_cpg, cw154_cell_cpg])


    

In [12]:

    

cpg_total.shape


    

    
Out[12]:





(513, 4)

In [13]:

    

merged = stats.merge(cpg_total, on="filename")


    

In [14]:

    

merged.shape


    

    
Out[14]:





(438, 18)

In [15]:

    

merged.head()


    

    
Out[15]:






  

    

      

      
filename
      
methylation
      
PDR_total
      
methylation_unweighted
      
PDR_unweighted
      
thisMeth
      
mixedReadCount
      
total_reads
      
type
      
bio
      
protocol
      
total_cpg_no_filter
      
total_cpg_gtrthan1
      
total_cpg_gtrthan38
      
bsRate
      
avgReadCpgs_nofilter
      
avgReadCpgs_lessthan1CpG
      
avgReadCpgs_gtreql3.8CpG
    
  
  

    

      
0
      
RRBS_NormalBCD19pCD27mcell67_88_CGTACTAG.CATGAC
      
0.529505
      
0.235795
      
0.632802
      
0.231878
      
2208325.0
      
983394.0
      
4170549.0
      
normal
      
CD19CD27m
      
NormalBCD19pCD27mcell67_88
      
525282.0
      
525251.0
      
435636.0
      
0.9975
      
5.354284
      
5.355660
      
7.019255
    
    

      
1
      
RRBS_NormalBCD19pCD27mcell67_88_CGTACTAG.CCTTCG
      
0.455550
      
0.177631
      
0.583859
      
0.175371
      
733064.0
      
285841.0
      
1609185.0
      
normal
      
CD19CD27m
      
NormalBCD19pCD27mcell67_88
      
221972.0
      
221962.0
      
186757.0
      
0.9975
      
5.587294
      
5.588449
      
7.302612
    
    

      
2
      
RRBS_NormalBCD19pCD27mcell67_88_CGTACTAG.CGGTAG
      
0.515269
      
0.177645
      
0.618578
      
0.174221
      
1452802.0
      
500870.0
      
2819500.0
      
normal
      
CD19CD27m
      
NormalBCD19pCD27mcell67_88
      
355730.0
      
355713.0
      
295624.0
      
0.9975
      
5.393199
      
5.394331
      
7.079288
    
    

      
3
      
RRBS_NormalBCD19pCD27mcell67_88_CGTACTAG.CTATTG
      
0.556175
      
0.176367
      
0.652727
      
0.172273
      
2279354.0
      
722800.0
      
4098270.0
      
normal
      
CD19CD27m
      
NormalBCD19pCD27mcell67_88
      
483179.0
      
483150.0
      
397812.0
      
0.9975
      
5.287116
      
5.288477
      
6.979525
    
    

      
4
      
RRBS_NormalBCD19pCD27mcell67_88_CGTACTAG.CTCAGC
      
0.528642
      
0.181331
      
0.640401
      
0.172287
      
1394208.0
      
478231.0
      
2637340.0
      
normal
      
CD19CD27m
      
NormalBCD19pCD27mcell67_88
      
356122.0
      
356100.0
      
294065.0
      
0.9975
      
5.314302
      
5.315719
      
6.995052
    
  

In [ ]:

    




    

In [ ]:

    




    

In [16]:

    

tritopool = merged[merged["protocol"] == 'trito_pool_1']       # select only "trito_pool_1" files
print(len(tritopool))
tritopool = tritopool.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
tritopoolA = tritopool.set_index("filename")
from itertools import combinations
cc = list(combinations(tritopool.filename,2)) # combines into all pairs
out = pd.DataFrame([tritopoolA.loc[c,:].mean() for c in cc], index=cc)  # covariates between pairs == mean
df_ex = pd.DataFrame(np.abs(np.subtract.outer(tritopool.PDR_unweighted, tritopool.PDR_unweighted)), tritopool.filename, tritopool.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs1 = pd.merge(out, PDR_differences, how='inner')
print(pairs1.shape)

pairs1 = pairs1.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs1.PDR_difference # dependent variable to predict

X = pairs1.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for CLL 'RRBS_trito_pool1', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




21
(210, 14)
(210,)
(210, 3)
Regression results for CLL 'RRBS_trito_pool1', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[16]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.046


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.033


  
Method:
             
Least Squares
  
  F-statistic:       
 
   3.346


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
0.0201
 


  
Time:
                 
21:41:58
     
  Log-Likelihood:    
 
  884.34


  
No. Observations:
      
   210
      
  AIC:               
 
  -1761.


  
Df Residuals:
          
   206
      
  BIC:               
 
  -1747.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
    8.9491
 
    3.463
 
    2.584
 
 0.010
 
    2.121    15.777


  
total_cpg_no_filter
 
 4.547e-09
 
 3.87e-09
 
    1.175
 
 0.241
 
-3.08e-09  1.22e-08


  
bsRate_mean
         
   -9.2084
 
    3.533
 
   -2.606
 
 0.010
 
  -16.175    -2.242


  
avgReadCpG_mean
     
    0.0147
 
    0.009
 
    1.574
 
 0.117
 
   -0.004     0.033




  
Omnibus:
       
11.112
 
  Durbin-Watson:     
 
   2.070


  
Prob(Omnibus):
 
 0.004
 
  Jarque-Bera (JB):  
 
  11.804


  
Skew:
          
 0.560
 
  Prob(JB):          
 
 0.00273


  
Kurtosis:
      
 2.692
 
  Cond. No.          
 
1.46e+10

In [17]:

    

print("PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_1")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs1, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_1")


    

    




PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_1






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[17]:





<matplotlib.text.Text at 0x10b8a0eb8>






    

In [18]:

    

print("PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_1")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs1, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_1")


    

    




PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_1






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[18]:





<matplotlib.text.Text at 0x10bc472b0>






    

In [19]:

    

print("PDR_difference vs. total unique CpGs, jointplot, RRBS_trito_pool CLL trito_1")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs1, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_1")


    

    




PDR_difference vs. total unique CpGs, jointplot, RRBS_trito_pool CLL trito_1






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[19]:





<matplotlib.text.Text at 0x10bdcab70>






    

In [20]:

    

tritopool = merged[merged["protocol"] == 'trito_pool_2']       # select only "trito_pool_1" files
print(len(tritopool))
tritopool = tritopool.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
tritopoolA = tritopool.set_index("filename")
from itertools import combinations
cc = list(combinations(tritopool.filename,2)) # combines into all pairs
out = pd.DataFrame([tritopoolA.loc[c,:].mean() for c in cc], index=cc)  # covariates between pairs == mean
df_ex = pd.DataFrame(np.abs(np.subtract.outer(tritopool.PDR_unweighted, tritopool.PDR_unweighted)), tritopool.filename, tritopool.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs1a = pd.merge(out, PDR_differences, how='inner')
print(pairs1a.shape)
pairs1a = pairs1a.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs1a.PDR_difference # dependent variable to predict

X = pairs1a.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for CLL 'RRBS_trito_pool1', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




21
(210, 14)
(210,)
(210, 3)
Regression results for CLL 'RRBS_trito_pool1', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[20]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.045


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.031


  
Method:
             
Least Squares
  
  F-statistic:       
 
   3.218


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
0.0238
 


  
Time:
                 
21:42:01
     
  Log-Likelihood:    
 
  861.73


  
No. Observations:
      
   210
      
  AIC:               
 
  -1715.


  
Df Residuals:
          
   206
      
  BIC:               
 
  -1702.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
    7.7698
 
    5.477
 
    1.419
 
 0.158
 
   -3.028    18.568


  
total_cpg_no_filter
 
 1.016e-08
 
 3.65e-09
 
    2.786
 
 0.006
 
 2.97e-09  1.73e-08


  
bsRate_mean
         
   -7.9514
 
    5.670
 
   -1.402
 
 0.162
 
  -19.129     3.226


  
avgReadCpG_mean
     
   -0.0023
 
    0.013
 
   -0.180
 
 0.858
 
   -0.027     0.023




  
Omnibus:
       
14.102
 
  Durbin-Watson:     
 
   2.436


  
Prob(Omnibus):
 
 0.001
 
  Jarque-Bera (JB):  
 
  10.064


  
Skew:
          
 0.417
 
  Prob(JB):          
 
 0.00653


  
Kurtosis:
      
 2.327
 
  Cond. No.          
 
2.05e+10

In [21]:

    

print("PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_2")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs1a, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_2")


    

    




PDR_difference vs. bsRate, jointplot, RRBS_trito_pool CLL trito_2






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[21]:





<matplotlib.text.Text at 0x10ca304a8>






    

In [22]:

    

print("PDR_difference vs. mean avgReadCpGs, jointplot, RRBS_trito_pool CLL trito_2")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs1a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_2")


    

    




PDR_difference vs. mean avgReadCpGs, jointplot, RRBS_trito_pool CLL trito_2






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[22]:





<matplotlib.text.Text at 0x10cc46128>






    

In [23]:

    

print("PDR_difference vs. total unique CpGs, jointplot, RRBS_trito_pool CLL trito_2")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs1a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, RRBS_trito_pool CLL trito_2")


    

    




PDR_difference vs. total unique CpGs, jointplot, RRBS_trito_pool CLL trito_2






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[23]:





<matplotlib.text.Text at 0x10cf51e80>






    

In [ ]:

    




    

In [24]:

    

cw154 = merged[merged["protocol"] == 'cw154_Tris_protease']
print(len(cw154))
cw154 = cw154.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
cw154 = cw154.reset_index(drop=True)
cw154A = cw154.set_index("filename")
from itertools import combinations
cc = list(combinations(cw154.filename,2))
out = pd.DataFrame([cw154A.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(cw154.PDR_unweighted, cw154.PDR_unweighted)), cw154.filename, cw154.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs2 = pd.merge(out, PDR_differences, how='inner')
print(pairs2.shape)
pairs2 = pairs2.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs2.PDR_difference # dependent variable to predict

X = pairs2.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for CLL 'cw154_Tris_protease', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




20
(190, 14)
(190,)
(190, 3)
Regression results for CLL 'cw154_Tris_protease', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[24]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.014


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
  -0.002


  
Method:
             
Least Squares
  
  F-statistic:       
 
  0.8881


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
 0.448
 


  
Time:
                 
21:42:05
     
  Log-Likelihood:    
 
  780.58


  
No. Observations:
      
   190
      
  AIC:               
 
  -1553.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1540.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
   -2.8995
 
    1.809
 
   -1.603
 
 0.111
 
   -6.468     0.669


  
total_cpg_no_filter
 
-1.407e-09
 
 2.99e-09
 
   -0.470
 
 0.639
 
-7.31e-09   4.5e-09


  
bsRate_mean
         
    3.0408
 
    1.891
 
    1.608
 
 0.110
 
   -0.691     6.772


  
avgReadCpG_mean
     
   -0.0030
 
    0.005
 
   -0.643
 
 0.521
 
   -0.012     0.006




  
Omnibus:
       
13.402
 
  Durbin-Watson:     
 
   2.209


  
Prob(Omnibus):
 
 0.001
 
  Jarque-Bera (JB):  
 
  11.455


  
Skew:
          
 0.518
 
  Prob(JB):          
 
 0.00326


  
Kurtosis:
      
 2.390
 
  Cond. No.          
 
3.35e+09

In [25]:

    

print("PDR_difference vs. bsRate, jointplot,  CLL cw154_Tris_protease")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs2, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot, cw154  CLL cw154_Tris_protease")


    

    




PDR_difference vs. bsRate, jointplot,  CLL cw154_Tris_protease






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[25]:





<matplotlib.text.Text at 0x10ca56e80>






    

In [26]:

    

print("PDR_difference vs. avgReadCpGs, jointplot, CLL cw154_Tris_protease")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs2, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154  CLL cw154_Tris_protease")


    

    




PDR_difference vs. avgReadCpGs, jointplot, CLL cw154_Tris_protease






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[26]:





<matplotlib.text.Text at 0x10d468f28>






    

In [27]:

    

print("PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_Tris_protease")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs2, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154 CLL cw154_Tris_protease")


    

    




PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_Tris_protease






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[27]:





<matplotlib.text.Text at 0x10d69e9e8>






    

In [28]:

    

cw154 = merged[merged["protocol"] == 'cw154_Tris_protease_GR']
print(len(cw154))
cw154 = cw154.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
cw154 = cw154.reset_index(drop=True)
cw154A = cw154.set_index("filename")
from itertools import combinations
cc = list(combinations(cw154.filename,2))
out = pd.DataFrame([cw154A.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(cw154.PDR_unweighted, cw154.PDR_unweighted)), cw154.filename, cw154.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs2a = pd.merge(out, PDR_differences, how='inner')
print(pairs2a.shape)
pairs2a = pairs2a.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs2a.PDR_difference # dependent variable to predict

X = pairs2a.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for CLL 'cw154_Tris_protease_GR', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




20
(190, 14)
(190,)
(190, 3)
Regression results for CLL 'cw154_Tris_protease_GR', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[28]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.134


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.120


  
Method:
             
Least Squares
  
  F-statistic:       
 
   9.604


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
6.30e-06


  
Time:
                 
21:42:08
     
  Log-Likelihood:    
 
  754.84


  
No. Observations:
      
   190
      
  AIC:               
 
  -1502.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1489.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
   -5.1998
 
    2.359
 
   -2.204
 
 0.029
 
   -9.853    -0.546


  
total_cpg_no_filter
 
-6.633e-09
 
 2.94e-09
 
   -2.254
 
 0.025
 
-1.24e-08 -8.27e-10


  
bsRate_mean
         
    5.3602
 
    2.470
 
    2.170
 
 0.031
 
    0.487    10.233


  
avgReadCpG_mean
     
    0.0108
 
    0.004
 
    2.451
 
 0.015
 
    0.002     0.019




  
Omnibus:
       
16.926
 
  Durbin-Watson:     
 
   1.740


  
Prob(Omnibus):
 
 0.000
 
  Jarque-Bera (JB):  
 
  18.834


  
Skew:
          
 0.756
 
  Prob(JB):          
 
8.13e-05


  
Kurtosis:
      
 3.304
 
  Cond. No.          
 
4.71e+09

In [29]:

    

print("PDR_difference vs. bsRate, jointplot,  CLL cw154_Tris_protease_GR")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs2a, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot, cw154  CLL cw154_Tris_protease_GR")


    

    




PDR_difference vs. bsRate, jointplot,  CLL cw154_Tris_protease_GR






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[29]:





<matplotlib.text.Text at 0x10d91db00>






    

In [30]:

    

print("PDR_difference vs. avgReadCpgs, jointplot, CLL cw154_Tris_protease_GR")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs2a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154  CLL cw154_Tris_protease_GR")


    

    




PDR_difference vs. avgReadCpgs, jointplot, CLL cw154_Tris_protease_GR






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[30]:





<matplotlib.text.Text at 0x10db24860>






    

In [31]:

    

print("PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_Tris_protease_GR")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs2a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154 CLL cw154_Tris_protease_GR")


    

    




PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_Tris_protease_GR






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[31]:





<matplotlib.text.Text at 0x10de36a90>






    

In [ ]:

    




    

In [ ]:

    




    

In [32]:

    

cw154 = merged[merged["protocol"] == 'cw154_CutSmart_proteinase_K']
print(len(cw154))
cw154 = cw154.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
cw154 = cw154.reset_index(drop=True)
cw154A = cw154.set_index("filename")
from itertools import combinations
cc = list(combinations(cw154.filename,2))
out = pd.DataFrame([cw154A.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(cw154.PDR_unweighted, cw154.PDR_unweighted)), cw154.filename, cw154.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs2b = pd.merge(out, PDR_differences, how='inner')
print(pairs2b.shape)
pairs2b = pairs2b.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs2b.PDR_difference # dependent variable to predict

X = pairs2b.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for CLL 'cw154_CutSmart_proteinase_K', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




20
(190, 14)
(190,)
(190, 3)
Regression results for CLL 'cw154_CutSmart_proteinase_K', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[32]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.131


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.117


  
Method:
             
Least Squares
  
  F-statistic:       
 
   9.334


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
8.86e-06


  
Time:
                 
21:42:12
     
  Log-Likelihood:    
 
  721.12


  
No. Observations:
      
   190
      
  AIC:               
 
  -1434.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1421.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
   -9.8580
 
    3.375
 
   -2.921
 
 0.004
 
  -16.516    -3.200


  
total_cpg_no_filter
 
-2.227e-08
 
 4.24e-09
 
   -5.251
 
 0.000
 
-3.06e-08 -1.39e-08


  
bsRate_mean
         
   10.1443
 
    3.506
 
    2.894
 
 0.004
 
    3.229    17.060


  
avgReadCpG_mean
     
    0.0183
 
    0.011
 
    1.598
 
 0.112
 
   -0.004     0.041




  
Omnibus:
       
13.999
 
  Durbin-Watson:     
 
   1.928


  
Prob(Omnibus):
 
 0.001
 
  Jarque-Bera (JB):  
 
  14.870


  
Skew:
          
 0.656
 
  Prob(JB):          
 
0.000590


  
Kurtosis:
      
 3.398
 
  Cond. No.          
 
6.97e+09

In [33]:

    

print("PDR_difference vs. bsRate, jointplot,  CLL cw154_CutSmart_proteinase_K")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs2b, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot, cw154 CLL cw154_CutSmart_proteinase_K")


    

    




PDR_difference vs. bsRate, jointplot,  CLL cw154_CutSmart_proteinase_K






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[33]:





<matplotlib.text.Text at 0x10e0cd5c0>






    

In [34]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot, CLL cw154_CutSmart_proteinase_K")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs2b, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154  CLL cw154_CutSmart_proteinase_K")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot, CLL cw154_CutSmart_proteinase_K






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[34]:





<matplotlib.text.Text at 0x10e2e5518>






    

In [35]:

    

print("PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_CutSmart_proteinase_K")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs2b, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, cw154 CLL cw154_CutSmart_proteinase_K")


    

    




PDR_difference vs. total unique CpGs, jointplot, cw154 CLL cw154_CutSmart_proteinase_K






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[35]:





<matplotlib.text.Text at 0x10e519198>






    

In [ ]:

    




    

In [36]:

    

pcell = merged[merged["protocol"] == 'NormalBCD19pCD27pcell1_22_']
print(len(pcell))
pcell = pcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
pcell = pcell.reset_index(drop=True)
pcellA = pcell.set_index("filename")
from itertools import combinations
cc = list(combinations(pcell.filename,2))
out = pd.DataFrame([pcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(pcell.PDR_unweighted, pcell.PDR_unweighted)), pcell.filename, pcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs3 = pd.merge(out, PDR_differences, how='inner')
print(pairs3.shape)
pairs3 = pairs3.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs3.PDR_difference # dependent variable to predict

X = pairs3.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




17
(136, 14)
(136,)
(136, 3)
Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[36]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.019


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
  -0.004


  
Method:
             
Least Squares
  
  F-statistic:       
 
  0.8348


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
 0.477
 


  
Time:
                 
21:42:15
     
  Log-Likelihood:    
 
  240.63


  
No. Observations:
      
   136
      
  AIC:               
 
  -473.3


  
Df Residuals:
          
   132
      
  BIC:               
 
  -461.6


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
    9.8844
 
   62.393
 
    0.158
 
 0.874
 
 -113.535   133.304


  
total_cpg_no_filter
 
-9.493e-08
 
 7.87e-08
 
   -1.207
 
 0.230
 
-2.51e-07  6.07e-08


  
bsRate_mean
         
   -9.0554
 
   62.449
 
   -0.145
 
 0.885
 
 -132.585   114.475


  
avgReadCpG_mean
     
   -0.1446
 
    0.093
 
   -1.550
 
 0.123
 
   -0.329     0.040




  
Omnibus:
       
 9.622
 
  Durbin-Watson:     
 
   1.838


  
Prob(Omnibus):
 
 0.008
 
  Jarque-Bera (JB):  
 
  10.406


  
Skew:
          
 0.665
 
  Prob(JB):          
 
 0.00550


  
Kurtosis:
      
 2.740
 
  Cond. No.          
 
7.73e+09

In [37]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs3, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[37]:





<matplotlib.text.Text at 0x10e767898>






    

In [38]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs3, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27pcell1_22_")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[38]:





<matplotlib.text.Text at 0x10e9d69b0>






    

In [39]:

    

print("PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs3, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_")


    

    




PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[39]:





<matplotlib.text.Text at 0x10ecdd940>






    

In [ ]:

    




    

In [40]:

    

pcell = merged[merged["protocol"] == 'NormalBCD19pCD27pcell23_44']
print(len(pcell))
pcell = pcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
pcell = pcell.reset_index(drop=True)
pcellA = pcell.set_index("filename")
from itertools import combinations
cc = list(combinations(pcell.filename,2))
out = pd.DataFrame([pcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(pcell.PDR_unweighted, pcell.PDR_unweighted)), pcell.filename, pcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs3a = pd.merge(out, PDR_differences, how='inner')
print(pairs3a.shape)
pairs3a = pairs3a.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs3a.PDR_difference # dependent variable to predict

X = pairs3a.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




22
(231, 14)
(231,)
(231, 3)
Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[40]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.086


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.074


  
Method:
             
Least Squares
  
  F-statistic:       
 
   7.148


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
0.000131


  
Time:
                 
21:42:18
     
  Log-Likelihood:    
 
  518.36


  
No. Observations:
      
   231
      
  AIC:               
 
  -1029.


  
Df Residuals:
          
   227
      
  BIC:               
 
  -1015.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
  -53.4329
 
   42.873
 
   -1.246
 
 0.214
 
 -137.913    31.047


  
total_cpg_no_filter
 
 7.974e-08
 
 3.37e-08
 
    2.367
 
 0.019
 
 1.34e-08  1.46e-07


  
bsRate_mean
         
   53.0258
 
   42.955
 
    1.234
 
 0.218
 
  -31.616   137.667


  
avgReadCpG_mean
     
    0.1018
 
    0.023
 
    4.504
 
 0.000
 
    0.057     0.146




  
Omnibus:
       
41.269
 
  Durbin-Watson:     
 
   1.288


  
Prob(Omnibus):
 
 0.000
 
  Jarque-Bera (JB):  
 
  58.104


  
Skew:
          
 1.109
 
  Prob(JB):          
 
2.41e-13


  
Kurtosis:
      
 4.058
 
  Cond. No.          
 
1.03e+10

In [41]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell23_44")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs3a, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell23_44")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[41]:





<matplotlib.text.Text at 0x10f037518>






    

In [42]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell23_44")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs3a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27pcell23_44")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[42]:





<matplotlib.text.Text at 0x10f271128>






    

In [43]:

    

print("PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell23_44")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs3a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27pcell23_44")


    

    




PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[43]:





<matplotlib.text.Text at 0x10f4a5f98>






    

In [ ]:

    




    

In [44]:

    

pcell = merged[merged["protocol"] == 'NormalBCD19pCD27pcell45_66']
print(len(pcell))
pcell = pcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
pcell = pcell.reset_index(drop=True)
pcellA = pcell.set_index("filename")
from itertools import combinations
cc = list(combinations(pcell.filename,2))
out = pd.DataFrame([pcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(pcell.PDR_unweighted, pcell.PDR_unweighted)), pcell.filename, pcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs3b = pd.merge(out, PDR_differences, how='inner')
print(pairs3b.shape)
pairs3b = pairs3b.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs3b.PDR_difference # dependent variable to predict

X = pairs3b.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




9
(36, 14)
(36,)
(36, 3)
Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[44]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.170


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.092


  
Method:
             
Least Squares
  
  F-statistic:       
 
   2.185


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
 0.109
 


  
Time:
                 
21:42:21
     
  Log-Likelihood:    
 
  75.676


  
No. Observations:
      
    36
      
  AIC:               
 
  -143.4


  
Df Residuals:
          
    32
      
  BIC:               
 
  -137.0


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
 -193.0638
 
  256.757
 
   -0.752
 
 0.458
 
 -716.061   329.933


  
total_cpg_no_filter
 
-2.863e-08
 
 3.73e-08
 
   -0.768
 
 0.448
 
-1.05e-07  4.73e-08


  
bsRate_mean
         
  192.8895
 
  256.980
 
    0.751
 
 0.458
 
 -330.561   716.340


  
avgReadCpG_mean
     
    0.1321
 
    0.094
 
    1.407
 
 0.169
 
   -0.059     0.323




  
Omnibus:
       
 1.669
 
  Durbin-Watson:     
 
   1.709


  
Prob(Omnibus):
 
 0.434
 
  Jarque-Bera (JB):  
 
   1.582


  
Skew:
          
 0.442
 
  Prob(JB):          
 
   0.453


  
Kurtosis:
      
 2.478
 
  Cond. No.          
 
3.28e+10

In [45]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell45_66")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs3b, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell45_66")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell45_66






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[45]:





<matplotlib.text.Text at 0x10f7d7438>






    

In [46]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell45_66")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs3b, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27pcell45_66")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27pcell45_66






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[46]:





<matplotlib.text.Text at 0x10f9ff7f0>






    

In [47]:

    

print("PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell45_66")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs3b, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27pcell45_66")


    

    




PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27pcell45_66






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[47]:





<matplotlib.text.Text at 0x10fbe4588>






    

In [ ]:

    




    

In [ ]:

    




    

In [48]:

    

pcell = merged[merged["protocol"] == 'NormalBCD19pCD27pcell67_88']
print(len(pcell))
pcell = pcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
pcell = pcell.reset_index(drop=True)
pcellA = pcell.set_index("filename")
from itertools import combinations
cc = list(combinations(pcell.filename,2))
out = pd.DataFrame([pcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(pcell.PDR_unweighted, pcell.PDR_unweighted)), pcell.filename, pcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs3c = pd.merge(out, PDR_differences, how='inner')
print(pairs3c.shape)
pairs3c = pairs3c.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs3c.PDR_difference # dependent variable to predict

X = pairs3c.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




20
(190, 14)
(190,)
(190, 3)
Regression results for Normal 'NormalBCD19pCD27pcell1_22_', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[48]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.005


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
  -0.011


  
Method:
             
Least Squares
  
  F-statistic:       
 
  0.2850


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
  
 0.836
 


  
Time:
                 
21:42:25
     
  Log-Likelihood:    
 
  400.81


  
No. Observations:
      
   190
      
  AIC:               
 
  -793.6


  
Df Residuals:
          
   186
      
  BIC:               
 
  -780.6


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
  -15.7505
 
   46.752
 
   -0.337
 
 0.737
 
 -107.984    76.483


  
total_cpg_no_filter
 
 2.689e-08
 
 3.69e-08
 
    0.728
 
 0.468
 
 -4.6e-08  9.98e-08


  
bsRate_mean
         
   15.8050
 
   46.744
 
    0.338
 
 0.736
 
  -76.412   108.022


  
avgReadCpG_mean
     
    0.0048
 
    0.059
 
    0.081
 
 0.936
 
   -0.112     0.121




  
Omnibus:
       
11.921
 
  Durbin-Watson:     
 
   1.586


  
Prob(Omnibus):
 
 0.003
 
  Jarque-Bera (JB):  
 
  13.004


  
Skew:
          
 0.628
 
  Prob(JB):          
 
 0.00150


  
Kurtosis:
      
 2.741
 
  Cond. No.          
 
8.67e+09

In [49]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell67_88")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs3c, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell67_88")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27pcell67_88






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[49]:





<matplotlib.text.Text at 0x10ff16c18>






    

In [50]:

    

print("PDR_difference vs. avgReadCpGs, jointplot,  Normal 'NormalBCD19pCD27pcell67_88")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs3c, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27pcell67_88")


    

    




PDR_difference vs. avgReadCpGs, jointplot,  Normal 'NormalBCD19pCD27pcell67_88






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[50]:





<matplotlib.text.Text at 0x11014f748>






    

In [51]:

    

print("PDR_difference vs. total unique CpGs per cell, jointplot,  Normal 'NormalBCD19pCD27pcell67_88")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs3c, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27pcell67_88")


    

    




PDR_difference vs. total unique CpGs per cell, jointplot,  Normal 'NormalBCD19pCD27pcell67_88






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[51]:





<matplotlib.text.Text at 0x1104879e8>






    

In [ ]:

    




    

In [ ]:

    




    

In [52]:

    

mcell = merged[merged["protocol"] == 'NormalBCD19pCD27mcell1_22_']
print(len(mcell))
mcell = mcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
mcell = mcell.reset_index(drop=True)
mcellA = mcell.set_index("filename")
from itertools import combinations
cc = list(combinations(mcell.filename,2))
out = pd.DataFrame([mcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(mcell.PDR_unweighted, mcell.PDR_unweighted)), mcell.filename, mcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs4 = pd.merge(out, PDR_differences, how='inner')
print(pairs4.shape)
pairs4 = pairs4.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs4.PDR_difference # dependent variable to predict

X = pairs4.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27mcell1_22_', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




18
(153, 14)
(153,)
(153, 3)
Regression results for Normal 'NormalBCD19pCD27mcell1_22_', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[52]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.195


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.178


  
Method:
             
Least Squares
  
  F-statistic:       
 
   12.00


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
4.41e-07


  
Time:
                 
21:42:28
     
  Log-Likelihood:    
 
  418.19


  
No. Observations:
      
   153
      
  AIC:               
 
  -828.4


  
Df Residuals:
          
   149
      
  BIC:               
 
  -816.3


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
  221.6233
 
   39.953
 
    5.547
 
 0.000
 
  142.676   300.571


  
total_cpg_no_filter
 
-4.771e-08
 
 2.51e-08
 
   -1.900
 
 0.059
 
-9.73e-08  1.92e-09


  
bsRate_mean
         
 -222.5182
 
   40.124
 
   -5.546
 
 0.000
 
 -301.805  -143.232


  
avgReadCpG_mean
     
    0.0595
 
    0.027
 
    2.207
 
 0.029
 
    0.006     0.113




  
Omnibus:
       
50.762
 
  Durbin-Watson:     
 
   1.626


  
Prob(Omnibus):
 
 0.000
 
  Jarque-Bera (JB):  
 
  92.796


  
Skew:
          
 1.610
 
  Prob(JB):          
 
7.07e-21


  
Kurtosis:
      
 5.045
 
  Cond. No.          
 
1.42e+10

In [53]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs4, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[53]:





<matplotlib.text.Text at 0x1106d2710>






    

In [54]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs4, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27mcell1_22_")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[54]:





<matplotlib.text.Text at 0x1109a1b00>






    

In [55]:

    

print("PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs4, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_")


    

    




PDR_difference vs. total unique CpGs, jointplot,  Normal 'NormalBCD19pCD27mcell1_22_






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[55]:





<matplotlib.text.Text at 0x110cc6400>






    

In [ ]:

    




    

In [56]:

    

mcell = merged[merged["protocol"] == 'NormalBCD19pCD27mcell23_44']
print(len(mcell))
mcell = mcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
mcell = mcell.reset_index(drop=True)
mcellA = mcell.set_index("filename")
from itertools import combinations
cc = list(combinations(mcell.filename,2))
out = pd.DataFrame([mcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(mcell.PDR_unweighted, mcell.PDR_unweighted)), mcell.filename, mcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs4a = pd.merge(out, PDR_differences, how='inner')
print(pairs4a.shape)
pairs4a = pairs4a.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs4a.PDR_difference # dependent variable to predict

X = pairs4a.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27mcell23_44', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




18
(153, 14)
(153,)
(153, 3)
Regression results for Normal 'NormalBCD19pCD27mcell23_44', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[56]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.108


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.090


  
Method:
             
Least Squares
  
  F-statistic:       
 
   6.029


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
0.000665


  
Time:
                 
21:42:32
     
  Log-Likelihood:    
 
  365.67


  
No. Observations:
      
   153
      
  AIC:               
 
  -723.3


  
Df Residuals:
          
   149
      
  BIC:               
 
  -711.2


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
  -24.0830
 
   68.026
 
   -0.354
 
 0.724
 
 -158.504   110.338


  
total_cpg_no_filter
 
  1.29e-07
 
 4.91e-08
 
    2.630
 
 0.009
 
 3.21e-08  2.26e-07


  
bsRate_mean
         
   24.2419
 
   68.201
 
    0.355
 
 0.723
 
 -110.524   159.008


  
avgReadCpG_mean
     
   -0.0234
 
    0.054
 
   -0.431
 
 0.667
 
   -0.131     0.084




  
Omnibus:
       
71.754
 
  Durbin-Watson:     
 
   0.220


  
Prob(Omnibus):
 
 0.000
 
  Jarque-Bera (JB):  
 
 178.319


  
Skew:
          
 2.080
 
  Prob(JB):          
 
1.90e-39


  
Kurtosis:
      
 6.266
 
  Cond. No.          
 
1.55e+10

In [57]:

    

print("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell23_44")
sns.jointplot(x="bsRate_mean", y="PDR_difference",  data=pairs4a, kind="reg")
plt.title("PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell23_44")


    

    




PDR_difference vs. bsRate, jointplot,  Normal 'NormalBCD19pCD27mcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[57]:





<matplotlib.text.Text at 0x11103a630>






    

In [58]:

    

print("PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27mcell23_44")
sns.jointplot(x="avgReadCpG_mean", y="PDR_difference",  data=pairs4a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot, Normal 'NormalBCD19pCD27mcell23_44")


    

    




PDR_difference vs. avgReadCpG_mean, jointplot,  Normal 'NormalBCD19pCD27mcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[58]:





<matplotlib.text.Text at 0x111325780>






    

In [59]:

    

print("PDR_difference vs. total unique CpGs per cell, jointplot,  Normal 'NormalBCD19pCD27mcell23_44")
sns.jointplot(x="total_cpg_no_filter", y="PDR_difference",  data=pairs4a, kind="reg")
plt.title("PDR_difference vs. avgReadCpGs_mean, jointplot,  Normal 'NormalBCD19pCD27mcell23_44")


    

    




PDR_difference vs. total unique CpGs per cell, jointplot,  Normal 'NormalBCD19pCD27mcell23_44






    




/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[59]:





<matplotlib.text.Text at 0x11169e550>






    

In [ ]:

    




    

In [60]:

    

mcell = merged[merged["protocol"] == 'NormalBCD19pCD27mcell45_66']
print(len(mcell))
mcell = mcell.drop(["thisMeth", "mixedReadCount", "bio", "protocol"], axis=1) # will not need these columns; efficiency
mcell = mcell.reset_index(drop=True)
mcellA = mcell.set_index("filename")
from itertools import combinations
cc = list(combinations(mcell.filename,2))
out = pd.DataFrame([mcellA.loc[c,:].mean() for c in cc], index=cc)
df_ex = pd.DataFrame(np.abs(np.subtract.outer(mcell.PDR_unweighted, mcell.PDR_unweighted)), mcell.filename, mcell.filename)
stacked = df_ex.stack()
PDR_differences = pd.DataFrame({'filename': stacked.index.to_series(), 'PDR_difference': stacked})[['filename', 'PDR_difference']].reset_index(drop=True)
out['filename'] = out.index
out = out.reset_index(drop=True)
pairs4b = pd.merge(out, PDR_differences, how='inner')
print(pairs4b.shape)
pairs4b = pairs4b.rename(columns = {'total_reads':'total_reads_mean', "bsRate":"bsRate_mean", "avgReadCpgs_nofilter":"avgReadCpG_mean"})

y = pairs4b.PDR_difference # dependent variable to predict

X = pairs4b.drop(["PDR_difference", "filename", "methylation", "PDR_total", "methylation_unweighted", "PDR_unweighted", "total_reads_mean",
               "total_cpg_gtrthan1", "total_cpg_gtrthan38", "avgReadCpgs_lessthan1CpG", "avgReadCpgs_gtreql3.8CpG"], axis=1)

print(y.shape)
print(X.shape)
X = sm.add_constant(X)
est = sm.OLS(y, X).fit()
print("Regression results for Normal 'NormalBCD19pCD27mcell45_66', predict \delta PDR_unweighted")
print("")
print("Variates used are " + str(X.columns))
print("")
est.summary()


    

    




17
(136, 14)
(136,)
(136, 3)
Regression results for Normal 'NormalBCD19pCD27mcell45_66', predict \delta PDR_unweighted

Variates used are Index(['const', 'total_cpg_no_filter', 'bsRate_mean', 'avgReadCpG_mean'], dtype='object')







    
Out[60]:





OLS Regression Results

  
Dep. Variable:
     
PDR_difference
  
  R-squared:         
 
   0.174


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.155


  
Method:
             
Least Squares
  
  F-statistic:       
 
   9.257


  
Date:
             
Tue, 09 Aug 2016
 
  Prob (F-statistic):
 
1.34e-05


  
Time:
                 
21:42:36
     
  Log-Likelihood:    
 
  607.92


  
No. Observations:
      
   136
      
  AIC:               
 
  -1208.


  
Df Residuals:
          
   132
      
  BIC:               
 
  -1196.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




           
              
coef
     
std err
      
t
      
P>|t|
 
[95.0% Conf. Int.]
 


  
const
               
    2.2489
 
    7.386
 
    0.304
 
 0.761
 
  -12.361    16.858


  
total_cpg_no_filter
 
 2.618e-08
 
 6.08e-09
 
    4.304
 
 0.000
 
 1.41e-08  3.82e-08


  
bsRate_mean
         
   -2.3824
 
    7.400
 
   -0.322
 
 0.748
 
  -17.020    12.255


  
avgReadCpG_mean
     
    0.0233
 
    0.005
 
    5.002
 
 0.000
 
    0.014     0.033




  
Omnibus:
       
21.471
 
  Durbin-Watson:     
 
   1.299


  
Prob(Omnibus):
 
 0.000
 
  Jarque-Bera (JB):  
 
  26.182


  
Skew:
          
 0.971
 
  Prob(JB):          
 
2.06e-06


  
Kurtosis:
      
 3.920
 
  Cond. No.          
 
1.27e+10