<jupyter_start><jupyter_text>0. required packages for h5py<jupyter_code>import h5py
from ImageAnalysis3.classes import _allowed_kwds
import ast<jupyter_output><empty_output><jupyter_text>1. Create field-of-view class<jupyter_code>reload(ia)
reload(classes)
reload(classes.batch_functions)
reload(classes.field_of_view)
reload(io_tools.load)

reload(visual_tools)
reload(ia.correction_tools)
reload(ia.correction_tools.alignment)
reload(ia.spot_tools.matching)
reload(ia.segmentation_tools.chromosome)
reload(ia.spot_tools.fitting)

fov_param = {'data_folder':r'\\10.245.74.158\Chromatin_NAS_6\20200707-IMR90_SI16-5kb',
             'save_folder':r'W:\Pu_Temp\20200707_IMR90_5kb_SI13',
             'experiment_type': 'DNA',
             'num_threads': 12,
             'correction_folder':r'\\10.245.74.158\Chromatin_NAS_0\Corrections\20200803-Corrections_3color',
             'shared_parameters':{
                     'single_im_size':[30,2048,2048],
                     'corr_channels':['750','647','561'],
                     'num_empty_frames': 0, 
                     'corr_hot_pixel':True,
                     'corr_Z_shift':False,
                     'min_num_seeds':200,
                     'max_num_seeds': 400,
                     'spot_seeding_th':200,
                     'normalize_intensity_local':True,
                     'normalize_intensity_background':False,
                 },             
             }

fov = classes.field_of_view.Field_of_View(fov_param, _fov_id=3,
                                          _color_info_kwargs={
                                              '_color_filename':'Color_Usage',
                                          }, 
                                          _prioritize_saved_attrs=False,
                                          )
#fov._load_correction_profiles()<jupyter_output>Get Folder Names: (ia.get_img_info.get_folders)
- Number of folders: 37
- Number of field of views: 81
- Importing csv file: \\10.245.74.158\Chromatin_NAS_6\20200707-IMR90_SI16-5kb\Analysis\Color_Usage.csv
- header: ['Hyb', '750', '647', '561', '488', '405']
-- Hyb H0R0 exists in this data
-- DAPI exists in hyb: H0R0
- 32 folders are found according to color-usage annotation.
++ load bleed correction profile from original file.
-- loading bleedthrough correction profile from file:bleedthrough_correction_750_647_561_2048_2048.npy
++ load chromatic correction profile from original file.
-- loading chromatic correction profile from file:
	 750 chromatic_correction_750_647_30_2048_2048.npy
	 647 None
	 561 chromatic_correction_561_647_30_2048_2048.npy
++ load chromatic_constants correction profile from original file.
-- loading chromatic_constants correction profile from file:
	 750 chromatic_correction_750_647_30_2048_2048_const.pkl
	 647 None
	 561 chromatic_correction_561_647_30_2048_20[...]<jupyter_text>2. Process image into candidate spots<jupyter_code>reload(io_tools.load)
reload(spot_tools.fitting)
reload(correction_tools.chromatic)
reload(classes.batch_functions)

# process image into spots
id_list, spot_list = fov._process_image_to_spots('unique', 
                                                 #_sel_ids=np.arange(41,47),
                                                 _load_common_reference=True,
                                                 _load_with_multiple=False,
                                                 _save_images=True,
                                                 _warp_images=False, 
                                                 _overwrite_drift=False,
                                                 _overwrite_image=False,
                                                 _overwrite_spot=False,
                                                 _verbose=True)<jupyter_output>-- No folder selected, allow processing all 32 folders
