<jupyter_start><jupyter_text>Get SMC positions data<jupyter_code>base_folder_string = 'output_0310_tiling600_2populationNonBypassingLEFs_GapMethod_SparseDSB_FullParamter_DenseDSBs_run'
folders = [foldername for foldername in os.listdir('Data/') if base_folder_string in foldername]

# Parameter combination giving more than 95% synapsis efficiency
par_combo = ['sep125_proc250_superportion0.2_superproc5000_ctcf16_dsb4_superloading1000_bs0.5']

time_point_list = [0]
box_size = 2500 # number of monomers on each side of the DSB that we will make the contact map for
short_list_beforeDSB = []
for m in range(len(folders)):
    for n in range(len(par_combo)):
        data_dir = 'Data/'+folders[m]
        par = par_combo[n]
        smc_file = data_dir + '/LEFcoordinates_' + par + '.npy'
        SMCs = np.load(smc_file)

        dsb_file = data_dir +'/DSB_boundary_coordinates.npy'
        with open(dsb_file, 'rb') as infile:
            DSBs = np.load(infile)
            BEs = np.load(infile)

        smcs_dict = {}
        for t in time_point_list:
            smcs_dict[t] = [(int(l),int(r)) for (l,r) in zip(SMCs[2*t,:],SMCs[2*t+1,:])]

        for d in DSBs[1::2]: # d is the right DSB end
            
            shift_factor = d-box_size

            for t in time_point_list: #smcs_dict.keys():
                smcs = [(x[0]-shift_factor,x[1]-shift_factor) for x in smcs_dict[t] if (np.abs(x[0]-d)<=box_size)  and (np.abs(x[1]-d)<box_size)]
                #short_list.append(smcs)
                has_constraint = False
                for x in smcs:
                    if x[0]<box_size and x[1]>box_size:
                        has_constraint = True
                        break
                if has_constraint==True:
                    short_list_beforeDSB.append(smcs) 
base = 50
nreps = 2
for si, smcs in enumerate(short_list_beforeDSB):
    dsb = [box_size-1,box_size]
    L_tot = 2*box_size + base
    
    for rep in range(nreps):

        m_samples1000_beforeDSB, c_samples1000_beforeDSB = contact_map_generator_from_SMC_list_and_DSBs(smcs,dsb,L_tot,base=base,nsamples=1000)

        if rep == 0 and si == 0:
            M_samples1000_beforeDSB = m_samples1000_beforeDSB
            C_samples1000_beforeDSB = c_samples1000_beforeDSB
        else:
            M_samples1000_beforeDSB += m_samples1000_beforeDSB
            C_samples1000_beforeDSB += c_samples1000_beforeDSB
        
M_new_samples1000_beforeDSB = np.zeros((100,100))
C_new_samples1000_beforeDSB = np.zeros((100,100))
for x in range(100):
    for y in range(100):
        M_new_samples1000_beforeDSB[x//1,y//1] += M_samples1000_beforeDSB[x,y]
        C_new_samples1000_beforeDSB[x//1,y//1] += C_samples1000_beforeDSB[x,y]
base_folder_string = 'output_0310_tiling600_2populationNonBypassingLEFs_GapMethod_SparseDSB_FullParamter_DenseDSBs_run'
folders = [foldername for foldername in os.listdir('Data/') if base_folder_string in foldername]

# Parameter combination giving more than 95% synapsis efficiency
par_combo = ['sep125_proc250_superportion0.2_superproc5000_ctcf16_dsb4_superloading1000_bs0.5']

time_point_list = [6]
box_size = 2500 # number of monomers on each side of the DSB that we will make the contact map for
short_list = []
for m in range(len(folders)):
    for n in range(len(par_combo)):
        data_dir = 'Data/'+folders[m]
        par = par_combo[n]
        smc_file = data_dir + '/LEFcoordinates_' + par + '.npy'
        SMCs = np.load(smc_file)

        dsb_file = data_dir +'/DSB_boundary_coordinates.npy'
        with open(dsb_file, 'rb') as infile:
            DSBs = np.load(infile)
            BEs = np.load(infile)

        smcs_dict = {}
        for t in time_point_list:
            smcs_dict[t] = [(int(l),int(r)) for (l,r) in zip(SMCs[2*t,:],SMCs[2*t+1,:])]

        for d in DSBs[1::2]: # d is the right DSB end

            shift_factor = d-box_size

            for t in time_point_list: #smcs_dict.keys():
                smcs = [(x[0]-shift_factor,x[1]-shift_factor) for x in smcs_dict[t] if (np.abs(x[0]-d)<=box_size)  and (np.abs(x[1]-d)<box_size)]

                has_constraint = False
                for x in smcs:
                    if x[0]<box_size and x[1]>box_size:
                        has_constraint = True
                        break
                if has_constraint==True:
                    short_list.append(smcs) 
# number of DSBs in the contact map
len(short_list)
base = 50
nreps = 2
for si, smcs in enumerate(short_list):
    dsb = [box_size-1,box_size] 
    L_tot = 2*box_size + base
    
    for rep in range(nreps):
        m_samples1000, c_samples1000 = contact_map_generator_from_SMC_list_and_DSBs(smcs,dsb,L_tot,base=base,nsamples=1000)
        if rep == 0 and si == 0:
            M_samples1000_postDSB = m_samples1000
            C_samples1000_postDSB = c_samples1000
        else:
            M_samples1000_postDSB += m_samples1000
            C_samples1000_postDSB += c_samples1000
        
M_new_samples1000 = np.zeros((100,100))
C_new_samples1000 = np.zeros((100,100))
for x in range(100):
    for y in range(100):
        M_new_samples1000[x//1,y//1] += M_samples1000_postDSB[x,y]
        C_new_samples1000[x//1,y//1] += C_samples1000_postDSB[x,y]
raw=(M_new_samples1000/C_new_samples1000)/(M_new_samples1000_beforeDSB/C_new_samples1000_beforeDSB)
corrected=iterative_correction_symmetric(raw)
plt.imshow(np.log2(corrected[0]),cmap='cividis',vmax=0.4,vmin=-0.4)
# plt.hlines(50,0,100,linestyles=':')
# plt.vlines(50,0,100,linestyles=':')
plt.xlim([0,100])
plt.ylim([0,100])

plt.colorbar() 
plt.title("With DSB")
plt.savefig('Figures/'+'HiC_120min_2rep_1000samples_DividedByBeforeBSB_allTADs_vmax0.4_vmin-0.4_3rdParameterCombo_IterativeCorrected.pdf',format='pdf',bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10961272475795297, 'fraction_numerical': 0.059301521438450897, 'mean_word_length': 3.065354884047786, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '218055', 'license': 'MIT', 'n_tokens_mistral': 2412, 'n_tokens_neox': 2114, 'n_words': 385, 'path': 'GeneratingContactMaps.ipynb', 'repo_name': 'ahansenlab/DNA_break_synapsis_models', 'size': '5784'}