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_total, tritopool.PDR_total)), 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_total")
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_total

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.015


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.001


  
Method:
             
Least Squares
  
  F-statistic:       
 
   1.038


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
 0.377
 


  
Time:
                 
00:35:47
     
  Log-Likelihood:    
 
  826.52


  
No. Observations:
      
   210
      
  AIC:               
 
  -1645.


  
Df Residuals:
          
   206
      
  BIC:               
 
  -1632.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
    0.0408
 
    4.568
 
    0.009
 
 0.993
 
   -8.966     9.048


  
total_cpg_no_filter
 
-5.724e-09
 
 4.51e-09
 
   -1.268
 
 0.206
 
-1.46e-08  3.17e-09


  
bsRate_mean
         
   -0.0226
 
    4.658
 
   -0.005
 
 0.996
 
   -9.206     9.161


  
avgReadCpG_mean
     
   -0.0011
 
    0.012
 
   -0.094
 
 0.925
 
   -0.024     0.022




  
Omnibus:
       
14.287
 
  Durbin-Watson:     
 
   1.921


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


  
Skew:
          
 0.618
 
  Prob(JB):          
 
0.000566


  
Kurtosis:
      
 2.575
 
  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 0x107628e80>






    

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






    

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






    

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_total, tritopool.PDR_total)), 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_total")
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_total

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.260


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
0.0225
 


  
Time:
                 
00:35:52
     
  Log-Likelihood:    
 
  801.55


  
No. Observations:
      
   210
      
  AIC:               
 
  -1595.


  
Df Residuals:
          
   206
      
  BIC:               
 
  -1582.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
    1.2956
 
    5.970
 
    0.217
 
 0.828
 
  -10.475    13.066


  
total_cpg_no_filter
 
 1.307e-08
 
 4.75e-09
 
    2.749
 
 0.007
 
  3.7e-09  2.24e-08


  
bsRate_mean
         
   -1.3146
 
    6.192
 
   -0.212
 
 0.832
 
  -13.522    10.892


  
avgReadCpG_mean
     
   -0.0020
 
    0.015
 
   -0.135
 
 0.893
 
   -0.032     0.028




  
Omnibus:
       
12.567
 
  Durbin-Watson:     
 
   2.501


  
Prob(Omnibus):
 
 0.002
 
  Jarque-Bera (JB):  
 
  10.852


  
Skew:
          
 0.478
 
  Prob(JB):          
 
 0.00440


  
Kurtosis:
      
 2.428
 
  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 0x10873f748>






    

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






    

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






    

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_total, cw154.PDR_total)), 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_total")
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_total

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.025


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.009


  
Method:
             
Least Squares
  
  F-statistic:       
 
   1.571


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
 0.198
 


  
Time:
                 
00:35:55
     
  Log-Likelihood:    
 
  719.16


  
No. Observations:
      
   190
      
  AIC:               
 
  -1430.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1417.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
    1.5147
 
    2.434
 
    0.622
 
 0.535
 
   -3.287     6.317


  
total_cpg_no_filter
 
-8.124e-09
 
 4.13e-09
 
   -1.969
 
 0.050
 
-1.63e-08  1.44e-11


  
bsRate_mean
         
   -1.5880
 
    2.548
 
   -0.623
 
 0.534
 
   -6.615     3.439


  
avgReadCpG_mean
     
    0.0032
 
    0.008
 
    0.411
 
 0.682
 
   -0.012     0.018




  
Omnibus:
       
14.219
 
  Durbin-Watson:     
 
   2.170


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


  
Skew:
          
 0.378
 
  Prob(JB):          
 
  0.0122


  
Kurtosis:
      
 2.264
 
  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 0x1087241d0>






    

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






    

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






    

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_total, cw154.PDR_total)), 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_total")
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_total

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.060


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.045


  
Method:
             
Least Squares
  
  F-statistic:       
 
   3.934


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
0.00942


  
Time:
                 
