<jupyter_start><jupyter_text>__Author__: <NAME>__Email__: <EMAIL>__Date__:3/4/2020__Platform__: Python 2.7<jupyter_code>import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['font.size']=15
matplotlib.rcParams['font.family']='Arial'
import matplotlib.pylab as plt
import numpy as np
import os

import cPickle as pickle
import matplotlib.pylab as plt
import numpy as np
from scipy.spatial.distance import pdist,cdist,squareform
from tqdm import tqdm_notebook as tqdm<jupyter_output><empty_output><jupyter_text>Please specifiy the data/save folders<jupyter_code>data_folder = r'\data' #This is the folder containing the .tsv data files
save_data = r'\data'<jupyter_output><empty_output><jupyter_text>1. Load the positions of the ~1000 chromatin loci imaged across the ~5400 cells<jupyter_code>folder = data_folder

experiment = []
fid = open(folder+os.sep+r'genomic-scale.tsv','r')
lines = np.array([ln[:-1].split('\t')for ln in fid if len(ln)>0])
head = list(lines[0])
experiment = np.concatenate([experiment,lines[1::2082,head.index('experiment number')].astype(int)])
zxy = np.array(lines[1:,:3][:],dtype=np.float)
dLAM = np.array(lines[1:,-1].astype(float))

fid = open(folder+os.sep+r'genomic-scale-with transcription and nuclear bodies.tsv','r')
lines = np.array([ln[:-1].split('\t')for ln in fid if len(ln)>0])
head = list(lines[0])
experiment = np.concatenate([experiment,lines[1::2082,head.index('experiment number')].astype(int)])
dLAM = np.concatenate([dLAM,np.array(lines[1:,-3].astype(float))])
zxy = np.concatenate([zxy,np.array(lines[1:,:3][:],dtype=np.float)])
zxy = zxy.reshape([-1,2082,3])/1000 #transform to um
dLAM = dLAM.reshape([-1,2082])/1000<jupyter_output><empty_output><jupyter_text>2. Calculate median distance matrices and proximity matrices for cis- and trans- interactions<jupyter_code>lens = [76, 80, 66, 63, 60, 55, 53, 48, 40, 43, 44, 44, 33, 30, 31, 30, 33, 33, 33, 33, 31, 31, 51]
edges = [0]+list(np.cumsum(lens))
ijs = []
for i in range(len(lens)):
    for j in range(len(lens)):
        ijs.append((i,j))
im_med = np.zeros([edges[-1],edges[-1]])
cut_offs = [0.25,0.5,0.75,1]
#cut_offs = np.arange(0.1,1.05,0.05)
im_fr = np.zeros([edges[-1],edges[-1],len(cut_offs)])
im_med_trans = []
im_med_cis = []
im_fr_trans = [[] for _ in cut_offs]
im_fr_cis = [[] for _ in cut_offs]

#Decide whether to add simulated noise
#rnd_noise = np.random.normal(scale=60/1.6,size=zxy.shape)/1000
zxy_ = zxy#+rnd_noise
for i,j in tqdm(ijs):
    arr = []
    for st1 in [0,edges[-1]]:
        for st2 in [0,edges[-1]]:
            zxy1 = zxy_[:,st1+edges[i]:st1+edges[i+1]]
            zxy2 = zxy_[:,st2+edges[j]:st2+edges[j+1]]
            arr =arr+[cdist(zxy1[k],zxy2[k]) for k in range(len(zxy1))]
    arr = np.array(arr)
    im_med[edges[i]:edges[i+1],edges[j]:edges[j+1]]=np.nanmedian(arr,axis=0)
    if i==j:
        im_med_cis.append(np.nanmedian(arr[::2],axis=0))
        im_med_trans.append(np.nanmedian(arr[1::2],axis=0))
    for ic,cutoff in enumerate(cut_offs):
        im_fr[edges[i]:edges[i+1],edges[j]:edges[j+1],ic] = 1.*np.sum(arr<cutoff,0)/np.sum(arr>-1,0)
        if i==j:
            im_fr_trans[ic].append(1.*np.sum(arr[1::2]<cutoff,0)/np.sum(arr[1::2]>-1,0))
            im_fr_cis[ic].append(1.*np.sum(arr[::2]<cutoff,0)/np.sum(arr[::2]>-1,0))
