# -*- coding: utf-8 -*-
"""hicstraw.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1V4CRdM_hOt4KcM7jBFWjQ6NQyWxECim3

# 1. Setup and Metadata
"""

import numpy as np
import plotly.graph_objs as go
from plotly.offline import iplot

"""
Straw module
Straw enables programmatic access to .hic files.
.hic files store the contact matrices from Hi-C experiments and the
normalization and expected vectors, along with meta-data in the header.
Usage: strawObj = straw <hicFile(s)>
       matrixObj = strawObj.getNormalizedMatrix <chr1> <chr2> <NONE/VC/VC_SQRT/KR> <BP/FRAG> <binsize>
       data = matrixObj.getDataFromBinRegion <x1,x2,y1,y2>
Example:
   import straw
   strawObj = straw(filename)
   matrixObj = strawObj.getNormalizedMatrix('5', '5', 'KR', 'BP', 5000)
   result = matrixObj.getDataFromBinRegion(0,500,0,500)
   for i in range(len(result[0])):
...   print("{0}\t{1}\t{2}".format(result[0][i], result[1][i], result[2][i]))
See https://github.com/theaidenlab/straw/wiki/Python for more documentation
"""
from __future__ import absolute_import, division, print_function, unicode_literals

__author__ = "<NAME>, <NAME>, <NAME>, <NAME>, <NAME>"
__license__ = "MIT"

import struct
import zlib
import requests
import io
import concurrent.futures
import math

def __readcstr(f):
    """ Helper function for reading in C-style string from file
    """
    buf = b""
    while True:
        b = f.read(1)
        if b is None or b == b"\0":
            return buf.decode("utf-8")
        elif b == "":
            raise EOFError("Buffer unexpectedly empty while trying to read null-terminated string")
        else:
            buf += b


"""
functions for chrom.sizes
internal representation is a dictionary with 
chromosome name as the key
value maps to a tuple containing the index and chromosome length
"""


class ChromDotSizes:
    def __init__(self, data):
        self.data = data

    def getLength(self, chrom):
        try:
            return int(self.data[chrom][1])
        except:
            raise ValueError(str(chrom) + " not in chrom.sizes. Check that the chromosome name matches the genome.\n")

    def getIndex(self, chrom):
        try:
            return int(self.data[chrom][0])
        except:
            raise ValueError(str(chrom) + " not in chrom.sizes. Check that the chromosome name matches the genome.\n")

    def figureOutEndpoints(self, chrAndPositions):
        chrAndPositionsArray = chrAndPositions.split(":")
        chrom = chrAndPositionsArray[0]

        indx1 = 0
        indx2 = self.getLength(chrom)

        if len(chrAndPositionsArray) == 3:
            indx1 = int(chrAndPositionsArray[1])
            indx2 = int(chrAndPositionsArray[2])

        return chrom, indx1, indx2
def read_metadata(infile,verbose=False):
    """
    Reads the metadata of HiC file from header.
    Args
    infile: str, path to the HiC file 
    verbose: bool
    
    Returns
    metadata: dict, containing the metadata. 
                Keys of the metadata: 
                HiC version, 
                Master index, 
                Genome ID (str), 
                Attribute dictionary (dict), 
                Chromosomes (dict), 
                Base pair-delimited resolutions (list), 
                Fragment-delimited resolutions (list). 
    """
    metadata={}
    import io
    import struct
    if (infile.startswith("http")):
        # try URL first. 100K should be sufficient for header
        headers={'range' : 'bytes=0-100000', 'x-amz-meta-requester' : 'straw'}
        s = requests.Session()
        r=s.get(infile, headers=headers)
        if (r.status_code >=400):
            print("Error accessing " + infile) 
            print("HTTP status code " + str(r.status_code))
            sys.exit(1)
        req=io.BytesIO(r.content)        
        myrange=r.headers['content-range'].split('/')
        totalbytes=myrange[1]
    else:
        req=open(infile, 'rb')
    magic_string = struct.unpack('<3s', req.read(3))[0]
    req.read(1)
    if (magic_string != b"HIC"):
        sys.exit('This does not appear to be a HiC file magic string is incorrect')
    version = struct.unpack('<i',req.read(4))[0]
    metadata['HiC version']=version 
    masterindex = struct.unpack('<q',req.read(8))[0]
    metadata['Master index']=masterindex
    genome = ""
    c=req.read(1).decode("utf-8") 
    while (c != '\0'):
        genome += c
        c=req.read(1).decode("utf-8") 
    metadata['Genome ID']=genome        
    if (version > 8):
        nvi = struct.unpack('<q',req.read(8))[0]
        nvisize = struct.unpack('<q',req.read(8))[0]
        metadata['NVI'] = nvi
        metadata['NVI size'] = nvisize
    ## read and throw away attribute dictionary (stats+graphs)
    nattributes = struct.unpack('<i',req.read(4))[0]
    d={}
    for x in range(0, nattributes):
        key = __readcstr(req)
        value = __readcstr(req)
        d[key]=value
    metadata['Attribute dictionary']=d
    nChrs = struct.unpack('<i',req.read(4))[0]
    d={}
    for x in range(0, nChrs):
        key = __readcstr(req)
        if (version > 8):
            value = struct.unpack('q',req.read(8))[0]
        else:
            value = struct.unpack('<i',req.read(4))[0]
        d[key]=value
    metadata["Chromosomes"]=d
    nBpRes = struct.unpack('<i',req.read(4))[0]
    l=[]
    for x in range(0, nBpRes):
        res = struct.unpack('<i',req.read(4))[0]
        l.append(res)
    metadata["Base pair-delimited resolutions"]=l        
    nFrag = struct.unpack('<i',req.read(4))[0]
    l=[]
    for x in range(0, nFrag):
        res = struct.unpack('<i',req.read(4))[0]
        l.append(res)
    metadata["Fragment-delimited resolutions"]=l 
    for k in metadata:
        if k!='Attribute dictionary':
            print(k,':',metadata[k])
    if verbose:
        print('Attribute dictionary',':',metadata['Attribute dictionary'])        
    return metadata

# Downloads Hi-C file
#!wget -cq https://www.dropbox.com/s/t3d3kmoerm54dlr/GSM1551620_HIC071.hic

# Define files and HiC read settings
path = '/content/'
hicname = 'GSM1551620_HIC071.hic'
hic_metadata = read_metadata(hicname)
chromosomes = list(hic_metadata["Chromosomes"].keys())
base_resolutions = hic_metadata["Base pair-delimited resolutions"]
fragment_resolutions = hic_metadata["Fragment-delimited resolutions"]

"""# 2. Hi-C Data Processing"""

# Commented out IPython magic to ensure Python compatibility.
# %%capture
# !pip install hic-straw
# import straw

# Saves all self-contact data in Hi-C file
def self_chromosome_interactions(chromosomes, hicname, res):
  genome = {}
  for i in chromosomes: 
    # Adds chromosome interaction to genome dictionary, 
    # interactions stored as list of contactRecord objects
    try:
      genome[i] = straw.straw('observed','KR', hicname, i, i, 'BP', res)
    except: 
      print(f"Interaction data for chromosome {i} not found")
  return genome

# Converts chromosome interactions from contactRecords into a list
def contactRecords_to_list(records):
  table = np.zeros((len(records), 3))
  index = 0 
  for contact in records:
    table[index][0] = int(contact.binX)
    table[index][1] = int(contact.binY)
    table[index][2] = contact.counts
    index += 1
  return table

# Converts data in resolution units into particle numbers
def res_correct(table, res):
    table[:, 0] //= res
    table[:, 1] //= res
    table[:, 0] += 1
    table[:, 1] += 1
    return table

# Keeps entries where contact counts are higher than user-given threshold
# and gets rid of atom self contacts
def lammps_correct(table, threshold):
    table = table[table[:,2] > threshold]
    table = table[table[:,0] != table[:,1]]
    return table

res = 250000
threshold = 3.3

#genome = self_chromosome_interactions(chromosomes, hicname, res)
genome = {}
for i in chromosomes: 
    # Adds chromosome interaction to genome dictionary, 
    # interactions stored as list of contactRecord objects
    try:
      genome[i] = straw.straw('observed','KR', hicname, i, i, 'BP', res)
    except: 
      print(f"Interaction data for chromosome {i} not found")

x_correct = lammps_correct(res_correct(contactRecords_to_list(genome["X"]), res), threshold)

"""# 3. LAMMPS Input & Data File Construction"""

import time

# function to plot 3D structure
def plot_3D(coords):
    x = np.zeros(len(coords))
    y = np.zeros(len(coords))
    z = np.zeros(len(coords))
    for i in range(len(coords)):
        x[i] = coords[i][0]
        y[i] = coords[i][1]
        z[i] = coords[i][2]
    trace = go.Scatter3d(
      x = x, y = y, z = z, mode = 'lines+markers', marker = dict(
          size = 5,
          colorscale = 'Viridis'
          )
      )
    layout = go.Layout(title = f'Initial Random Structure')
    fig = go.Figure(data = [trace], layout = layout)
    iplot(fig)
    return None

# function to check if next site is already occupied
def check_if_free(lattice_coords, next_coords, index):
    for i in range(index):
        if lattice_coords[i][0] == next_coords[0] and lattice_coords[i][1] == next_coords[1] \
                and lattice_coords[i][2] == next_coords[2]:
            return False
    return True

# function to create random 3D walk on lattice
def random_walk(n):
    backtrack = 10
    lattice_coords = np.zeros([n, 3])
    steps = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],
                      [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
    i = 1
    while i < n:
        issue = False
        s = time.time()
        #if i % 100 == 0:
            #print(i, s-start, 's')
        rand = np.random.randint(0, 6)
        next_coords = lattice_coords[i - 1] + steps[rand]
        while check_if_free(lattice_coords, next_coords, i) == False:
            rand = np.random.randint(0, 6)
            next_coords = lattice_coords[i-1] + steps[rand]
            e = time.time()
            if e - s > 0.1:
                issue = True
                #print('Stuck! Go back and find a new way... %s' % i)
                for k in range(1, backtrack + 1):
                    lattice_coords[i-k] = np.zeros(3)
                i -= backtrack + 1
                break
        if issue == False:
            lattice_coords[i] = next_coords
            i += 1
    return lattice_coords

# function to create molecule tags
def create_molecule_tags(n, lengths):
    tags = []
    tag = 1
    cumlength = np.cumsum(lengths)
    for i in range(1, n+1):
        if i - 1 in cumlength:
            tag += 1
        tags.append(tag)

    return tags

# function to create bonds
def create_bonds(n, lengths):
    bonds = []
    cumlength = np.cumsum(lengths)
    for i in range(1, n):
        if i not in cumlength:
            bonds.append([i, i+1])
    return bonds


# function to create angles
def create_angles(n, lengths):
    angles = []
    cumlength = np.cumsum(lengths)
    for i in range(1, n-1):
        if (i not in cumlength) and (i+1 not in cumlength):
            angles.append([i, i + 1, i + 2])
    return angles

# function to create data file
def create_datafile(n, lengths, spacing, lattice_numbers, dimensions):
    chains = int(len(lengths))  # number of chains
    bond_number = int(sum(lengths) - chains) # number of bonds

    angle_number = 0
    for length in lengths: 
        if length > 2.0:
            angle_number += int(length - 2)  # number of bond angles

    lattice_coords = random_walk(n) * spacing  # coordinates of lattice points
    tags = create_molecule_tags(n, lengths)  # molecule tags
    bonds = create_bonds(n, lengths)  # indicates bonds between particles
    angles = create_angles(n, lengths)  # indicates angles between particles

    # open datafile to write to
    datafile = open(f'random_coil_N{n}.dat', 'w')
    datafile.write(f'LAMMPS data file for random 3D walk on lattice: N = {n}, Chain length = {length}\n\n')
    datafile.write(f'{n} atoms\n1 atom types\n{bond_number} bonds\n2 bond types\n1000 extra bond per atom\n{angle_number} angles\n1 angle types\n\n')
    datafile.write(f'{-dimensions[0] / 2} {dimensions[0] / 2} xlo xhi\n{-dimensions[1] / 2} {dimensions[1] / 2} ylo yhi\n{-dimensions[2] / 2} {dimensions[2] / 2} zlo zhi\n\n')
    datafile.write('Masses\n\n1 1\n\nAtoms\n')
    for i in range(n):
        datafile.write(f'\n{i + 1}\t{tags[i]}\t1\t{lattice_coords[i][0]}\t{lattice_coords[i][1]}\t{lattice_coords[i][2]}\t0\t0\t0')
    if bond_number > 0:
        datafile.write('\n\nBonds\n')
        for i in range(len(bonds)):
            datafile.write(f'\n{i + 1}\t1\t{bonds[i][0]}\t{bonds[i][1]}')
    if angle_number > 0:
        datafile.write('\n\nAngles\n')
        for i in range(len(angles)):
            datafile.write(f'\n{i + 1}\t1\t{angles[i][0]}\t{angles[i][1]}\t{angles[i][2]}')
    datafile.close()
    plot_3D(lattice_coords)
    return None

# function to create input file
def create_inputfile(n, timesteps, bondconnect):
    # opens input file to write to
    datafile = open('in.input', 'w')
    dataname = f'random_coil_N{n}.dat' # data file name
    lang = np.random.randint(1,1000000) # generates noise term for langevin
    
    datafile.write('log sim.log\nunits lj\n\n')
    datafile.write('atom_style angle\nboundary        p p p\n\n')
    datafile.write('neighbor 4 bin\nneigh_modify every 1 delay 1 check yes\n\n')
    datafile.write('atom_modify sort 0 0\n\n')
    datafile.write(f'#restart 1000000 N{n}.restart\n\n')
    datafile.write(f'read_data {dataname}\n')
    datafile.write('reset_timestep 0\n\n')
    datafile.write(f'write_data equilibrated_N{n}.dat\n\n')
    datafile.write('group all type 1\n\n')
    datafile.write('dump   1   all   custom   1000   sim.dump  id  x y z  ix iy iz\n')
    datafile.write('dump_modify   1   format line "%d %.5f %.5f %.5f %d %d %d"\n\n')
    datafile.write('angle_style   cosine\nangle_coeff   1 0.0\n\n')
    datafile.write('pair_style      lj/cut 1.12246152962189\n')
    datafile.write('pair_modify     shift yes\n')
    datafile.write('pair_coeff      * * 1.0 1.0\n\n')
    datafile.write('bond_style hybrid harmonic fene\n')
    datafile.write('bond_coeff 1 fene 30.0 10 1.0 1.0\n')
    datafile.write('bond_coeff 2 harmonic  1.0 2.2\n')
    datafile.write('special_bonds fene\n\n')
    datafile.write('fix 1 all nve\n')
    datafile.write(f'fix 2 all langevin   1.0 1.0   1.0   {lang}\n\n')
    datafile.write('thermo 50000\n\n')
    
    for i in range(len(bondconnect)):
        datafile.write(f'create_bonds single/bond 2 {int(bondconnect[i][0])} {int(bondconnect[i][1])} special yes\n')
    
    datafile.write('thermo_style   custom   step  temp etotal epair  emol  press pxx pyy pzz lx ly lz pe ke ebond evdwl\n\n')
    datafile.write('timestep 0.00001\n')
    datafile.write(f'run {timesteps}')

# Defining LAMMPS properties
n = int(max(x_correct.max(axis = 0)[0], x_correct.max(axis = 0)[1]))  # total number of particles
lengths = [n]  # length of chains
spacing = 3.0  # lattice spacing
lattice_numbers = np.array([200, 200, 200])
dimensions = lattice_numbers * 2  # dimensions of box
timesteps = 1000000

# Creating input and data files for simulation
start = time.time()
create_datafile(n, lengths, spacing, lattice_numbers, dimensions)
create_inputfile(n, timesteps, x_correct)
end = time.time()

"""# 4. LAMMPS Simulation"""

# Commented out IPython magic to ensure Python compatibility.
# %%capture
# # Creates LAMMPS executable
# !wget https://github.com/lmhale99/atomman-demo/raw/main/lmp.gz
# !gunzip lmp.gz
# !chmod 755 lmp
# 
# import subprocess
# import shlex

# Runs LAMMPS simulation based on input file
def run_lammps(infile):
  command = shlex.split(f'./lmp -in {infile}')
  subprocess.run(command, check=True, capture_output=False, text=True)

# Reads dump file and extracts coordinate data
def readdump(file):
  f = open(file).read().splitlines()
  dump = {}
  
  #finds coordinates using "ITEM:"
  #creates a dictionary with timestep and coordinates
  i = 0
  while i < len(f): 
    timestep = int(f[f.index("ITEM: TIMESTEP", i) + 1].strip())
    
    begind = f.index("ITEM: ATOMS id x y z ix iy iz", i)
    try:
      endind = f.index("ITEM: TIMESTEP", i+1)
    except ValueError:
      endind = len(f)
    coords = f[begind + 1:endind]
  
    data = np.zeros((len(coords), 7))
    index = 0
    for j in coords:
      line = j.split()
      data[index][0] = int(line[0])
      data[index][1] = float(line[1])
      data[index][2] = float(line[2])
      data[index][3] = float(line[3])
      data[index][4] = float(line[4])
      data[index][5] = float(line[5])
      data[index][6] = float(line[6])
      index += 1

    dump[timestep] = data
    i = endind
  
  return dump

# Displays simulation at timestep
def dumpplot(dump, timestep):
  trace = go.Scatter3d(
    x = dump[timestep][:,1], y = dump[timestep][:,2], z = dump[timestep][:,3], mode = 'lines+markers', marker = dict(
        size = 5,
        color = dump[timestep][:,0],
        colorscale = 'Viridis'
        )
    )
  layout = go.Layout(title = f'Simulation at timestep {timestep}')
  fig = go.Figure(data = [trace], layout = layout)
  iplot(fig)

run_lammps("in.input")

dump = readdump("sim.dump")
dumpplot(dump, 100000)

"""# 5. Analysis"""

#!wget -cq https://www.dropbox.com/s/nkbuklgq2c9ewvw/nWTXa
from scipy.stats import pearsonr
from scipy.spatial.distance import pdist, squareform

#Finds all contacts
def find_contacts(data, timestep, dist):
    coords = data[timestep][:,1:4]
    IDs = data[timestep][:,0]
    distances = np.triu(squareform(pdist(coords)))
    contact_IDarrays = np.where((distances<dist) & (distances>0))
    C = np.ma.size(contact_IDarrays[0])
    contact_IDs = np.hstack((np.reshape(np.take(IDs,contact_IDarrays[0]), (C,1)), np.reshape(np.take(IDs,contact_IDarrays[1]), (C,1))))   
    contacts = np.hstack((np.reshape(contact_IDs[:,0], (C,1)), np.reshape(contact_IDs[:,1], (C,1)), distances[np.reshape(contact_IDarrays[0], (C,1)), np.reshape(contact_IDarrays[1], (C,1))])) 
    return contacts

def find_nwtxa_contacts(file, dist):
    rawdata = open(file).read().splitlines()
    table = np.zeros((len(rawdata), 3))
    index = 0 
    for i in rawdata:
      line = i.split()
      if(float(line[0]) > float(line[1]) and line[2] != 'NA' and float(line[2]) > dist):
        table[index][0] = int(line[0])
        table[index][1] = int(line[1])
        table[index][2] = float(line[2])
        index += 1
    return table[0:index]

def plotcontactmap(contacts):
  fig = go.Figure(data=[go.Scatter(
    x = contacts[:,0],
    y = contacts[:,1],
    mode = 'markers',
    marker=dict(
        opacity=0.9,
        color='Black',
        size=1,
    ))
    ])

  fig.update_layout(
      width = 600,
      height = 600,
      title = "Contact Map"
  )

  fig.update_yaxes(
      scaleanchor = "x",
      scaleratio = 1,
    )

  fig.show()

def combine(mapone, maptwo):
  onemax = int(max(mapone.max(axis = 0)[0], mapone.max(axis = 0)[1]))
  twomax = int(max(maptwo.max(axis = 0)[0], maptwo.max(axis = 0)[1]))
  maptwo[:, 0] *= onemax/twomax
  maptwo[:, 1] *= onemax/twomax
  combined = np.vstack((mapone, maptwo))
  return combined

contacts = find_contacts(dump, 100000, 3.3)
nwtxa = find_nwtxa_contacts("nWTXa", 5)
combined = combine(nwtxa, contacts)
plotcontactmap(combined)
# -*- coding: utf-8 -*-
"""hicstraw.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1V4CRdM_hOt4KcM7jBFWjQ6NQyWxECim3

# 1. Setup and Metadata
"""

import numpy as np
import plotly.graph_objs as go
from plotly.offline import iplot

"""
Straw module
Straw enables programmatic access to .hic files.
.hic files store the contact matrices from Hi-C experiments and the
normalization and expected vectors, along with meta-data in the header.
Usage: strawObj = straw <hicFile(s)>
       matrixObj = strawObj.getNormalizedMatrix <chr1> <chr2> <NONE/VC/VC_SQRT/KR> <BP/FRAG> <binsize>
       data = matrixObj.getDataFromBinRegion <x1,x2,y1,y2>
Example:
   import straw
   strawObj = straw(filename)
   matrixObj = strawObj.getNormalizedMatrix('5', '5', 'KR', 'BP', 5000)
   result = matrixObj.getDataFromBinRegion(0,500,0,500)
   for i in range(len(result[0])):
...   print("{0}\t{1}\t{2}".format(result[0][i], result[1][i], result[2][i]))
See https://github.com/theaidenlab/straw/wiki/Python for more documentation
"""
from __future__ import absolute_import, division, print_function, unicode_literals

__author__ = "<NAME>, <NAME>, <NAME>, <NAME>, <NAME>"
__license__ = "MIT"

import struct
import zlib
import requests
import io
import concurrent.futures
import math

def __readcstr(f):
    """ Helper function for reading in C-style string from file
    """
    buf = b""
    while True:
        b = f.read(1)
        if b is None or b == b"\0":
            return buf.decode("utf-8")
        elif b == "":
            raise EOFError("Buffer unexpectedly empty while trying to read null-terminated string")
        else:
            buf += b


"""
functions for chrom.sizes
internal representation is a dictionary with 
chromosome name as the key
value maps to a tuple containing the index and chromosome length
"""


class ChromDotSizes:
    def __init__(self, data):
        self.data = data

    def getLength(self, chrom):
        try:
            return int(self.data[chrom][1])
        except:
            raise ValueError(str(chrom) + " not in chrom.sizes. Check that the chromosome name matches the genome.\n")

    def getIndex(self, chrom):
        try:
            return int(self.data[chrom][0])
        except:
            raise ValueError(str(chrom) + " not in chrom.sizes. Check that the chromosome name matches the genome.\n")

    def figureOutEndpoints(self, chrAndPositions):
        chrAndPositionsArray = chrAndPositions.split(":")
        chrom = chrAndPositionsArray[0]

        indx1 = 0
        indx2 = self.getLength(chrom)

        if len(chrAndPositionsArray) == 3:
            indx1 = int(chrAndPositionsArray[1])
            indx2 = int(chrAndPositionsArray[2])

        return chrom, indx1, indx2
def read_metadata(infile,verbose=False):
    """
    Reads the metadata of HiC file from header.
    Args
    infile: str, path to the HiC file 
    verbose: bool
    
    Returns
    metadata: dict, containing the metadata. 
                Keys of the metadata: 
                HiC version, 
                Master index, 
                Genome ID (str), 
                Attribute dictionary (dict), 
                Chromosomes (dict), 
                Base pair-delimited resolutions (list), 
                Fragment-delimited resolutions (list). 
    """
    metadata={}
    import io
    import struct
    if (infile.startswith("http")):
        # try URL first. 100K should be sufficient for header
        headers={'range' : 'bytes=0-100000', 'x-amz-meta-requester' : 'straw'}
        s = requests.Session()
        r=s.get(infile, headers=headers)
        if (r.status_code >=400):
            print("Error accessing " + infile) 
            print("HTTP status code " + str(r.status_code))
            sys.exit(1)
        req=io.BytesIO(r.content)        
        myrange=r.headers['content-range'].split('/')
        totalbytes=myrange[1]
    else:
        req=open(infile, 'rb')
    magic_string = struct.unpack('<3s', req.read(3))[0]
    req.read(1)
    if (magic_string != b"HIC"):
        sys.exit('This does not appear to be a HiC file magic string is incorrect')
    version = struct.unpack('<i',req.read(4))[0]
    metadata['HiC version']=version 
    masterindex = struct.unpack('<q',req.read(8))[0]
    metadata['Master index']=masterindex
    genome = ""
    c=req.read(1).decode("utf-8") 
    while (c != '\0'):
        genome += c
        c=req.read(1).decode("utf-8") 
    metadata['Genome ID']=genome        
    if (version > 8):
        nvi = struct.unpack('<q',req.read(8))[0]
        nvisize = struct.unpack('<q',req.read(8))[0]
        metadata['NVI'] = nvi
        metadata['NVI size'] = nvisize
    ## read and throw away attribute dictionary (stats+graphs)
    nattributes = struct.unpack('<i',req.read(4))[0]
    d={}
    for x in range(0, nattributes):
        key = __readcstr(req)
        value = __readcstr(req)
        d[key]=value
    metadata['Attribute dictionary']=d
    nChrs = struct.unpack('<i',req.read(4))[0]
    d={}
    for x in range(0, nChrs):
        key = __readcstr(req)
        if (version > 8):
            value = struct.unpack('q',req.read(8))[0]
        else:
            value = struct.unpack('<i',req.read(4))[0]
        d[key]=value
    metadata["Chromosomes"]=d
    nBpRes = struct.unpack('<i',req.read(4))[0]
    l=[]
    for x in range(0, nBpRes):
        res = struct.unpack('<i',req.read(4))[0]
        l.append(res)
    metadata["Base pair-delimited resolutions"]=l        
    nFrag = struct.unpack('<i',req.read(4))[0]
    l=[]
    for x in range(0, nFrag):
        res = struct.unpack('<i',req.read(4))[0]
        l.append(res)
    metadata["Fragment-delimited resolutions"]=l 
    for k in metadata:
        if k!='Attribute dictionary':
            print(k,':',metadata[k])
    if verbose:
        print('Attribute dictionary',':',metadata['Attribute dictionary'])        
    return metadata

# Downloads Hi-C file
#!wget -cq https://www.dropbox.com/s/t3d3kmoerm54dlr/GSM1551620_HIC071.hic

# Define files and HiC read settings
path = '/content/'
hicname = 'GSM1551620_HIC071.hic'
hic_metadata = read_metadata(hicname)
chromosomes = list(hic_metadata["Chromosomes"].keys())
base_resolutions = hic_metadata["Base pair-delimited resolutions"]
fragment_resolutions = hic_metadata["Fragment-delimited resolutions"]

"""# 2. Hi-C Data Processing"""

# Commented out IPython magic to ensure Python compatibility.
# %%capture
# !pip install hic-straw
# import straw

# Saves all self-contact data in Hi-C file
def self_chromosome_interactions(chromosomes, hicname, res):
  genome = {}
  for i in chromosomes: 
    # Adds chromosome interaction to genome dictionary, 
    # interactions stored as list of contactRecord objects
    try:
      genome[i] = straw.straw('observed','KR', hicname, i, i, 'BP', res)
    except: 
      print(f"Interaction data for chromosome {i} not found")
  return genome

# Converts chromosome interactions from contactRecords into a list
def contactRecords_to_list(records):
  table = np.zeros((len(records), 3))
  index = 0 
  for contact in records:
    table[index][0] = int(contact.binX)
    table[index][1] = int(contact.binY)
    table[index][2] = contact.counts
    index += 1
  return table

# Converts data in resolution units into particle numbers
def res_correct(table, res):
    table[:, 0] //= res
    table[:, 1] //= res
    table[:, 0] += 1
    table[:, 1] += 1
    return table

# Keeps entries where contact counts are higher than user-given threshold
# and gets rid of atom self contacts
def lammps_correct(table, threshold):
    table = table[table[:,2] > threshold]
    table = table[table[:,0] != table[:,1]]
    return table

res = 250000
threshold = 3.3

#genome = self_chromosome_interactions(chromosomes, hicname, res)
genome = {}
for i in chromosomes: 
    # Adds chromosome interaction to genome dictionary, 
    # interactions stored as list of contactRecord objects
    try:
      genome[i] = straw.straw('observed','KR', hicname, i, i, 'BP', res)
    except: 
      print(f"Interaction data for chromosome {i} not found")

x_correct = lammps_correct(res_correct(contactRecords_to_list(genome["X"]), res), threshold)

"""# 3. LAMMPS Input & Data File Construction"""

import time

# function to plot 3D structure
def plot_3D(coords):
    x = np.zeros(len(coords))
    y = np.zeros(len(coords))
    z = np.zeros(len(coords))
    for i in range(len(coords)):
        x[i] = coords[i][0]
        y[i] = coords[i][1]
        z[i] = coords[i][2]
    trace = go.Scatter3d(
      x = x, y = y, z = z, mode = 'lines+markers', marker = dict(
          size = 5,
          colorscale = 'Viridis'
          )
      )
    layout = go.Layout(title = f'Initial Random Structure')
    fig = go.Figure(data = [trace], layout = layout)
    iplot(fig)
    return None

# function to check if next site is already occupied
def check_if_free(lattice_coords, next_coords, index):
    for i in range(index):
        if lattice_coords[i][0] == next_coords[0] and lattice_coords[i][1] == next_coords[1] \
                and lattice_coords[i][2] == next_coords[2]:
            return False
    return True

# function to create random 3D walk on lattice
def random_walk(n):
    backtrack = 10
    lattice_coords = np.zeros([n, 3])
    steps = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],
                      [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
    i = 1
    while i < n:
        issue = False
        s = time.time()
        #if i % 100 == 0:
            #print(i, s-start, 's')
        rand = np.random.randint(0, 6)
        next_coords = lattice_coords[i - 1] + steps[rand]
        while check_if_free(lattice_coords, next_coords, i) == False:
            rand = np.random.randint(0, 6)
            next_coords = lattice_coords[i-1] + steps[rand]
            e = time.time()
            if e - s > 0.1:
                issue = True
                #print('Stuck! Go back and find a new way... %s' % i)
                for k in range(1, backtrack + 1):
                    lattice_coords[i-k] = np.zeros(3)
                i -= backtrack + 1
                break
        if issue == False:
            lattice_coords[i] = next_coords
            i += 1
    return lattice_coords

# function to create molecule tags
def create_molecule_tags(n, lengths):
    tags = []
    tag = 1
    cumlength = np.cumsum(lengths)
    for i in range(1, n+1):
        if i - 1 in cumlength:
            tag += 1
        tags.append(tag)

    return tags

# function to create bonds
def create_bonds(n, lengths):
    bonds = []
    cumlength = np.cumsum(lengths)
    for i in range(1, n):
        if i not in cumlength:
            bonds.append([i, i+1])
    return bonds


# function to create angles
def create_angles(n, lengths):
    angles = []
    cumlength = np.cumsum(lengths)
    for i in range(1, n-1):
        if (i not in cumlength) and (i+1 not in cumlength):
            angles.append([i, i + 1, i + 2])
    return angles

# function to create data file
def create_datafile(n, lengths, spacing, lattice_numbers, dimensions):
    chains = int(len(lengths))  # number of chains
    bond_number = int(sum(lengths) - chains) # number of bonds

    angle_number = 0
    for length in lengths: 
        if length > 2.0:
            angle_number += int(length - 2)  # number of bond angles

    lattice_coords = random_walk(n) * spacing  # coordinates of lattice points
    tags = create_molecule_tags(n, lengths)  # molecule tags
    bonds = create_bonds(n, lengths)  # indicates bonds between particles
    angles = create_angles(n, lengths)  # indicates angles between particles

    # open datafile to write to
    datafile = open(f'random_coil_N{n}.dat', 'w')
    datafile.write(f'LAMMPS data file for random 3D walk on lattice: N = {n}, Chain length = {length}\n\n')
    datafile.write(f'{n} atoms\n1 atom types\n{bond_number} bonds\n2 bond types\n1000 extra bond per atom\n{angle_number} angles\n1 angle types\n\n')
    datafile.write(f'{-dimensions[0] / 2} {dimensions[0] / 2} xlo xhi\n{-dimensions[1] / 2} {dimensions[1] / 2} ylo yhi\n{-dimensions[2] / 2} {dimensions[2] / 2} zlo zhi\n\n')
    datafile.write('Masses\n\n1 1\n\nAtoms\n')
    for i in range(n):
        datafile.write(f'\n{i + 1}\t{tags[i]}\t1\t{lattice_coords[i][0]}\t{lattice_coords[i][1]}\t{lattice_coords[i][2]}\t0\t0\t0')
    if bond_number > 0:
        datafile.write('\n\nBonds\n')
        for i in range(len(bonds)):
            datafile.write(f'\n{i + 1}\t1\t{bonds[i][0]}\t{bonds[i][1]}')
    if angle_number > 0:
        datafile.write('\n\nAngles\n')
        for i in range(len(angles)):
            datafile.write(f'\n{i + 1}\t1\t{angles[i][0]}\t{angles[i][1]}\t{angles[i][2]}')
    datafile.close()
    plot_3D(lattice_coords)
    return None

# function to create input file
def create_inputfile(n, timesteps, bondconnect):
    # opens input file to write to
    datafile = open('in.input', 'w')
    dataname = f'random_coil_N{n}.dat' # data file name
    lang = np.random.randint(1,1000000) # generates noise term for langevin
    
    datafile.write('log sim.log\nunits lj\n\n')
    datafile.write('atom_style angle\nboundary        p p p\n\n')
    datafile.write('neighbor 4 bin\nneigh_modify every 1 delay 1 check yes\n\n')
    datafile.write('atom_modify sort 0 0\n\n')
    datafile.write(f'#restart 1000000 N{n}.restart\n\n')
    datafile.write(f'read_data {dataname}\n')
    datafile.write('reset_timestep 0\n\n')
    datafile.write(f'write_data equilibrated_N{n}.dat\n\n')
    datafile.write('group all type 1\n\n')
    datafile.write('dump   1   all   custom   1000   sim.dump  id  x y z  ix iy iz\n')
    datafile.write('dump_modify   1   format line "%d %.5f %.5f %.5f %d %d %d"\n\n')
    datafile.write('angle_style   cosine\nangle_coeff   1 0.0\n\n')
    datafile.write('pair_style      lj/cut 1.12246152962189\n')
    datafile.write('pair_modify     shift yes\n')
    datafile.write('pair_coeff      * * 1.0 1.0\n\n')
    datafile.write('bond_style hybrid harmonic fene\n')
    datafile.write('bond_coeff 1 fene 30.0 10 1.0 1.0\n')
    datafile.write('bond_coeff 2 harmonic  1.0 2.2\n')
    datafile.write('special_bonds fene\n\n')
    datafile.write('fix 1 all nve\n')
    datafile.write(f'fix 2 all langevin   1.0 1.0   1.0   {lang}\n\n')
    datafile.write('thermo 50000\n\n')
    
    for i in range(len(bondconnect)):
        datafile.write(f'create_bonds single/bond 2 {int(bondconnect[i][0])} {int(bondconnect[i][1])} special yes\n')
    
    datafile.write('thermo_style   custom   step  temp etotal epair  emol  press pxx pyy pzz lx ly lz pe ke ebond evdwl\n\n')
    datafile.write('timestep 0.00001\n')
    datafile.write(f'run {timesteps}')

# Defining LAMMPS properties
n = int(max(x_correct.max(axis = 0)[0], x_correct.max(axis = 0)[1]))  # total number of particles
lengths = [n]  # length of chains
spacing = 3.0  # lattice spacing
lattice_numbers = np.array([200, 200, 200])
dimensions = lattice_numbers * 2  # dimensions of box
timesteps = 1000000

# Creating input and data files for simulation
start = time.time()
create_datafile(n, lengths, spacing, lattice_numbers, dimensions)
create_inputfile(n, timesteps, x_correct)
end = time.time()

"""# 4. LAMMPS Simulation"""

# Commented out IPython magic to ensure Python compatibility.
# %%capture
# # Creates LAMMPS executable
# !wget https://github.com/lmhale99/atomman-demo/raw/main/lmp.gz
# !gunzip lmp.gz
# !chmod 755 lmp
# 
# import subprocess
# import shlex

# Runs LAMMPS simulation based on input file
def run_lammps(infile):
  command = shlex.split(f'./lmp -in {infile}')
  subprocess.run(command, check=True, capture_output=False, text=True)

# Reads dump file and extracts coordinate data
def readdump(file):
  f = open(file).read().splitlines()
  dump = {}
  
  #finds coordinates using "ITEM:"
  #creates a dictionary with timestep and coordinates
  i = 0
  while i < len(f): 
    timestep = int(f[f.index("ITEM: TIMESTEP", i) + 1].strip())
    
    begind = f.index("ITEM: ATOMS id x y z ix iy iz", i)
    try:
      endind = f.index("ITEM: TIMESTEP", i+1)
    except ValueError:
      endind = len(f)
    coords = f[begind + 1:endind]
  
    data = np.zeros((len(coords), 7))
    index = 0
    for j in coords:
      line = j.split()
      data[index][0] = int(line[0])
      data[index][1] = float(line[1])
      data[index][2] = float(line[2])
      data[index][3] = float(line[3])
      data[index][4] = float(line[4])
      data[index][5] = float(line[5])
      data[index][6] = float(line[6])
      index += 1

    dump[timestep] = data
    i = endind
  
  return dump

# Displays simulation at timestep
def dumpplot(dump, timestep):
  trace = go.Scatter3d(
    x = dump[timestep][:,1], y = dump[timestep][:,2], z = dump[timestep][:,3], mode = 'lines+markers', marker = dict(
        size = 5,
        color = dump[timestep][:,0],
        colorscale = 'Viridis'
        )
    )
  layout = go.Layout(title = f'Simulation at timestep {timestep}')
  fig = go.Figure(data = [trace], layout = layout)
  iplot(fig)

run_lammps("in.input")

dump = readdump("sim.dump")
dumpplot(dump, 100000)

"""# 5. Analysis"""

#!wget -cq https://www.dropbox.com/s/nkbuklgq2c9ewvw/nWTXa
from scipy.stats import pearsonr
from scipy.spatial.distance import pdist, squareform

#Finds all contacts
def find_contacts(data, timestep, dist):
    coords = data[timestep][:,1:4]
    IDs = data[timestep][:,0]
    distances = np.triu(squareform(pdist(coords)))
    contact_IDarrays = np.where((distances<dist) & (distances>0))
    C = np.ma.size(contact_IDarrays[0])
    contact_IDs = np.hstack((np.reshape(np.take(IDs,contact_IDarrays[0]), (C,1)), np.reshape(np.take(IDs,contact_IDarrays[1]), (C,1))))   
    contacts = np.hstack((np.reshape(contact_IDs[:,0], (C,1)), np.reshape(contact_IDs[:,1], (C,1)), distances[np.reshape(contact_IDarrays[0], (C,1)), np.reshape(contact_IDarrays[1], (C,1))])) 
    return contacts

def find_nwtxa_contacts(file, dist):
    rawdata = open(file).read().splitlines()
    table = np.zeros((len(rawdata), 3))
    index = 0 
    for i in rawdata:
      line = i.split()
      if(float(line[0]) > float(line[1]) and line[2] != 'NA' and float(line[2]) > dist):
        table[index][0] = int(line[0])
        table[index][1] = int(line[1])
        table[index][2] = float(line[2])
        index += 1
    return table[0:index]

def plotcontactmap(contacts):
  fig = go.Figure(data=[go.Scatter(
    x = contacts[:,0],
    y = contacts[:,1],
    mode = 'markers',
    marker=dict(
        opacity=0.9,
        color='Black',
        size=1,
    ))
    ])

  fig.update_layout(
      width = 600,
      height = 600,
      title = "Contact Map"
  )

  fig.update_yaxes(
      scaleanchor = "x",
      scaleratio = 1,
    )

  fig.show()

def combine(mapone, maptwo):
  onemax = int(max(mapone.max(axis = 0)[0], mapone.max(axis = 0)[1]))
  twomax = int(max(maptwo.max(axis = 0)[0], maptwo.max(axis = 0)[1]))
  maptwo[:, 0] *= onemax/twomax
  maptwo[:, 1] *= onemax/twomax
  combined = np.vstack((mapone, maptwo))
  return combined

contacts = find_contacts(dump, 100000, 3.3)
nwtxa = find_nwtxa_contacts("nWTXa", 5)
combined = combine(nwtxa, contacts)
plotcontactmap(combined)
# -*- coding: utf-8 -*- hicstraw.ipynb  Automatically generated by Colaboratory.  Original file is located at     https://colab.research.google.com/drive/1V4CRdM_hOt4KcM7jBFWjQ6NQyWxECim3  # 1. Setup and Metadata Straw module Straw enables programmatic access to .hic files. .hic files store the contact matrices from Hi-C experiments and the normalization and expected vectors, along with meta-data in the header. Usage: strawObj = straw <hicFile(s)>        matrixObj = strawObj.getNormalizedMatrix <chr1> <chr2> <NONE/VC/VC_SQRT/KR> <BP/FRAG> <binsize>        data = matrixObj.getDataFromBinRegion <x1,x2,y1,y2> Example:    import straw    strawObj = straw(filename)    matrixObj = strawObj.getNormalizedMatrix('5', '5', 'KR', 'BP', 5000)    result = matrixObj.getDataFromBinRegion(0,500,0,500)    for i in range(len(result[0])): ...   print("{0}\t{1}\t{2}".format(result[0][i], result[1][i], result[2][i])) See https://github.com/theaidenlab/straw/wiki/Python for more documentation Helper function for reading in C-style string from file functions for chrom.sizes internal representation is a dictionary with  chromosome name as the key value maps to a tuple containing the index and chromosome length Reads the metadata of HiC file from header.     Args     infile: str, path to the HiC file      verbose: bool          Returns     metadata: dict, containing the metadata.                  Keys of the metadata:                  HiC version,                  Master index,                  Genome ID (str),                  Attribute dictionary (dict),                  Chromosomes (dict),                  Base pair-delimited resolutions (list),                  Fragment-delimited resolutions (list). # try URL first. 100K should be sufficient for header ## read and throw away attribute dictionary (stats+graphs) # Downloads Hi-C file #!wget -cq https://www.dropbox.com/s/t3d3kmoerm54dlr/GSM1551620_HIC071.hic # Define files and HiC read settings # 2. Hi-C Data Processing # Commented out IPython magic to ensure Python compatibility. # %%capture # !pip install hic-straw # import straw # Saves all self-contact data in Hi-C file # Adds chromosome interaction to genome dictionary, # interactions stored as list of contactRecord objects # Converts chromosome interactions from contactRecords into a list # Converts data in resolution units into particle numbers # Keeps entries where contact counts are higher than user-given threshold # and gets rid of atom self contacts #genome = self_chromosome_interactions(chromosomes, hicname, res) # Adds chromosome interaction to genome dictionary, # interactions stored as list of contactRecord objects # 3. LAMMPS Input & Data File Construction # function to plot 3D structure # function to check if next site is already occupied # function to create random 3D walk on lattice #if i % 100 == 0: #print(i, s-start, 's') #print('Stuck! Go back and find a new way... %s' % i) # function to create molecule tags # function to create bonds # function to create angles # function to create data file # number of chains # number of bonds # number of bond angles # coordinates of lattice points # molecule tags # indicates bonds between particles # indicates angles between particles # open datafile to write to # function to create input file # opens input file to write to # data file name # generates noise term for langevin # Defining LAMMPS properties # total number of particles # length of chains # lattice spacing # dimensions of box # Creating input and data files for simulation # 4. LAMMPS Simulation # Commented out IPython magic to ensure Python compatibility. # %%capture # # Creates LAMMPS executable # !wget https://github.com/lmhale99/atomman-demo/raw/main/lmp.gz # !gunzip lmp.gz # !chmod 755 lmp # # import subprocess # import shlex # Runs LAMMPS simulation based on input file # Reads dump file and extracts coordinate data #finds coordinates using "ITEM:" #creates a dictionary with timestep and coordinates # Displays simulation at timestep # 5. Analysis #!wget -cq https://www.dropbox.com/s/nkbuklgq2c9ewvw/nWTXa #Finds all contacts