<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
Comparison of Hi-C experiments
Between conditions
HindIII
NcoI
Between replicates with different restriction enzymes
T0
T60
Merge Hi-C experiments
Normalizing merged data
Comparison of Hi-C experiments
Some settings
Step1:  Write a little function to load HiCData obeject
Step2:  Load data
ICE normalized
Step3:  DRY normalized
Step4:  Plot correlations
ICE
Step5:  DRY
Step6:  Repeat at 1 Mb resolution
Step7:  ICE
Step8:  DRY
Step9:  Compare eigenvectors
ICE
Step10:  DRY
Step11:  Merge Hi-C experiments
Once agreed that experiments are similar, they can be merged.
Here is a simple way to merge valid pairs. Arguably we may want to merge unfiltered data but the difference would be minimal specially with non-replicates.
Step12:  Normalizing merged data
Step13:  All in one loop to
<ASSISTANT_TASK:>
Python Code:
from pytadbit.mapping.analyze import eig_correlate_matrices, correlate_matrices
from pytadbit import load_hic_data_from_reads
from cPickle import load
from matplotlib import pyplot as plt
reso = 1000000
base_path = 'results/{0}/03_filtering/valid_reads12_{0}.tsv'
bias_ice_path = 'results/{1}/04_normalizing/biases_{0}_{1}.pick'
bias_dry_path = 'results/{1}/04_normalizing/biases_dryhic_{0}_{1}.tsv'
bads_path = 'results/{1}/04_normalizing/bad_columns_{0}_{1}.pick'
Explanation: Table of Contents
Comparison of Hi-C experiments
Between conditions
HindIII
NcoI
Between replicates with different restriction enzymes
T0
T60
Merge Hi-C experiments
Normalizing merged data
Comparison of Hi-C experiments
Some settings
End of explanation
def my_load_hic_data(renz, reso, which='ice'):
    hic_data = load_hic_data_from_reads(base_path.format(renz), resolution=reso)
    if which=='ice':
        hic_data.bias = load(open(bias_ice_path.format(reso, renz)))
    else:
        hic_data.bias = dict([(int(l.split()[0]), float(l.split()[1])) 
                              for l in open(bias_dry_path.format(reso, renz))])
    hic_data.bads = load(open(bads_path.format(reso, renz)))
    return hic_data
renz1 = 'HindIII'
renz2 = 'MboI'
Explanation: Write a little function to load HiCData obeject
End of explanation
hic_data1_ice = my_load_hic_data(renz1, reso, which='ice')
hic_data2_ice = my_load_hic_data(renz2, reso, which='ice')
Explanation: Load data
ICE normalized
End of explanation
hic_data1_dry = my_load_hic_data(renz1, reso, which='dry')
hic_data2_dry = my_load_hic_data(renz2, reso, which='dry')
%matplotlib inline
Explanation: DRY normalized
End of explanation
## this part is to "tune" the plot ##
plt.figure(figsize=(9, 6))
axe = plt.subplot()
axe.grid()
axe.set_xticks(range(0, 55, 5))
axe.set_xticklabels(['%d Mb' % int(i * 0.2) if i else '' for i in range(0, 55, 5)], rotation=-45)
#####################################
_ = correlate_matrices(hic_data1_ice, hic_data2_ice, max_dist=50, show=False, axe=axe, normalized=True)
Explanation: Plot correlations
ICE
End of explanation
## this part is to "tune" the plot ##
plt.figure(figsize=(9, 6))
axe = plt.subplot()
axe.grid()
axe.set_xticks(range(0, 55, 5))
axe.set_xticklabels(['%d Mb' % int(i * 0.2) if i else '' for i in range(0, 55, 5)], rotation=-45)
#####################################
_ = correlate_matrices(hic_data1_dry, hic_data2_dry, max_dist=50, show=False, axe=axe, normalized=True)
Explanation: DRY
End of explanation
reso = 1000000
hic_data1_ice = my_load_hic_data(renz1, reso, which='ice')
hic_data2_ice = my_load_hic_data(renz2, reso, which='ice')
hic_data1_dry = my_load_hic_data(renz1, reso, which='dry')
hic_data2_dry = my_load_hic_data(renz2, reso, which='dry')
Explanation: Repeat at 1 Mb resolution
End of explanation
## this part is to "tune" the plot ##
plt.figure(figsize=(9, 6))
axe = plt.subplot()
axe.grid()
axe.set_xticks(range(0, 55, 5))
axe.set_xticklabels(['%d Mb' % int(i * 0.2) if i else '' for i in range(0, 55, 5)], rotation=-45)
#####################################
_ = correlate_matrices(hic_data1_ice, hic_data2_ice, max_dist=50, show=False, axe=axe, normalized=True)
Explanation: ICE
End of explanation
## this part is to "tune" the plot ##
plt.figure(figsize=(9, 6))
axe = plt.subplot()
axe.grid()
axe.set_xticks(range(0, 55, 5))
axe.set_xticklabels(['%d Mb' % int(i * 0.2) if i else '' for i in range(0, 55, 5)], rotation=-45)
#####################################
_ = correlate_matrices(hic_data1_dry, hic_data2_dry, max_dist=50, show=False, axe=axe, normalized=True)
Explanation: DRY
End of explanation
corrs = eig_correlate_matrices(hic_data1_ice, hic_data2_ice, show=True, aspect='auto', normalized=True)
for cor in corrs:
    print ' '.join(['%5.3f' % (c) for c in cor]) + '\n'
