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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.047


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.033


  
Method:
             
Least Squares
  
  F-statistic:       
 
   3.374


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


  
Time:
                 
23:49:50
     
  Log-Likelihood:    
 
  884.39


  
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.7620
 
    3.468
 
    2.526
 
 0.012
 
    1.925    15.599


  
total_cpg_no_filter
 
  6.31e-09
 
 3.43e-09
 
    1.842
 
 0.067
 
-4.45e-10  1.31e-08


  
bsRate_mean
         
   -9.0390
 
    3.536
 
   -2.556
 
 0.011
 
  -16.011    -2.067


  
avgReadCpG_mean
     
    0.0140
 
    0.009
 
    1.600
 
 0.111
 
   -0.003     0.031




  
Omnibus:
       
11.066
 
  Durbin-Watson:     
 
   2.056


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


  
Skew:
          
 0.562
 
  Prob(JB):          
 
 0.00275


  
Kurtosis:
      
 2.707
 
  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 0x10bc6e940>






    

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 0x10c81d710>






    

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 0x10ca78908>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.050


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.036


  
Method:
             
Least Squares
  
  F-statistic:       
 
   3.605


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


  
Time:
                 
23:49:53
     
  Log-Likelihood:    
 
  862.29


  
No. Observations:
      
   210
      
  AIC:               
 
  -1717.


  
Df Residuals:
          
   206
      
  BIC:               
 
  -1703.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
    5.4407
 
    4.471
 
    1.217
 
 0.225
 
   -3.373    14.255


  
total_cpg_no_filter
 
   9.8e-09
 
 3.56e-09
 
    2.753
 
 0.006
 
 2.78e-09  1.68e-08


  
bsRate_mean
         
   -5.4892
 
    4.636
 
   -1.184
 
 0.238
 
  -14.630     3.652


  
avgReadCpG_mean
     
   -0.0120
 
    0.011
 
   -1.069
 
 0.286
 
   -0.034     0.010




  
Omnibus:
       
13.995
 
  Durbin-Watson:     
 
   2.461


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


  
Skew:
          
 0.406
 
  Prob(JB):          
 
 0.00744


  
Kurtosis:
      
 2.322
 
  Cond. No.          
 
1.68e+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 0x10cdc5358>






    

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 0x10cff3c88>






    

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 0x10d3014a8>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.029


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.014


  
Method:
             
Least Squares
  
  F-statistic:       
 
   1.879


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


  
Time:
                 
23:49:56
     
  Log-Likelihood:    
 
  782.07


  
No. Observations:
      
   190
      
  AIC:               
 
  -1556.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1543.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
   -3.3290
 
    1.748
 
   -1.904
 
 0.058
 
   -6.778     0.120


  
total_cpg_no_filter
 
-1.609e-09
 
 2.96e-09
 
   -0.543
 
 0.588
 
-7.45e-09  4.24e-09


  
bsRate_mean
         
    3.5459
 
    1.830
 
    1.938
 
 0.054
 
   -0.064     7.156


  
avgReadCpG_mean
     
   -0.0102
 
    0.006
 
   -1.831
 
 0.069
 
   -0.021     0.001




  
Omnibus:
       
13.706
 
  Durbin-Watson:     
 
   2.228


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


  
Skew:
          
 0.469
 
  Prob(JB):          
 
 0.00536


  
Kurtosis:
      
 2.335
 
  Cond. No.          
 
3.27e+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 0x10d06d198>






    

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 0x10d81b208>






    

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 0x10da56048>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.131


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.117


  
Method:
             
Least Squares
  
  F-statistic:       
 
   9.323


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


  
Time:
                 
23:50:00
     
  Log-Likelihood:    
 
  754.47


  
No. Observations:
      
   190
      
  AIC:               
 
  -1501.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1488.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
   -5.8954
 
    2.242
 
   -2.630
 
 0.009
 
  -10.318    -1.472


  
total_cpg_no_filter
 
-5.239e-09
 
 2.77e-09
 
   -1.891
 
 0.060
 
-1.07e-08  2.26e-10


  
bsRate_mean
         
    6.0663
 
    2.352
 
    2.579
 
 0.011
 
    1.425    10.707


  
avgReadCpG_mean
     
    0.0105
 
    0.005
 
    2.292
 
 0.023
 
    0.001     0.020




  
Omnibus:
       
16.408
 
  Durbin-Watson:     
 
   1.723


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


  
Skew:
          
 0.740
 
  Prob(JB):          
 
0.000117


  
Kurtosis:
      
 3.314
 
  Cond. No.          
 
4.47e+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 0x10dcd8ef0>






    

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 0x10dedb9e8>






    

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 0x10e120198>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.149


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.136


  
Method:
             
Least Squares
  
  F-statistic:       
 
   10.89


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


  
Time:
                 
23:50:03
     
  Log-Likelihood:    
 
  723.17


  
No. Observations:
      
   190
      
  AIC:               
 
  -1438.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1425.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -10.5371
 
    3.337
 
   -3.158
 
 0.002
 
  -17.120    -3.954


  
total_cpg_no_filter
 
-1.983e-08
 
 4.15e-09
 
   -4.782
 
 0.000
 
 -2.8e-08 -1.16e-08


  
bsRate_mean
         
   10.7370
 
    3.460
 
    3.104
 
 0.002
 
    3.912    17.562


  
avgReadCpG_mean
     
    0.0291
 
    0.011
 
    2.584
 
 0.011
 
    0.007     0.051




  
Omnibus:
       
11.003
 
  Durbin-Watson:     
 
   1.935


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


  
Skew:
          
 0.584
 
  Prob(JB):          
 
 0.00351


  
Kurtosis:
      
 3.252
 
  Cond. No.          
 
6.96e+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 0x10e489c18>






    

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 0x10e699dd8>






    

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 0x10e8c8780>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.012


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
  -0.010


  
Method:
             
Least Squares
  
  F-statistic:       
 
  0.5533


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


  
Time:
                 
23:50:06
     
  Log-Likelihood:    
 
  240.20


  
No. Observations:
      
   136
      
  AIC:               
 
  -472.4


  
Df Residuals:
          
   132
      
  BIC:               
 
  -460.7


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
   -9.8182
 
   61.115
 
   -0.161
 
 0.873
 
 -130.709   111.073


  
total_cpg_no_filter
 
-7.492e-08
 
 7.65e-08
 
   -0.979
 
 0.329
 
-2.26e-07  7.64e-08


  
bsRate_mean
         
   10.6164
 
   61.269
 
    0.173
 
 0.863
 
 -110.580   131.812


  
avgReadCpG_mean
     
   -0.0992
 
    0.079
 
   -1.249
 
 0.214
 
   -0.256     0.058




  
Omnibus:
       
10.120
 
  Durbin-Watson:     
 
   1.840


  
Prob(Omnibus):
 
 0.006
 
  Jarque-Bera (JB):  
 
  11.038


  
Skew:
          
 0.688
 
  Prob(JB):          
 
 0.00401


  
Kurtosis:
      
 2.769
 
  Cond. No.          
 
7.56e+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 0x10eb25048>






    

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 0x10ed857b8>






    

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 0x10f085f28>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.078


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.066


  
Method:
             
Least Squares
  
  F-statistic:       
 
   6.400


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


  
Time:
                 
23:50:10
     
  Log-Likelihood:    
 
  517.32


  
No. Observations:
      
   231
      
  AIC:               
 
  -1027.


  
Df Residuals:
          
   227
      
  BIC:               
 
  -1013.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -58.7101
 
   43.375
 
   -1.354
 
 0.177
 
 -144.179    26.759


  
total_cpg_no_filter
 
  7.09e-08
 
 3.33e-08
 
    2.130
 
 0.034
 
  5.3e-09  1.36e-07


  
bsRate_mean
         
   58.2610
 
   43.449
 
    1.341
 
 0.181
 
  -27.355   143.877


  
avgReadCpG_mean
     
    0.0857
 
    0.020
 
    4.249
 
 0.000
 
    0.046     0.125




  
Omnibus:
       
41.518
 
  Durbin-Watson:     
 
   1.321


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


  
Skew:
          
 1.115
 
  Prob(JB):          
 
1.93e-13


  
Kurtosis:
      
 4.054
 
  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 0x10f3bff98>






    

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 0x10f61b128>






    

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 0x10f918240>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.200


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.125


  
Method:
             
Least Squares
  
  F-statistic:       
 
   2.660


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


  
Time:
                 
23:50:13
     
  Log-Likelihood:    
 
  76.330


  
No. Observations:
      
    36
      
  AIC:               
 
  -144.7


  
Df Residuals:
          
    32
      
  BIC:               
 
  -138.3


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
 -267.9734
 
  252.270
 
   -1.062
 
 0.296
 
 -781.831   245.884


  
total_cpg_no_filter
 
-2.207e-08
 
  3.6e-08
 
   -0.614
 
 0.544
 
-9.54e-08  5.12e-08


  
bsRate_mean
         
  267.6661
 
  252.429
 
    1.060
 
 0.297
 
 -246.514   781.847


  
avgReadCpG_mean
     
    0.1455
 
    0.081
 
    1.799
 
 0.081
 
   -0.019     0.310




  
Omnibus:
       
 1.832
 
  Durbin-Watson:     
 
   1.728


  
Prob(Omnibus):
 
 0.400
 
  Jarque-Bera (JB):  
 
   1.706


  
Skew:
          
 0.450
 
  Prob(JB):          
 
   0.426


  