-- checking unique, region:[1 2 3] in 0.011s.
-- checking unique, region:[4 5 6] in 0.010s.
-- checking unique, region:[7 8 9] in 0.008s.
-- checking unique, region:[10 11 12] in 0.008s.
-- checking unique, region:[13 14 15] in 0.009s.
-- checking unique, region:[16 17 18] in 0.009s.
-- checking unique, region:[19 20 21] in 0.009s.
-- checking unique, region:[22 23 24] in 0.009s.
-- checking unique, region:[25 26 27] in 0.009s.
-- checking unique, region:[28 29 30] in 0.009s.
-- checking unique, region:[31 32 33] in 0.009s.
-- checking unique, region:[34 35 36] in 0.008s.
-- checking unique, region:[37 38 39] in 0.008s.
-- checking unique, region:[40 41 42] in 0.008s.
-- checking unique, region:[43 44 45] in 0.009s.
-- checking unique, region:[46 47 48] in 0.008s.
-- checking unique, region:[49 50 51] in 0.009s.
-- checking unique, region:[52 53 54] in 0.009s.
-- checking unique, region:[55 56 57] in 0.009s.
-- checking unique, reg[...]<jupyter_text>3. Find chromosomes  3.1 load chromosome image<jupyter_code>chrom_im = fov._load_chromosome_image(_type='forward', _overwrite=False)
visual_tools.imshow_mark_3d_v2([chrom_im])<jupyter_output><empty_output><jupyter_text>3.2 find candidate chromosomes<jupyter_code>chrom_coords = fov._find_candidate_chromosomes_by_segmentation(_binary_per_th=99.5, 
                                                               _overwrite=True)<jupyter_output>-- adjust seed image with filter size=4
-- binarize image with threshold: 99.5%
-- erosion and dialation with size=1.
-- find close objects.
-- random walk segmentation, beta=10.<jupyter_text>3.3 select among candidate chromosomes<jupyter_code>fov._load_from_file('unique')

chrom_coords = fov._select_chromosome_by_candidate_spots(_good_chr_loss_th=0.3,
                                                         _cand_spot_intensity_th=1,
                                                         _save=True, 
                                                         _overwrite=True)<jupyter_output>+ loading unique from file: W:\Pu_Temp\20200707_IMR90_5kb_SI13\Conv_zscan_03.hdf5
- start select from 144 chromosomes with loss threshold=0.3
-- remove chr id 35, percentage of lost rounds:0.968.
-- remove chr id 44, percentage of lost rounds:0.946.
-- remove chr id 76, percentage of lost rounds:0.946.
-- remove chr id 138, percentage of lost rounds:0.946.
-- remove chr id 130, percentage of lost rounds:0.935.
-- remove chr id 138, percentage of lost rounds:0.935.
-- remove chr id 136, percentage of lost rounds:0.925.
-- remove chr id 107, percentage of lost rounds:0.914.
-- remove chr id 3, percentage of lost rounds:0.903.
-- remove chr id 18, percentage of lost rounds:0.903.
-- remove chr id 56, percentage of lost rounds:0.903.
-- remove chr id 94, percentage of lost rounds:0.903.
-- remove chr id 50, percentage of lost rounds:0.892.
-- remove chr id 125, percentage of lost rounds:0.892.
-- remove chr id 117, percentage of lost rounds:0.882.
-- remove chr id 128, percentage of lost r[...]<jupyter_text>visualize chromosomes selections<jupyter_code>%matplotlib notebook
%matplotlib notebook
## visualize
coord_dict = {'coords':[np.flipud(_coord) for _coord in fov.chrom_coords],
              'class_ids':list(np.zeros(len(fov.chrom_coords),dtype=np.int)),
             }

visual_tools.imshow_mark_3d_v2([fov.chrom_im], 
                               given_dic=coord_dict,
                               save_file=None,
                               )<jupyter_output><empty_output><jupyter_text>select spots based on chromosomes<jupyter_code>fov._load_from_file('unique')
intensity_th = 0.25
from ImageAnalysis3.spot_tools.picking import assign_spots_to_chromosomes

kept_spots_list = []
for _spots in fov.unique_spots_list:
    kept_spots_list.append(_spots[_spots[:,0] > intensity_th])
# finalize candidate spots
cand_chr_spots_list = [[] for _ct in fov.chrom_coords]
for _spots in kept_spots_list:
    _cands_list = assign_spots_to_chromosomes(_spots, fov.chrom_coords)
    for _i, _cands in enumerate(_cands_list):
        cand_chr_spots_list[_i].append(_cands)
print(f"kept chromosomes: {len(fov.chrom_coords)}")
reload(spot_tools.picking)
from ImageAnalysis3.spot_tools.picking import convert_spots_to_hzxys

dna_cand_hzxys_list = [convert_spots_to_hzxys(_spots, fov.shared_parameters['distance_zxy'])
                       for _spots in cand_chr_spots_list]
