demo.py

#!/usr/bin/env python
# coding: utf-8

# # imports

# In[24]:


import pandas as pd
import matplotlib.pyplot as plt 
import seaborn as sns 
import os
import numpy as np
from sklearn.preprocessing import normalize
import scipy.cluster.hierarchy as shc
import numpy.ma as ma
np.seterr(divide='ignore', invalid='ignore')


# In[25]:


# import R to be used in python 
# sorry R people that may be reading this code 
import rpy2
get_ipython().run_line_magic('load_ext', 'rpy2.ipython')


# In[26]:


# set white seaborn plotting style 
sns.set(context='talk', style='white', rc={'figure.facecolor':'white'}, font_scale=1)
sns.set_style('ticks')


# # read in data

# In[27]:


# tissue proportion estimates generated by CelFiE
tissue_proportions = pd.read_csv("cfDNA/ONT_playing/clean_dorado/bed/celfie_output/1_tissue_proportions.txt", delimiter="\t")


# In[28]:


# rename column for nice plotting 
tissue_proportions.rename(columns={"Unnamed: 0": "samples"}, inplace=True)


# In[29]:


# change the orientation of the data so that it is easy to plot with seaborn 
tissue_proportions = tissue_proportions.melt("samples",  var_name="tissue", value_name="estimate")




# In[31]:


plt.hist(tissue_proportions, )


# In[32]:


# mean values of placenta 
placenta_estimates = tissue_proportions[tissue_proportions["tissue"] == "placenta"]
placenta_estimates.groupby("status")["estimate"].mean()


# # calculate p-value in R 

# In[33]:


# save it to be read by R; can also just do this stand alone in R 
tissue_proportions.to_csv("celfie_demo/sample_output/tissue_proportions.csv", index=False)


# In[34]:


get_ipython().run_cell_magic('R', '', '# R code to calculate the p-value between the estimates of preg/not preg women\n\ndf = read.csv("celfie_demo/sample_output/tissue_proportions.csv")  # re-read in the dataframe for this R code \ndf$status <- as.character(df$status)\ndf$status[df$status == "not pregnant"] <- 1\ndf$status[df$status == "pregnant"] <- 0\ndf$status <- as.numeric(df$status)\n\nplacenta = subset(df, df$tissue=="placenta")\n\nsummary(glm(status ~ estimate, data = placenta))\n')


# # methylation proportions 

# ## reference methylation 

# In[35]:


reference = pd.read_csv("celfie_demo/sample_data.txt", delimiter="\t").iloc[:, 33:]  # read in reference only 


# In[36]:


reference_meth = reference.iloc[:, 3::2].values  # methylation values
reference_depth = reference.iloc[:, 4::2].values  # depth values


# In[37]:


reference_prop = reference_meth/reference_depth  # convert to proportion 


# ## estimated methylation by CelFiE

# In[38]:


estimated = pd.read_csv("celfie_demo/sample_output/1_methylation_proportions.txt", delimiter="\t")


# In[39]:


estimated_prop = estimated.iloc[:, 1:].values 


# In[40]:


# calculate the correlation between the reference data and that estimated by CelFiE

correlation = []

for i in range(reference_prop.shape[1]): 
        corr = ma.corrcoef(ma.masked_invalid(reference_prop[:, i]), 
                          ma.masked_invalid(estimated_prop[:, i]))[0, 1]  # mask Nans
        
        correlation.append(corr)
        
    


# # plot correlation 

# In[41]:


sns.barplot(x=list(estimated)[1:-1], y=correlation, color="#61c2d3")
plt.xticks(rotation=90)
plt.xlabel("tissue")
plt.ylabel("correlation")
plt.show()


# # hierarchical clustering 

# In[42]:


# append the unknown to the reference data so that they can be clustered together
reference_unknown = np.append(reference_prop, estimated_prop[:, -1].reshape(-1, 1), axis=1)


# In[43]:


# drop any missing values 
reference_unknown = pd.DataFrame(reference_unknown)
reference_unknown = reference_unknown.dropna()


# In[44]:


reference_unknown_scaled = normalize(reference_unknown)  # normalize data 
reference_unknown_scaled = pd.DataFrame(reference_unknown_scaled, columns=list(estimated)[1:]) # make dataframe again 


# In[45]:


sns.set_context("poster")
plt.figure(figsize=(10, 7))  
plt.title("Dendrograms")  
shc.set_link_color_palette(['black', 'black'])
dend = shc.dendrogram(shc.linkage(reference_unknown_scaled.T, method='ward'), orientation='left', 
                      labels=reference_unknown_scaled.columns, leaf_font_size = 20)
plt.title("ALS unknown")
plt.show()


# In[ ]: