<jupyter_start><jupyter_text>just trying to take a glance at Hi-C biases associated with allele1/2 to try to understand 2x difference in coverage<jupyter_code>%%javascript
require(["codemirror/keymap/sublime", "notebook/js/cell", "base/js/namespace"],
    function(sublime_keymap, cell, IPython) {
        cell.Cell.options_default.cm_config.keyMap = 'sublime';    
        var cells = IPython.notebook.get_cells();    
        for(var cl=0; cl< cells.length ; cl++){
            cells[cl].code_mirror.setOption('keyMap', 'sublime');
        }
    }
);
%load_ext autoreload
%autoreload 2
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib import gridspec
from matplotlib import colors

import os
import re
import numpy as np
import pandas as pd
import cooler

import itertools

import seaborn as sns
# cooltools used dekkerlab-fork, develop-branch, commit-af508bba9c6577b8271b716ba27c330bead15981
from cooler.tools import split, partition
from functools import partial
from cooltools import expected

import bioframe
from bioframe.io import resources
from bioframe import io

from cooltools import numutils
from cooltools.eigdecomp import cooler_cis_eig
from helper_module import fillcolor_compartment_style, to_uscs, bar_signal, random_test
cooler_names = [
    'all-allele-spe-2percent-highcis.galGal5.mapq_30.1000.mcool',
#     'all-allele-spe-1percent-highcis.galGal5.mapq_30.1000.mcool',
    'MNP-DT40-1-3-16-2p-R1-T1__galGal5.galGal5.mapq_30.1000.mcool',
#     'MNP-DT40-1-3-16-R2-T1__galGal5.galGal5.mapq_30.1000.mcool',
    'MNP-DT40-1-3-17-2p-R1-T1__galGal5.galGal5.mapq_30.1000.mcool',
#     'MNP-DT40-1-3-17-R1-T1__galGal5.galGal5.mapq_30.1000.mcool',
    'MNP-DT40-WT1-R1-T1__galGal5.galGal5.mapq_30.1000.mcool',
    'MNP-DT40-WT2-R1-T1__galGal5.galGal5.mapq_30.1000.mcool']<jupyter_output><empty_output><jupyter_text>generate some nicer sample names from coolers ...<jupyter_code># multiqc learnings here, btw!...
renamings = {"all-allele-spe":"pool",
            ".galGal5.mapq_30.1000.mcool":"",
            "MNP-DT40-":"",
            "-T1__galGal5":"",
            "-R1":"",
            "1-3-":"clone-"}

sample_names = []
for c in cooler_names:
    for rf,rt in renamings.items():
        c = re.sub(rf,rt,c)
    sample_names.append(c)

# sample_names

# form a df of samples here:
samples = pd.DataFrame({"sample":sample_names,"mcool":cooler_names})

samples<jupyter_output><empty_output><jupyter_text>read some coolers and extract raw number of interactions between allele1/2 and other chromosomes ...<jupyter_code># READ SOME COOLERS ...
# oldref-20200617
cool_path = "/home/venevs/DOTS_TESTING/alv-chicken/oldref-20200617"
get_cpath = lambda c,b: os.path.join(cool_path,c+f"::/resolutions/{b}")

insertion_site = 72_500_000

dfs = {}
for s in samples["sample"]:
    _bsize = 1_000_000
    s_cname, = samples["mcool"][samples["sample"]==s]
    _c = cooler.Cooler(get_cpath(s_cname,_bsize))
    _df = pd.DataFrame(_c.chromsizes).reset_index().rename({"name":"chrom"},axis=1)
    a1_idx = _c.offset("allele1")
    a2_idx = _c.offset("allele2")
    chr1_idx = slice(*_c.extent("chr1"))
    # chr1_idx = _c.offset(("chr1",insertion_site,insertion_site+_bsize))
    _df["allele1"] = [np.sum(_c.matrix(balance=False)[a1_idx,slice(*_c.extent(chrom))]) for chrom in _c.chromnames]
    _df["allele2"] = [np.sum(_c.matrix(balance=False)[a2_idx,slice(*_c.extent(chrom))]) for chrom in _c.chromnames]
    _df["chr1"] = [np.sum(_c.matrix(balance=False)[chr1_idx,slice(*_c.extent(chrom))]) for chrom in _c.chromnames]
    dfs[s] = _df<jupyter_output><empty_output><jupyter_text>test GC content of allele1 and allele2, and "density" of DpnII cut-sites<jupyter_code>gc_dict = {}
from Bio import SeqIO, SeqUtils
for seq in SeqIO.parse("galGal5.fa","fasta"):
    if seq.id in ["allele1","allele2"]:
        print(seq.id, SeqUtils.GC(seq.seq), len(seq))
        cuts = len(SeqUtils.nt_search(str(seq.seq),"GATC"))-1
        print(seq.id, "cutsites", cuts/len(seq) )
        gc_dict[seq.id] = {"GC":SeqUtils.GC(seq.seq),"length":len(seq),"cut":cuts/len(seq)}
    elif seq.id == "chr1":
        # print(seq.id, SeqUtils.GC(seq.seq[insertion_site-1_000_000:insertion_site+1_000_000]), len(seq))
        # only an estimate of cut sites in chr1 - only in the first 1Mb
        # otherwise it takes forever to calculate ...
        cuts = len(SeqUtils.nt_search(str(seq.seq)[:1_000_000],"GATC"))-1
        print(seq.id, SeqUtils.GC(seq.seq), len(seq))
        print(seq.id, "cutsites", cuts/1_000_000 )
        gc_dict[seq.id] = {"GC":SeqUtils.GC(seq.seq),"length":len(seq),"cut":cuts/1_000_000}
    else:
        pass