dna_reg_ids = fov.unique_ids
# select_hzxys close to the chromosome center
dist_th = 3000 # upper limit is 5000nm
sel_dna_cand_hzxys_list = []
for _cand_hzxys, _chrom_coord in zip(dna_cand_hzxys_list, fov.chrom_coords):
    _sel_cands_list = []
    
    for _cands in _cand_hzxys:
        if len(_cands) == 0:
            _sel_cands_list.append([])
        else:
            _dists = np.linalg.norm(_cands[:,1:4] - _chrom_coord*np.array([200,108,108]), axis=1)
            _sel_cands_list.append(_cands[_dists < dist_th])
            
    # append
    sel_dna_cand_hzxys_list.append(_sel_cands_list)<jupyter_output><empty_output><jupyter_text>EM pick spots<jupyter_code>reload(ia.spot_tools.picking)
# load functions
from ImageAnalysis3.spot_tools.picking import Pick_spots_by_intensity, EM_pick_scores_in_population, generate_reference_from_population,evaluate_differences

%matplotlib inline
niter= 10
nkeep = len(sel_dna_cand_hzxys_list)
num_threads = 12
# initialize
init_dna_hzxys = Pick_spots_by_intensity(sel_dna_cand_hzxys_list[:nkeep])
# set save list
sel_dna_hzxys_list, sel_dna_scores_list, all_dna_scores_list = [init_dna_hzxys], [], []
for _iter in range(niter):
    print(f"- iter:{_iter}")
    # generate reference
    ref_ct_dists, ref_local_dists, ref_ints = generate_reference_from_population(
        sel_dna_hzxys_list[-1], dna_reg_ids, 
        sel_dna_hzxys_list[-1][:nkeep], dna_reg_ids,
        num_threads=num_threads,
        collapse_regions=True,
    )
    plt.figure(figsize=(4,2))
    plt.hist(np.ravel(ref_ints), bins=np.arange(0,20,0.5))
    plt.figure(figsize=(4,2))
    plt.hist(np.ravel(ref_ct_dists), bins=np.arange(0,5000,100))
    plt.figure(figsize=(4,2))
    plt.hist(np.ravel(ref_local_dists), bins=np.arange(0,5000,100))
    plt.show()
    # scoring
    sel_hzxys, sel_scores, all_scores = EM_pick_scores_in_population(
            sel_dna_cand_hzxys_list[:nkeep], dna_reg_ids, sel_dna_hzxys_list[-1], 
            ref_ct_dists, ref_local_dists, ref_ints,
            sel_dna_hzxys_list[-1], dna_reg_ids, num_threads=num_threads,
        )
    update_rate = evaluate_differences(sel_hzxys, sel_dna_hzxys_list[-1])
    print(f"-- region kept: {update_rate:.4f}")
    sel_dna_hzxys_list.append(sel_hzxys)
    sel_dna_scores_list.append(sel_scores)
    all_dna_scores_list.append(all_scores)
    
    if update_rate > 0.995:
        break
plt.figure()
plt.hist(np.log(sel_dna_scores_list[-1][5]), 40)
plt.show()
sel_iter = -1

final_dna_hzxys_list = []
distmap_list = []
score_th = np.exp(-8)
bad_spot_percentage = 0.6
for _hzxys, _scores in zip(sel_dna_hzxys_list[sel_iter], sel_dna_scores_list[sel_iter]):
    _kept_hzxys = np.array(_hzxys).copy()
    _kept_hzxys[_scores < score_th] = np.nan
    if np.mean(np.isnan(_kept_hzxys).sum(1)>0)<bad_spot_percentage:
        final_dna_hzxys_list.append(_kept_hzxys)
        distmap_list.append(squareform(pdist(_kept_hzxys[:,1:4])))
from scipy.spatial.distance import pdist, squareform
distmap_list = np.array(distmap_list)
median_distmap = np.nanmedian(distmap_list, axis=0)
loss_rates = np.mean(np.sum(np.isnan(final_dna_hzxys_list), axis=2)>0, axis=0)
print(np.mean(loss_rates))
fig, ax = plt.subplots(figsize=(4,2),dpi=200)
ax.plot(loss_rates, '.-')
ax.set_xticks(np.arange(0,len(fov.unique_ids),20))
plt.show()
kept_inds = np.where(loss_rates<0.3)[0]
cy7_zs = np.array(final_dna_hzxys_list)[:,0::3,1]
cy5_zs = np.array(final_dna_hzxys_list)[:,1::3,1]
cy3_zs = np.array(final_dna_hzxys_list)[:,2::3,1]