pickle.dump([im_med,im_fr,im_med_trans,im_med_cis,im_fr_trans,im_fr_cis,len(zxy)],
        open(save_data+r'/mat_contact_IMR90_untreated.pkl','wb'))<jupyter_output><empty_output><jupyter_text>Load saved data<jupyter_code>im_med,im_fr,im_med_trans,im_med_cis,im_fr_trans,im_fr_cis,lzxy = pickle.load(
            open(save_data+r'/mat_contact_IMR90_untreated.pkl','rb'))
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['font.size']=15
matplotlib.rcParams['font.family']='Arial'
import matplotlib.pylab as plt
import numpy as np
import os<jupyter_output><empty_output><jupyter_text>Display matrices<jupyter_code>#number of regions per chormosome
lens = [76, 80, 66, 63, 60, 55, 53, 48, 40, 43, 44, 44, 33, 30, 31, 30, 33, 33, 33, 33, 31, 31, 51]
edges = np.cumsum([0]+lens)
colors_all = np.array([1.*il/len(lens)for il,ln in enumerate(lens) for _ in range(ln)])

from mpl_toolkits.axes_grid1 import make_axes_locatable
fig = plt.figure(figsize=(10,15))
plt.title('Median spatial distance\n('+str(lzxy)+' cells)')
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", size="5%", pad=0.)
cax2 = divider.append_axes("left", size="5%", pad=0.)
ax.imshow(im_med,cmap='coolwarm_r',vmin=3,vmax=10)
deltah = 50
cax.imshow(np.array([colors_all]*deltah),cmap='rainbow_r',interpolation='nearest')
cax2.imshow(np.array([colors_all]*deltah).T,cmap='rainbow_r',interpolation='nearest')
xticks = np.where(np.diff(colors_all)!=0)[0]
for tk in xticks: cax.plot([tk,tk],[0.5,deltah-0.5],'k')
for tk in xticks: cax2.plot([0.5,deltah-0.5],[tk,tk],'k')
centers = (np.cumsum([0]+lens)[:-1]+np.cumsum(lens))/2.
nms = map(str,np.arange(len(centers))+1)
nms[-1]='X'
cax.set_xticks(centers)
cax.set_xticklabels(nms)
cax.set_xlabel('Chromosomes')
cax2.set_yticks(centers)
cax2.set_yticklabels(nms)
cax2.set_ylabel('Chromosomes')
ax.set_xticks(xticks)
ax.set_xticklabels([''])
ax.set_yticks([])
cax.set_yticks([])
cax2.set_xticks([])

#fig.savefig(r'C:\Users\Bogdan\Dropbox\ScienceWGCT_Figures\baseimages\Fig1_med_map_all_nocbar.pdf')
#number of regions per chormosome
lens = [76, 80, 66, 63, 60, 55, 53, 48, 40, 43, 44, 44, 33, 30, 31, 30, 33, 33, 33, 33, 31, 31, 51]
edges = np.cumsum([0]+lens)
colors_all = np.array([1.*il/len(lens)for il,ln in enumerate(lens) for _ in range(ln)])