00:35:58
     
  Log-Likelihood:    
 
  702.58


  
No. Observations:
      
   190
      
  AIC:               
 
  -1397.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1384.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
   -8.6575
 
    2.946
 
   -2.939
 
 0.004
 
  -14.469    -2.846


  
total_cpg_no_filter
 
-5.206e-10
 
 3.64e-09
 
   -0.143
 
 0.886
 
 -7.7e-09  6.66e-09


  
bsRate_mean
         
    9.0310
 
    3.091
 
    2.922
 
 0.004
 
    2.933    15.129


  
avgReadCpG_mean
     
   -0.0016
 
    0.006
 
   -0.270
 
 0.788
 
   -0.014     0.010




  
Omnibus:
       
 9.844
 
  Durbin-Watson:     
 
   1.839


  
Prob(Omnibus):
 
 0.007
 
  Jarque-Bera (JB):  
 
   8.185


  
Skew:
          
 0.420
 
  Prob(JB):          
 
  0.0167


  
Kurtosis:
      
 2.427
 
  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 0x1096935c0>






    

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






    

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






    

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_total, cw154.PDR_total)), 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_total")
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_total

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.089


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.074


  
Method:
             
Least Squares
  
  F-statistic:       
 
   6.049


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
 
0.000596


  
Time:
                 
00:36:01
     
  Log-Likelihood:    
 
  651.12


  
No. Observations:
      
   190
      
  AIC:               
 
  -1294.


  
Df Residuals:
          
   186
      
  BIC:               
 
  -1281.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
   -3.8803
 
    4.876
 
   -0.796
 
 0.427
 
  -13.500     5.739


  
total_cpg_no_filter
 
 -2.13e-08
 
 6.06e-09
 
   -3.515
 
 0.001
 
-3.33e-08 -9.34e-09


  
bsRate_mean
         
    3.8863
 
    5.055
 
    0.769
 
 0.443
 
   -6.087    13.859


  
avgReadCpG_mean
     
    0.0224
 
    0.016
 
    1.362
 
 0.175
 
   -0.010     0.055




  
Omnibus:
       
10.262
 
  Durbin-Watson:     
 
   1.833


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


  
Skew:
          
 0.582
 
  Prob(JB):          
 
 0.00401


  
Kurtosis:
      
 2.805
 
  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 0x109e0cf98>






    

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






    

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






    

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_total, pcell.PDR_total)), 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_total")
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_total

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.030


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.008


  
Method:
             
Least Squares
  
  F-statistic:       
 
   1.357


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
 0.259
 


  
Time:
                 
00:36:04
     
  Log-Likelihood:    
 
  247.33


  
No. Observations:
      
   136
      
  AIC:               
 
  -486.7


  
Df Residuals:
          
   132
      
  BIC:               
 
  -475.0


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -30.6966
 
   57.991
 
   -0.529
 
 0.597
 
 -145.409    84.015


  
total_cpg_no_filter
 
-1.001e-07
 
 7.26e-08
 
   -1.378
 
 0.170
 
-2.44e-07  4.35e-08


  
bsRate_mean
         
   31.8723
 
   58.137
 
    0.548
 
 0.584
 
  -83.129   146.874


  
avgReadCpG_mean
     
   -0.1449
 
    0.075
 
   -1.923
 
 0.057
 
   -0.294     0.004




  
Omnibus:
       
 9.591
 
  Durbin-Watson:     
 
   1.827


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


  
Skew:
          
 0.671
 
  Prob(JB):          
 
 0.00561


  
Kurtosis:
      
 2.830
 
  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 0x10a4a8cc0>






    

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






    

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






    

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_total, pcell.PDR_total)), 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_total")
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_total

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.081


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.069


  
Method:
             
Least Squares
  
  F-statistic:       
 
   6.678


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
 
0.000244


  
Time:
                 
00:36:07
     
  Log-Likelihood:    
 
  520.12


  
No. Observations:
      
   231
      
  AIC:               
 
  -1032.


  