Kurtosis:
      
 2.428
 
  Cond. No.          
 
3.29e+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 0x10fb77a90>






    

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 0x10fd9fe48>






    

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 0x10febed30>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.3295


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


  
Time:
                 
23:50:16
     
  Log-Likelihood:    
 
  400.88


  
No. Observations:
      
   190
      
  AIC:               
 
  -793.8


  
Df Residuals:
          
   186
      
  BIC:               
 
  -780.8


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -21.1774
 
   45.736
 
   -0.463
 
 0.644
 
 -111.406    69.051


  
total_cpg_no_filter
 
  3.07e-08
 
 3.76e-08
 
    0.816
 
 0.415
 
-4.35e-08  1.05e-07


  
bsRate_mean
         
   21.1448
 
   45.725
 
    0.462
 
 0.644
 
  -69.062   111.352


  
avgReadCpG_mean
     
    0.0183
 
    0.049
 
    0.373
 
 0.709
 
   -0.078     0.115




  
Omnibus:
       
11.949
 
  Durbin-Watson:     
 
   1.590


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


  
Skew:
          
 0.627
 
  Prob(JB):          
 
 0.00149


  
Kurtosis:
      
 2.733
 
  Cond. No.          
 
8.48e+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 0x1102ab3c8>






    

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 0x1104eccc0>






    

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 0x110822588>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.179


  
Method:
             
Least Squares
  
  F-statistic:       
 
   12.01


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


  
Time:
                 
23:50:20
     
  Log-Likelihood:    
 
  418.20


  
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
               
  234.9377
 
   42.901
 
    5.476
 
 0.000
 
  150.165   319.710


  
total_cpg_no_filter
 
-5.267e-08
 
 2.36e-08
 
   -2.228
 
 0.027
 
-9.94e-08 -5.96e-09


  
bsRate_mean
         
 -235.9340
 
   43.113
 
   -5.472
 
 0.000
 
 -321.126  -150.742


  
avgReadCpG_mean
     
    0.0550
 
    0.025
 
    2.213
 
 0.028
 
    0.006     0.104




  
Omnibus:
       
50.644
 
  Durbin-Watson:     
 
   1.624


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


  
Skew:
          
 1.609
 
  Prob(JB):          
 
8.81e-21


  
Kurtosis:
      
 5.034
 
  Cond. No.          
 
1.52e+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 0x110a723c8>






    

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 0x110d3c8d0>






    

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 0x111066400>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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:       
 
   5.987


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


  
Time:
                 
23:50:24
     
  Log-Likelihood:    
 
  365.61


  
No. Observations:
      
   153
      
  AIC:               
 
  -723.2


  
Df Residuals:
          
   149
      
  BIC:               
 
  -711.1


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -24.3938
 
   68.215
 
   -0.358
 
 0.721
 
 -159.187   110.400


  
total_cpg_no_filter
 
  1.36e-07
 
 4.57e-08
 
    2.978
 
 0.003
 
 4.58e-08  2.26e-07


  
bsRate_mean
         
   24.5140
 
   68.405
 
    0.358
 
 0.721
 
 -110.655   159.683


  
avgReadCpG_mean
     
   -0.0126
 
    0.046
 
   -0.272
 
 0.786
 
   -0.104     0.079




  
Omnibus:
       
71.909
 
  Durbin-Watson:     
 
   0.216


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


  
Skew:
          
 2.083
 
  Prob(JB):          
 
1.23e-39


  
Kurtosis:
      
 6.279
 
  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 0x1113257b8>






    

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 0x1115f3d30>






    

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 0x1118aae48>






    

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_gtreql3.8CpG":"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_nofilter", "avgReadCpgs_lessthan1CpG"], 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.150


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.131


  
Method:
             
Least Squares
  
  F-statistic:       
 
   7.772


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


  
Time:
                 
23:50:28
     
  Log-Likelihood:    
 
  606.00


  
No. Observations:
      
   136
      
  AIC:               
 
  -1204.


  
Df Residuals:
          
   132
      
  BIC:               
 
  -1192.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
    2.0897
 
    7.525
 
    0.278
 
 0.782
 
  -12.796    16.975


  
total_cpg_no_filter
 
 2.636e-08
 
 6.51e-09
 
    4.049
 
 0.000
 
 1.35e-08  3.92e-08


  
bsRate_mean
         
   -2.2537
 
    7.537
 
   -0.299
 
 0.765
 
  -17.162    12.655


  
avgReadCpG_mean
     
    0.0222
 
    0.005
 
    4.544
 
 0.000
 
    0.013     0.032




  
Omnibus:
       
22.529
 
  Durbin-Watson:     
 
   1.289


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


  
Skew:
          
 0.995
 
  Prob(JB):          
 
8.24e-07


  
Kurtosis:
      
 3.991
 
  Cond. No.          
 
1.28e+10