ipynb

tutorials/optimization/chr_nuclear_landmarks/calculate_contact_prob.py → openNucleome/utility/DamID_TSASeq_calculation.py

@@ -9,27 +9,9 @@
 import sys
 import os
 
-iter_num        = int(sys.argv[1])
-sigma_tanh          =  4
-rc_tanh             =  0.75
-
-N_chr_beads         = 60642
-
-first_frame                 = 500
-N_nucleolus_particles       = 300
-N_speckles_particles        = 1600
-N_lamina_particles          = 8000
-radius_nucleus              = 13.0
-
-"""Info files"""
-gLength             = np.loadtxt("mol_info/gLengthFile.txt",dtype=int)
-maternalIdx         = np.loadtxt("mol_info/maternalIdxFile.txt",dtype=int)
-paternalIdx         = np.loadtxt("mol_info/paternalIdxFile.txt",dtype=int)
-""""""
-
 def process_traj_win(win_num):
 
-    traj_data                   = mda.coordinates.LAMMPS.DCDReader("../../../tutorials/HFF_100KB/DUMP_FILE.dcd")
+    traj_data                   = mda.coordinates.LAMMPS.DCDReader("../HFF_100KB/SphericalNucleus/DUMP_FILE.dcd")
     N_frame                     = len(traj_data)-first_frame
     exp_tsa_seq                 = np.zeros(N_chr_beads)
     exp_damid                   = np.zeros(N_chr_beads)
@@ -102,5 +84,3 @@ def process_traj_win(win_num):
     np.save("contact_prob/cvInd.txt_%d"%win_num,np.sum(damid_all_frames_haploid,axis=0))
     np.save("contact_prob/cvInd_tsa.txt_%d"%win_num,np.sum(tsaseq_all_frames_haploid,axis=0))
 
-win_to_run = int(sys.argv[2])
-process_traj_win(win_to_run)

tutorials/analysis_code/.gitignore

@@ -0,0 +1,2 @@
+chr_chr/
+chr_nuclear_landmarks/

tutorials/analysis_code/Compute_contact_prob.ipynb

@@ -0,0 +1,119 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "3ef4e4c3-bbff-4170-9390-b51feae5efb6",
+   "metadata": {},
+   "source": [
+    "## Import the necessary packages"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be874cb0-b8ff-4af7-ac51-7fd3c914291a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "sys.path.insert(1, '../../openNucleome/utility')\n",
+    "import numpy as np\n",
+    "import MDAnalysis as mda\n",
+    "import sklearn.cluster as sk\n",
+    "import matplotlib.pyplot as plt\n",
+    "import copy\n",
+    "import subprocess\n",
+    "from joblib import Parallel, delayed\n",
+    "import scipy.stats\n",
+    "import os\n",
+    "from DamID_TSASeq_calculation import process_traj_win"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b3ae3b2-9dc7-4d74-a860-0b35631f3c28",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_tanh          =  4\n",
+    "rc_tanh             =  0.75\n",
+    "\n",
+    "N_chr_beads         = 60642\n",
+    "\n",
+    "first_frame                 = 500\n",
+    "N_nucleolus_particles       = 300\n",
+    "N_speckles_particles        = 1600\n",
+    "N_lamina_particles          = 8000\n",
+    "radius_nucleus              = 13.0\n",
+    "\n",
+    "\"\"\"Info files\"\"\"\n",
+    "gLength             = np.loadtxt(\"mol_info/gLengthFile.txt\",dtype=int)\n",
+    "maternalIdx         = np.loadtxt(\"mol_info/maternalIdxFile.txt\",dtype=int)\n",
+    "paternalIdx         = np.loadtxt(\"mol_info/paternalIdxFile.txt\",dtype=int)\n",
+    "\"\"\"\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d5c1c187-a093-40c5-b67f-94ffc8c8806e",
+   "metadata": {},
+   "source": [
+    "## Compute the DamID and TSASeq"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29222bbb-8773-4571-ac97-5756c1570785",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "iter_num = int(sys.argv[1])\n",
+    "win_to_run = int(sys.argv[2])\n",
+    "\n",
+    "process_traj_win(win_to_run)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15f5efce-885f-4716-80fa-7e66a8333457",
+   "metadata": {},
+   "source": [
+    "## Compute the chromosome contact probabilities"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7dd7e0df-e4aa-491f-842a-65f44879c955",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "subprocess.call([\"gfortran -o chromsome_contact_calculation ../../openNucleome/utility/chromsome_contact_calculation.f90\"],shell=True,stdout=subprocess.PIPE)\n",
+    "subprocess.call([\"./chromsome_contact_calculation ../HFF_100KB/SphericalNucleus/DUMP_FILE.dcd 501 -1 ./contact_prob/ 1 counter.txt\"],shell=True,stdout=subprocess.PIPE)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}