from mpl_toolkits.axes_grid1 import make_axes_locatable
fig = plt.figure(figsize=(10,15))
plt.title('Proximity frequency matrix\n('+str(lzxy)+' cells)')
ax = fig.add_subplot(111)
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", size="5%", pad=0.)
cax2 = divider.append_axes("left", size="5%", pad=0.)
ax.imshow(np.log(im_fr[:,:,1]),cmap='seismic',vmin=-8,vmax=-2.5)
deltah = 50
cax.imshow(np.array([colors_all]*deltah),cmap='rainbow_r',interpolation='nearest')
cax2.imshow(np.array([colors_all]*deltah).T,cmap='rainbow_r',interpolation='nearest')
xticks = np.where(np.diff(colors_all)!=0)[0]
for tk in xticks: cax.plot([tk,tk],[0.5,deltah-0.5],'k')
for tk in xticks: cax2.plot([0.5,deltah-0.5],[tk,tk],'k')
centers = (np.cumsum([0]+lens)[:-1]+np.cumsum(lens))/2.
nms = map(str,np.arange(len(centers))+1)
nms[-1]='X'
cax.set_xticks(centers)
cax.set_xticklabels(nms)
cax.set_xlabel('Chromosomes')
cax2.set_yticks(centers)
cax2.set_yticklabels(nms)
cax2.set_ylabel('Chromosomes')
ax.set_xticks(xticks)
ax.set_xticklabels([''])
ax.set_yticks([])
cax.set_yticks([])
cax2.set_xticks([])<jupyter_output><empty_output><jupyter_text>3. Compare population average maps deverived from Hi-C to those derived from imaging at the genomic-scale  Import the corresponding Hi-C data to the 1041 regionsNote: The Hi-C, from the combined IMR90 Rao et al 2019 paper, reports the number of raw contacts dectected between regions 500kb in size and centered at the coordinates of each imaged region.<jupyter_code>#This data has been downloaded from Rao et al. Cell 2014
hic_fl = data_folder+os.sep+r'\Hi-C matrices\Hi-C_contacts_genome-scale.tsv'
hic = np.array([ln[:-1].split('\t') for ln in open(hic_fl,'r')])
hic = hic[1:,1:].astype(float)<jupyter_output><empty_output><jupyter_text>Calculate cis cross-correlation with HiC<jupyter_code>#number of regions per chormosome
lens = [76, 80, 66, 63, 60, 55, 53, 48, 40, 43, 44, 44, 33, 30, 31, 30, 33, 33, 33, 33, 31, 31, 51]
edges = np.cumsum([0]+lens)
im_med,im_fr,im_med_trans,im_med_cis,im_fr_trans,im_fr_cis,nlen=  pickle.load(
    open(save_data+r'/mat_contact_IMR90_untreated.pkl','rb'))
xhic = []
ymed = []
yfr = []
for ichr in range(len(edges)-1):
    xhic.extend(hic[edges[ichr]:edges[ichr+1],edges[ichr]:edges[ichr+1]].ravel())
    ymed.extend(im_med_cis[ichr].ravel())
    im_ = im_fr_cis[1][ichr].copy()
    im_[np.arange(len(im_)),np.arange(len(im_))]=np.nan
    yfr.extend(im_.ravel())
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['font.size']=50
matplotlib.rcParams['font.family']='Arial'
import matplotlib.pylab as plt
import numpy as np
import os
def nan_corr_coef(x_,y_):
    x=np.ravel(x_)
    y=np.ravel(y_)
    keep=(np.isinf(x)==False)&(np.isinf(y)==False)&(np.isnan(x)==False)&(np.isnan(y)==False)
    x=x[keep]
    y=y[keep]
    return np.corrcoef([x,y])[0,1]
#fig,ax = plt.subplots(figsize=(24,24))
fig,ax = plt.subplots(figsize=(24,24))
rho = nan_corr_coef(xhic,yfr)
plt.loglog(xhic,yfr,'o',color='k',markersize=25,markeredgecolor='k',label='p='+str(np.round(rho,2)),alpha=0.005);
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
ax.set_ylabel('Imaging proximity frequency')
ax.set_xlabel('Hi-C number of contacts')
#plt.axis('equal')
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18790976775197243, 'fraction_numerical': 0.05167240804533837, 'mean_word_length': 6.41962077493817, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '209089', 'license': 'MIT', 'n_tokens_mistral': 4216, 'n_tokens_neox': 3840, 'n_words': 618, 'path': 'combinatorial_tracing/PostAnalysis/Part1__MedianDistanceAndProximityMaps.ipynb', 'repo_name': 'ZhuangLab/Chromatin_Analysis_2020_cell', 'size': '9018'}