Df Residuals:
          
   227
      
  BIC:               
 
  -1018.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -67.7274
 
   42.852
 
   -1.580
 
 0.115
 
 -152.167    16.712


  
total_cpg_no_filter
 
 7.637e-08
 
 3.29e-08
 
    2.322
 
 0.021
 
 1.16e-08  1.41e-07


  
bsRate_mean
         
   67.2942
 
   42.926
 
    1.568
 
 0.118
 
  -17.290   151.879


  
avgReadCpG_mean
     
    0.0865
 
    0.020
 
    4.340
 
 0.000
 
    0.047     0.126




  
Omnibus:
       
39.953
 
  Durbin-Watson:     
 
   1.261


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


  
Skew:
          
 1.099
 
  Prob(JB):          
 
1.04e-12


  
Kurtosis:
      
 3.947
 
  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 0x10ad586a0>






    

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






    

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






    

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_total, pcell.PDR_total)), 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_total")
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_total

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.169


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.091


  
Method:
             
Least Squares
  
  F-statistic:       
 
   2.174


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
 0.110
 


  
Time:
                 
00:36:10
     
  Log-Likelihood:    
 
  77.381


  
No. Observations:
      
    36
      
  AIC:               
 
  -146.8


  
Df Residuals:
          
    32
      
  BIC:               
 
  -140.4


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
 -204.2335
 
  245.012
 
   -0.834
 
 0.411
 
 -703.306   294.839


  
total_cpg_no_filter
 
-1.981e-08
 
 3.49e-08
 
   -0.567
 
 0.575
 
 -9.1e-08  5.14e-08


  
bsRate_mean
         
  203.9084
 
  245.166
 
    0.832
 
 0.412
 
 -295.478   703.295


  
avgReadCpG_mean
     
    0.1252
 
    0.079
 
    1.594
 
 0.121
 
   -0.035     0.285




  
Omnibus:
       
 1.836
 
  Durbin-Watson:     
 
   1.647


  
Prob(Omnibus):
 
 0.399
 
  Jarque-Bera (JB):  
 
   1.718


  
Skew:
          
 0.484
 
  Prob(JB):          
 
   0.424


  
Kurtosis:
      
 2.541
 
  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 0x10b50c940>






    

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






    

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






    

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_total, pcell.PDR_total)), 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_total")
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_total

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.006


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
  -0.010


  
Method:
             
Least Squares
  
  F-statistic:       
 
  0.4045


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
  
 0.750
 


  
Time:
                 
00:36:14
     
  Log-Likelihood:    
 
  397.02


  
No. Observations:
      
   190
      
  AIC:               
 
  -786.0


  
Df Residuals:
          
   186
      
  BIC:               
 
  -773.1


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -13.3122
 
   46.675
 
   -0.285
 
 0.776
 
 -105.392    78.768


  
total_cpg_no_filter
 
 3.897e-08
 
 3.84e-08
 
    1.015
 
 0.311
 
-3.68e-08  1.15e-07


  
bsRate_mean
         
   13.2619
 
   46.664
 
    0.284
 
 0.777
 
  -78.796   105.320


  
avgReadCpG_mean
     
    0.0174
 
    0.050
 
    0.348
 
 0.728
 
   -0.081     0.116




  
Omnibus:
       
12.462
 
  Durbin-Watson:     
 
   1.572


  
Prob(Omnibus):
 
 0.002
 
  Jarque-Bera (JB):  
 
  13.739


  
Skew:
          
 0.650
 
  Prob(JB):          
 
 0.00104


  
Kurtosis:
      
 2.788
 
  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 0x10bc3ea90>






    

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






    

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






    

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_total, mcell.PDR_total)), 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_total")
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_total

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.249


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.234


  
Method:
             
Least Squares
  
  F-statistic:       
 
   16.45


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
 
2.77e-09


  
Time:
                 
00:36:17
     
  Log-Likelihood:    
 
  435.37


  