Explanation: Compare eigenvectors
ICE
End of explanation
corrs = eig_correlate_matrices(hic_data1_dry, hic_data2_dry, show=True, aspect='auto', normalized=True)
for cor in corrs:
    print ' '.join(['%5.3f' % (c) for c in cor]) + '\n'
Explanation: DRY
End of explanation
from pytadbit.mapping import merge_2d_beds
! mkdir -p results/both
! mkdir -p results/both/
! mkdir -p results/both03_filtering/
! mkdir -p results/both/03_filtering/
renz1 = 'HindIII'
renz2 = 'MboI'
hic_data1 = 'results/{0}/03_filtering/valid_reads12_{0}.tsv'.format(renz1)
hic_data2 = 'results/{0}/03_filtering/valid_reads12_{0}.tsv'.format(renz2)
hic_data  = 'results/both/03_filtering/valid_reads12.tsv'
merge_2d_beds(hic_data1, hic_data2, hic_data)
Explanation: Merge Hi-C experiments
Once agreed that experiments are similar, they can be merged.
Here is a simple way to merge valid pairs. Arguably we may want to merge unfiltered data but the difference would be minimal specially with non-replicates.
End of explanation
from pytadbit.mapping.analyze import hic_map
from cPickle import dump
! mkdir -p results/both/04_normalizing
Explanation: Normalizing merged data
End of explanation
for reso in [1000000, 300000, 100000]:
    print '   *', reso
    # load hic_data
    hic_data = load_hic_data_from_reads(
        'results/both/03_filtering/valid_reads12.tsv', reso)
    # filter columns
    hic_data.filter_columns(draw_hist=False, min_count=10, by_mean=True)
    # normalize
    hic_data.normalize_hic(iterations=0)
    # save biases to reconstruct normalization
    out = open('results/both/04_normalizing/biases_{0}.pick'.format(reso), 'w')
    dump(hic_data.bias, out)
    out.close()
    # save filtered out columns
    out = open('results/both/04_normalizing/bad_columns_{0}.pick'.format(reso), 'w')
    dump(hic_data.bads, out)
    out.close()
    # save data as raw matrix per chromsome
    hic_map(hic_data, by_chrom='intra', normalized=False,
            savedata='results/both/04_normalizing/{0}_raw'.format(reso))
    # save data as normalized matrix per chromosome
    hic_map(hic_data, by_chrom='intra', normalized=True,
            savedata='results/both/04_normalizing/{0}_norm'.format(reso))
    # if the resolution is low save the full genomic matrix
    if reso > 500000:
        hic_map(hic_data, by_chrom=False, normalized=False, 
                savefig ='results/both/04_normalizing/{0}_raw.png'.format(reso),
                savedata='results/both/04_normalizing/{0}_raw.mat'.format(reso))
        hic_map(hic_data, by_chrom=False, normalized=True,
                savefig ='results/both/04_normalizing/{0}_norm.png'.format(reso) ,
                savedata='results/both/04_normalizing/{0}_norm.mat'.format(reso))
Explanation: All in one loop to:
 - filter
 - normalize
 - generate intra-chromosome and genomic matrices
all at diferent resolutions, and for both time points
End of explanation
<END_TASK>