# allele1 is a bit more GC-rich
s = samples["sample"][0]
_bsize = 50_000
s_cname, = samples["mcool"][samples["sample"]==s]
clr = cooler.Cooler(get_cpath(s_cname,_bsize))
print(s,clr.info)
binss = []
for seqrec in SeqIO.parse("galGal5.fa","fasta"):
    chrom = seqrec.id
    if chrom in clr.chromnames[:10]:
        bins = clr.bins().fetch(chrom)
        bins["GC"] = [ SeqUtils.GC(seqrec.seq[r.start:r.end])  for r in bins.itertuples()]
        bins["cut"] = [ len(SeqUtils.nt_search(str(seqrec.seq[r.start:r.end]),"GATC"))/(r.end-r.start) for r in bins.itertuples()]
        binss.append(bins) 
binss = pd.concat(binss)
# sns.jointplot(data=binss,y="weight",x="GC",kind="hex", color="k");
g = sns.PairGrid(binss[["GC","weight","cut"]], diag_sharey=False)
g.map_upper(sns.scatterplot,**{"alpha":0.1})
g.map_lower(plt.hexbin ,gridsize=1000)#, **{"kind":"hex"})
g.map_diag(plt.hist, bins=50)<jupyter_output><empty_output><jupyter_text>let's estimate fraction of `cis` interactions for allele1/2this is a rough estimate , calculated by summing raw allele1/2 interactions with chr1 and dividing by the total number of interactions detected with allele1/2<jupyter_code># construct a little table to characterize each sample allele1/allele2 interactions
df_dict = {}
for s in samples["sample"]:
    if "WT" not in s:
        df_dict[s] = {}
        df_dict[s]["allele1"] = {}
        df_dict[s]["allele2"] = {}
        df_dict[s]["chr1"] = {}

df_combo = []
for s in samples["sample"]:
    if "WT" not in s:
        dd = dfs[s]
        for ccc in ["allele1","allele2","chr1"]:
            df_dict[s][ccc]["total interactions"] = str(dd[ccc].sum())
            df_dict[s][ccc]["interactions w. chr1"] = str(dd[dd["chrom"]=="chr1"][ccc].sum())
            df_dict[s][ccc]["cis_percent"] = str(f'{100*dd[dd["chrom"].isin([ccc,"chr1"])][ccc].sum()/dd[ccc].sum():.2f}')
            df_dict[s][ccc]["GC"] = gc_dict[ccc]["GC"]
            df_dict[s][ccc]["DpnII sites/bp"] = gc_dict[ccc]["cut"]
            df_dict[s][ccc]["length"] = gc_dict[ccc]["length"]
            #
            _df = pd.DataFrame(df_dict[s]).T
            _df["sample"] = s
            df_combo.append(_df)
            
        
# allele1 seem to have more inter-chromosomal interactions
# allele2 have average level of interactions with other chromosomes - roughly as a whole chr1

dtypes = {"total interactions":int,
        "interactions w. chr1":int,
        "cis_percent":float,
        "GC":float,
        "length":int,
        }
# nice summary table about the samples/alleles :
final_df = pd.concat(df_combo) \
    .dropna() \
    .drop_duplicates() \
    .reset_index() \
    .rename({"index":"name"},axis=1) \
    .sort_values("name") \
    .reset_index(drop=True) \
    .astype(dtypes)
cols = ["name",
"length",
"GC",
"DpnII sites/bp",
"total interactions",
# "interactions w. chr1",
"cis_percent"]


print("clone-16-2p")
display( final_df[final_df["sample"] == "clone-16-2p"].reset_index(drop=True)[cols] )
print("clone-17-2p")
display( final_df[final_df["sample"] == "clone-17-2p"].reset_index(drop=True)[cols] )
cols = ["name",
"length",
"GC",
"DpnII sites/bp", "sample",
"total interactions",
# "interactions w. chr1",
"cis_percent"]


print("clone-16-2p")
display( final_df[final_df["sample"] == "clone-16-2p"].reset_index(drop=True)[cols] )
print("clone-17-2p")
display( final_df[final_df["sample"] == "clone-17-2p"].reset_index(drop=True)[cols] )
final_df[final_df["sample"] != "pool-2percent-highcis"][cols].to_csv("xxx.csv")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1668896727891984, 'fraction_numerical': 0.03855978591412237, 'mean_word_length': 4.103662321539416, 'pattern_counts': {'":': 25, '<': 22, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '631891', 'license': 'MIT', 'n_tokens_mistral': 3383, 'n_tokens_neox': 3137, 'n_words': 545, 'path': 'explore_bias.ipynb', 'repo_name': 'dekkerlab/allelic-rt-supplement', 'size': '8221'}