No. Observations:
      
   153
      
  AIC:               
 
  -862.7


  
Df Residuals:
          
   149
      
  BIC:               
 
  -850.6


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  233.6229
 
   38.347
 
    6.092
 
 0.000
 
  157.850   309.396


  
total_cpg_no_filter
 
-6.274e-08
 
 2.11e-08
 
   -2.969
 
 0.003
 
-1.04e-07  -2.1e-08


  
bsRate_mean
         
 -234.6229
 
   38.536
 
   -6.088
 
 0.000
 
 -310.771  -158.475


  
avgReadCpG_mean
     
    0.0566
 
    0.022
 
    2.548
 
 0.012
 
    0.013     0.101




  
Omnibus:
       
50.820
 
  Durbin-Watson:     
 
   1.593


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


  
Skew:
          
 1.586
 
  Prob(JB):          
 
2.51e-21


  
Kurtosis:
      
 5.197
 
  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 0x10c3fde10>






    

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






    

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






    

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_total, mcell.PDR_total)), 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_total")
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_total

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.107


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.089


  
Method:
             
Least Squares
  
  F-statistic:       
 
   5.958


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
 
0.000729


  
Time:
                 
00:36:21
     
  Log-Likelihood:    
 
  383.38


  
No. Observations:
      
   153
      
  AIC:               
 
  -758.8


  
Df Residuals:
          
   149
      
  BIC:               
 
  -746.6


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -16.9739
 
   60.736
 
   -0.279
 
 0.780
 
 -136.989   103.041


  
total_cpg_no_filter
 
 1.259e-07
 
 4.07e-08
 
    3.098
 
 0.002
 
 4.56e-08  2.06e-07


  
bsRate_mean
         
   17.0169
 
   60.905
 
    0.279
 
 0.780
 
 -103.333   137.366


  
avgReadCpG_mean
     
   -0.0038
 
    0.041
 
   -0.093
 
 0.926
 
   -0.085     0.078




  
Omnibus:
       
71.167
 
  Durbin-Watson:     
 
   0.230


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


  
Skew:
          
 2.063
 
  Prob(JB):          
 
6.18e-39


  
Kurtosis:
      
 6.253
 
  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 0x10cc2e0f0>






    

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






    

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






    

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_total, mcell.PDR_total)), 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_total")
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_total

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.192


  
Model:
                   
OLS
       
  Adj. R-squared:    
 
   0.173


  
Method:
             
Least Squares
  
  F-statistic:       
 
   10.43


  
Date:
             
Wed, 10 Aug 2016
 
  Prob (F-statistic):
 
3.30e-06


  
Time:
                 
00:36:24
     
  Log-Likelihood:    
 
  586.97


  
No. Observations:
      
   136
      
  AIC:               
 
  -1166.


  
Df Residuals:
          
   132
      
  BIC:               
 
  -1154.


  
Df Model:
              
     3
      
                     
     
 
   


  
Covariance Type:
      
nonrobust
    
                     
     
 
   




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


  
const
               
  -39.1096
 
    8.655
 
   -4.519
 
 0.000
 
  -56.230   -21.989


  
total_cpg_no_filter
 
-4.442e-10
 
 7.49e-09
 
   -0.059
 
 0.953
 
-1.53e-08  1.44e-08


  
bsRate_mean
         
   39.1438
 
    8.669
 
    4.516
 
 0.000
 
   21.997    56.291


  
avgReadCpG_mean
     
    0.0105
 
    0.006
 
    1.870
 
 0.064
 
   -0.001     0.022




  
Omnibus:
       
 3.164
 
  Durbin-Watson:     
 
   1.788


  
Prob(Omnibus):
 
 0.206
 
  Jarque-Bera (JB):  
 
   3.203


  
Skew:
          
 0.356
 
  Prob(JB):          
 
   0.202


  
Kurtosis:
      
 2.759
 
  Cond. No.          
 
1.28e+10