plt.figure(dpi=100)
plt.hist(np.ravel(cy7_zs), bins=np.arange(0,6000, 200), 
         alpha=0.5, color='r', label='750')
plt.hist(np.ravel(cy5_zs), bins=np.arange(0,6000, 200), 
         alpha=0.5, color='y', label='647')
plt.hist(np.ravel(cy3_zs), bins=np.arange(0,6000, 200), 
         alpha=0.5, color='g', label='561')
plt.legend()
plt.show()
cy7_ints = np.array(final_dna_hzxys_list)[:,fd00:c2b6:b24b:be67:2827:688d:e6a1:6a3b,0]
cy5_ints = np.array(final_dna_hzxys_list)[:,fdf8:f53e:61e4::18,0]
cy3_ints = np.array(final_dna_hzxys_list)[:,fc00:db20:35b:7399::5,0]

plt.figure(dpi=100)
plt.hist(np.ravel(cy7_ints), bins=np.arange(0,20, 0.5), 
         alpha=0.5, color='r', label='750')
plt.hist(np.ravel(cy5_ints), bins=np.arange(0,20, 0.5), 
         alpha=0.5, color='y', label='647')
plt.hist(np.ravel(cy3_ints), bins=np.arange(0,20, 0.5), 
         alpha=0.5, color='g', label='561')
plt.legend()
plt.show()
fig, ax = plt.subplots(figsize=(4,3),dpi=200)
ax = ia.figure_tools.distmap.plot_distance_map(#median_distmap,
                                               median_distmap[2::3,2::3], 
                                               color_limits=[0,400],
                                               ax=ax,
                                               ticks=np.arange(0,len(fov.unique_ids),20), 
                                               figure_dpi=200)
ax.set_title(f"SI13-5kb IMR90, n={len(distmap_list)}", fontsize=7.5)
plt.gcf().subplots_adjust(bottom=0.1)
plt.show()
fig, ax = plt.subplots(figsize=(4,3),dpi=200)
ax = ia.figure_tools.distmap.plot_distance_map(#median_distmap,
                                               median_distmap[kept_inds][:,kept_inds], 
                                               color_limits=[0,350],
                                               ax=ax,
                                               ticks=np.arange(0,len(fov.unique_ids),20), 
                                               figure_dpi=200)
ax.set_title(f"SI13-5kb IMR90, n={len(distmap_list)}", fontsize=7.5)
plt.gcf().subplots_adjust(bottom=0.1)
plt.show()<jupyter_output><empty_output><jupyter_text>visualize single example<jupyter_code>%matplotlib inline
chrom_id = 0

valid_inds = np.where(np.isnan(final_dna_hzxys_list[chrom_id]).sum(1) == 0)[0]

fig, ax = plt.subplots(figsize=(4,3),dpi=200)
ax = ia.figure_tools.distmap.plot_distance_map(#distmap_list[chrom_id],
                                               distmap_list[chrom_id][valid_inds][:,valid_inds], 
                                               color_limits=[0,400],
                                               ax=ax,
                                               ticks=np.arange(0,150,20), 
                                               figure_dpi=200)
ax.set_title(f"proB bone marrow IgH+/+ chrom: {chrom_id}", fontsize=7.5)
plt.gcf().subplots_adjust(bottom=0.1)
plt.show()



reload(figure_tools.image)
ax3d = figure_tools.image.chromosome_structure_3d_rendering(#final_dna_hzxys_list[chrom_id][:,1:], 
                                                     final_dna_hzxys_list[chrom_id][valid_inds, 1:], 
                                                     marker_edge_line_width=0,
                                                     reference_bar_length=100, image_radius=200, 
                                                     line_width=0.5, figure_dpi=300, depthshade=True)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11440817672271678, 'fraction_numerical': 0.059281239696669964, 'mean_word_length': 2.6712660372791093, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 73, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '503493', 'license': 'MIT', 'n_tokens_mistral': 5917, 'n_tokens_neox': 5159, 'n_words': 870, 'path': '5kb_DNA_analysis/20200707_single_fov_SI13-5kb_IMR90_3color_Crick.ipynb', 'repo_name': 'shiwei23/Chromatin_Analysis_Scripts', 'size': '15099'}