fast forward to version 0.3.3

Author: elgw

CHANGELOG.md

@@ -0,0 +1,20 @@
+# CHANGELOG
+
+## 0.3.4
+- Cleaned up the makefile
+
+## 0.3.3
+- New behaviour, instead of `chromflock`, do `chromflock init`.
+- `cc2cpm` is baked into the binary `chromflock` and renamed to
+  `hic2cpm`
+- `sprite2cpm`, `any2string` and `string2any` also moved to `chromflock`.
+- documentation is based on markdown instead of troff.
+
+## 0.3.2
+ - Cleaned up the makefile.
+ - Added **--version** argument to **aflock**, **mflock**, and **cc2cpm**
+ - Fixed typos in the **cc2cpm** documentation.
+
+## 0.3.1
+ - Added script to make a deb package for Ubuntu. That is the
+   preferred way to install from now.

HAP1_100k.md

@@ -0,0 +1,59 @@
+# chromflock at 100kb on HAP1 using GPSeq
+
+This can be done on a system with 32 GB but will kill a system with 16 GB of ram (or use a lot of swapping).
+
+## Get Hi-C data
+```
+# 5,032,452,195 bytes
+wget https://data.4dnucleome.org/files-processed/4DNFI1E6NJQJ/@@download/4DNFI1E6NJQJ.hic
+```
+
+## Convert Hi-C matrix to raw matrices
+This step requires that straw is installed and can be found in the PATH.
+```
+# A 30877x30877 matrix
+~/code/chromflock_dev/util/python/chromflock_hic2H.py 4DNFI1E6NJQJ.hic 100000
+```
+
+## Translocate Hi-C data
+```
+~/code/chromflock_dev/util/python/chromflock_translocate.py  4DNFI1E6NJQJ.hic.H.double 4DNFI1E6NJQJ.hic.L0.uint8 100000
+```
+
+## Convert Hi-C (cc) data to contact probability matrix (cpm)
+ChrY is also removed at this stage.
+```
+mkdir 1000
+cd 1000
+cc2cpm --hfile ../4DNFI1E6NJQJ.hic.H_trans.double --lfile ../4DNFI1E6NJQJ.hic.H_trans.L.uint8 --nCont 8 --nStruct 1000
+```
+
+## Get GPSeq data (matlab)
+```
+L = read_raw('8.000000_1000_MAX.L.uint8', 'uint8');
+L = L(L<24); % Excluded Y if not already excluded
+
+G = loadGPSeq(L, ...
+    '/mnt/bicroserver2/projects/GPSeq/centrality_by_seq/SeqNoGroup/B170_transCorrected/all/B170_transCorrected.asB165.rescaled.bins.size100000.step100000.csm3.rmOutliers_chi2.rmAllOutliers.tsv', ...
+    100000);
+
+GN = log2(G);
+GN(GN>1) = 1;
+GN(GN<0) = 0;
+% GG: That's GPSeq score, ergo "centrality". I.e., greater at the center.
+GN = 1-GN;
+
+figure, plot(GN)
+write_raw(GN, 'G.double', 'double');
+```
+
+## Run chromflock
+```
+chromflock
+vim chromflock_gen
+# Add gpseq force
+vim chromflock.lua
+./chromflock_gen
+screen
+./chromflock_run
+```

INSTALL.md

@@ -0,0 +1,64 @@
+# INSTALLATION
+
+On Ubuntu 22.04
+
+``` shell
+cd src/lua-5.3.5/
+make linux # or pick another suitable target
+cd ../../
+make
+./makedeb-ubuntu_2204
+sudo apt install ./chromflock_*.deb
+```
+
+Check the dependencies section if there are any missing libraries. If
+you have another OS, please check the section **Custom install**.
+
+## Dependencies
+
+For compilation the following is needed:
+
+ * gcc
+ * libreadline -- to compile lua
+ * zlib -- to write `.gz` files
+ * pkg-config -- for the makefile
+ * libcairo -- optional, to write out the colour map
+ * libsdl2 -- optional, to compile mflock with live view (`mflock_sdl`)
+ * parallel
+ * liblua5.4.5 (included in src/lua-5.3.5/ )
+
+ On Ubuntu 22.04 these packages can be installed by
+
+  ``` shell
+  # Ubuntu
+  sudo apt-get install libcairo-dev
+  sudo apt-get install libreadline-dev
+  sudo apt install pkg-config
+  sudo apt-get install zlib1g-dev
+  sudo apt-get install libsdl2-dev
+  sudo apt-get install parallel
+  ```
+
+* lua, included but has to be built, on a Linux machine that corresponds to:
+  ```
+  cd src/lua-5.3.5/
+  # All possible platforms: aix bsd c89 freebsd generic linux macosx mingw posix solaris
+  # e.g., on linux use:
+  make linux
+  ```
+  Please check the Lua documentation if you are on another platform.
+
+
+## Custom install
+
+If your system does not support deb files you can still install
+chromflock. Typically by
+
+``` shell
+make
+sudo make install
+```
+
+please check what the **makefile** does before doing this. If you
+change the **DESTDIR** OR **PREFIX** variables please update the
+script called **chromflock** to reflect these changes.

NOTES.md

@@ -0,0 +1,8 @@
+## Files:
+```
+# Initial MATLAB prototype for the M-step
+src/mstep.m
+
+# Implementation of the 'error' function and it's gradient:
+src/functional.*
+```

README.md

@@ -1,44 +1,42 @@
 # CHROMFLOCK
 
-![example structure](cf_000038.png)
+![example structure](doc/cf_000038.png)
 
+## Introduction
 
-## Usage:
-Compile, if gcc, pkg-config and zlib are installed, compilation is as simple as:
+**chromflock** is a software package to deconvolve bulk Hi-C data into
+putative single-cell structures. The overall scheme is borrowed from
+[PGS](https://www.github.com/alberlab/PGS) although it is not a direct
+fork.
 
-```
-make all
-```
+Chromflock will:
+ - Work with haploid as well as diploid structures.
+ - Can integrate GPSeq data for radial preferences.
+ - Supports spherical as well as an ellipsoidal domain for the beads.
 
-To make the binaries and then man pages is visible, use:
+The documentation is not complete at the moment and it is suggested
+that anyone interested in chromflock start by reading [Kalhor et al,
+2012](https://doi.org/10.1038/nbt.2057), especially the supplementary
+materials since much of the terminology used here can be traced back
+to that paper.
 
-```
-here=$(pwd)
-export PATH=$PATH:$here/bin/
-export MANPATH=$MANPATH:$here/man/
-```
-possibly add that to your `.bashrc` or whatever is used on your system.
+Installation instructions can be found in [INSTALL.md](INSTALL.md) and
+some example usage in [USAGE.md](USAGE.md).
 
-To get started, read the man page for `chromflock`, i.e.:
+In addition to csv files, chromflock will use chimera marker files
+(.cmm) to write out structures for viewing. These can be opened by
+[UCFS chimera](https://www.cgl.ucsf.edu/chimera/) and also by
+[nua](https://www.github.com/elgw/nua/).
 
-```
-man chromflock
-```
+## References
 
-in general the sequence of command needed to start an experiment are:
-```
-chromflock
-vim chromflock_gen
-./chromflock_run
-```
-
-If `chromflock_run` was aborted for some reason, inspect `status.txt` to see the last thing that was finished and continue from any line, L, by
-```
-bash < (sed -n 'L,$p' chromflock_run)
-```
-
-# References
-For the RNG and Normal distribution:
- -  [A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers](http://www.tandfonline.com/doi/abs/10.1080/00949655.2015.1060234).
+chromflock was used in
+ - [Girelli et al. 2020](https://www.nature.com/articles/s41587-020-0519-y).
 
+The ideas behind it can be found in the following papers:
+ - [Kalhor et al, 2012](https://doi.org/10.1038/nbt.2057)
+ - [Tjong et al, 2016](http://dx.doi.org/10.1073/pnas.1512577113)
+ - [Hua et al, 2018](http://dx.doi.org/10.1038/nprot.2018.008) [github](https://www.github.com/alberlab/PGS)
 
+For random numbers chromflock uses
+ -  [McFarland, 2014](http://www.tandfonline.com/doi/abs/10.1080/00949655.2015.1060234) [code](https://github.com/cd-mcfarland/fast_prng).

TODO.md

@@ -0,0 +1,140 @@
+# Roadmap / to do
+
+## General todo:
+ - [ ] Fix so that `chromflock_gen` passes the gpseq file properly.
+
+## Molecular dynamics, M-step
+
+The M-step is where the model is updated. Initially the model is built from random positions, then at each following step:
+ 1. The contact indication matrix, W, for each structure is updated (A-step)
+ 2. New positions, X, are found in the structure that minimized the error functional (which includes W).
+
+Our implementation of the M-step uses Simulated Annealing with a Brownian force, much like what Peter Cook and Davide Marenduzzo does. Alber does not provide any details on what they are doing. However they refer to Russel2012, but not much details can be found their either. Since they published the Nat. Prot. paper the details are at least available in their source code.
+
+### TODO
+ - [ ] a lua script should be the only way to set the dynamics parameters.
+ - [ ] consider different stopping criteria, for example based on the radius of gyration.
+ - [ ] Fix so that mflock crashes gracefully when fRad is specified but no `-r` file is supplied.
+ - [ ] Append to log file if it already exists?
+ - [ ] Warn if largest `max(R) > 1-r0`. Or convert `R' = R*(1-r0)`? Same goes of the output profiles, should `r` be output or `r*(1-r0)`? (`r0` is the bead radius).
+ - [ ] Alber uses k=10 for consecutive beads (and k=1 for non consecutive).
+ - [ ] Fix so that the supplied random seed is used.
+ - [x] Put live view in main thread and optimization in child in order for SDL to be able to capture events on OSX.
+ - [x] Implement ellipsoidal geometry.
+ - [x] See if dInteraction = 4 is better than dInteraction = 3. -- seems slightly better. Changing default. See `fconf.dInteraction` in `md.c`.
+ - [x] Chimera script to add extra annotations and save automatically -- see `util/chromflock_images.py`
+
+ - [x] Removed any dependencies of GSL -- not necessary any more since only MD is used.
+ - [x] Switched RNG and Normal Distribution generator to [A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers](http://www.tandfonline.com/doi/abs/10.1080/00949655.2015.1060234).
+( Previously used Zignor, copied from Jurgen A. Doornik (2005). "An Improved Ziggurat Method to Generate Normal Random Samples" (PDF). Nuffield College, Oxford, together with [PCG](http://www.pcg-random.org) which should be better than a Marsenne Twister in terms of statistical properties and has about the same speed.)
+ - [x] Use 3D Gaussian for the Brownian force. (Used a force vector drawn uniformly from a 3D ball before).
+ - [x] Long option names, i.e. not just single letter names for the parameters.
+ - [x] Consider playing with the radius of the simulation (equivalent with the bead radius), one of the *tricks* used by Alber. -- Faster convergence when the volume quotient of the beads is smaller. However when doing that the beads does not use the full domain very well and the structures look very unrealistic.
+ - [ ] Consider to change from double to single precision. The instructions might not be faster for SP, but the memory locality will be better. -- Probably not worthwhile.
+ - [x] Any reason to use the varying `k_{con}` described on p. 64? -- Has possibly something todo with the *entropic forces* and convergence speed. -- Se below.
+ - [x] Use a decreasing volume exclusion force during the simulations. By hypothesis this reduces the bias towards *old* contacts, i.e., contacts introduced at an early stage. To be specific, the old contacts are the ones closest to the diagonal. This means that the diagonal is well represented among the contacts in the structures but non the off-diagnoal elements. This can be seen as a lack of of sharpness at the edges of the blocks.
+ - [x] Write error log to separate file, or use other mechanism to track crashes. -- Added `--joblog parallel.log` to `parallel`.
+ - [x] Is there a bias in the brownian force? (that brings the structure to the side) -- there was. Issue resolved and unit test in place.
+- [x] Implement Verlet integration. No significant changes in performance.
+- [x] How many "brownian" steps and what falloff?
+ - [x] Delta t, switched from 1 to 5 2019-04-17. What would be a good value?
+ - [x] Implement live view of current state. -- available from SDL with the -a swithch (when compiled with -DSDL).
+ - [x] Add a COM force (per chr).
+ - [x] Switch to radial data in the "structure summary"
+ - [x] Timing does not work in the -D mode.
+ - [x] Use a label vector to generate cmms.
+ - [x] Include L - label, both for cmm creation and for various per chromosome things.
+ - [x] Integrate with GPSeq data.
+ - [x] Add comparison of all `err` and all `grad` in the unit tests.
+ - [x] Make bead size a parameter or set so that the volume occupancy is 20% (or some other number).
+ - [x] Output coordinates for chimera.
+ - [x] HASH table to find points in contact.
+ - [x] Use a list of pairwise interactions rather than the full matrix.
+
+And for the MD implemenation
+ - [ ] Keep track of Ek, and potentially set it to a specific value as in IMP.
+ - [ ] For the above point, use 1-diagonal only.
+ - [x] Implement a force that brings chrs together to their centre of mass as Alber does.
+
+### Parallelization
+
+ - [x] Structure optimizations can be carried out in parallel, one thread per structure. The jobs are distributed with GNU parallel.
+
+### Memory
+
+## A-step
+ - [ ] Consider re-running the m-step until the number of contacts falls below some threshold.
+ - [ ] Consider writing coords as raw double data.
+ - [ ] Add counting of failed contacts to the `-F` mode.
+ - [ ] Add a flag to ignore inter contacts, then there won't be any need to generate a separate A-matrix when trying without inter contacts.
+ - [x] Allow diploid experiments.
+ - [x] Use zlib to compress the W matrices.
+ - [x] Consider *even* assignment in the sense that each structure get the same number of contact initially.
+ - [x] Implement re-assignment, and improve the assignment.
+ - [x] Possible error when setting the bead radius from volume quotient?!
+ - [x] Use zlib to write the cmms.
+ - [x] Enable chromflock to use X from previous simulation with updated A.
+ - [x] Implement the HiC-normalization from p. 53. The transform was actually already performed on the data!
+ - [x] Implement the sphere contact probability matrix from p. 65.
+ - [x] Determine W, see p. 66. Rather assign the contacts as in Tjong2016, i.e., the structures that already have the shortest distances will have to do the job.
+And for flock:
+ - [x] Read sphere contact probability matrix A.
+ - [x] Initialize new experiment here or separate? -- in bash script
+ - [x] Read and write the contact indicator tensor W.
+ - [x] Write radial profile at -F
+
+### Memory
+ - Coordinates: Surprisingly not that much memory is required. Let's get an estimate:
+ 10,000 structure x 5,000 points x 3 dimensions = 1,200 Mb (64 bit floats). Currently aflock is using single precision to represent the coordinates. Using 10,000 structures x 3,300 points requires about .5 Gb.
+
+ - Constraints: Each constraint requires 8 bytes `(uint_32, uint_32)`. By keeping separate lists for separate types of constraints there is no need to story the type.
+ 10,000 structures x 5,000 constraints = 400 MB.
+
+### Parallelization
+ - [x] The assignment process could be parallelized over each contact to be handed out, i.e., the qsorts. -- The generation of the contact distances is now parallelized using `pthread`. Initially all available cps are used.
+
+## General and Random notes
+
+ - [ ] Can we estimate the entropy based on the gradient norm?
+
+ - [ ] Does the simulation tool yield end-to-end distances of L^{3/5} for unperturbated chains?
+
+ - [ ] Update makefile and folder layout, have a look at [stackoverflow](https://stackoverflow.com/questions/7004702/how-can-i-create-a-makefile-for-c-projects-with-src-obj-and-bin-subdirectories)
+
+ - [x] Make Hi-C look like TCC!
+
+ - [ ] Estimate number of failed contacts after each `mflock`-run in order to decide if another round is needed of `aflock` before going to smaller thetas.
+
+ - [ ] ~~libexpat for settings?~~
+
+ - [x] lua for interactions/setting md/sa schedules?
+
+ - [x] The chr4 trick?! -- No
+
+ - [x] Rewrite `chromflock_run` so that it is easier to continue/restart failed jobs. New pipeline: `chromflock` generates `chromflock_gen` which should be edited and run. That script does in turn create `chromflock_run` which is a linear batch of jobs to be run. If anything in `chromflock_run` fails it should be easier to restart from that position.
+
+ - [ ] Are beads of different type handled differently by Alber?
+
+ - [ ] What are the results if only the structures with average number of contact restraints are studied?
+
+ - [x] How many inter vs intra contacts are there when theta = 1, .2, ... , 0.01. (Is the algorithm agnostic to anything but the intra contacts? That would explain random orientation of the chromosomes).
+
+ - [x] What if delta t is increased? -- I think that I got that parameter quite well.
+
+ - PGP is only for diploid cells.
+
+ - Maybe it makes sense to 'kill your darlings', i.e., open up for the possibility to not use all of the initial structures (discard or replace). We don't expect the nucleus to be completely random so random initializations might go wrong.
+
+ - chromflock is snappier than chromoflock. chrof, crof.
+
+ - Convergence critera shouldn't have to be that strict until final theta.
+
+  - A very limited number of theta values, i.e., only few assignment steps (A) are used. This could potentially lead to the inclusion of conflicting constraints. However, it is extremely costly to increase this number...
+
+ - The Hi-C maps does probably implicitly suggest radial positions. However what radial positions that are assigned to each loci depends on the reconstruction method. No reconstruction method is perfect, rather the opposite. This means that GPSeq adds extra constraints, i.e., reduces the space of probable states/configurations.
+
+ - There is an upper limit on the number of neighbours each bead can have. Geometry tells us 12, see the wikipedia page on sphere packings.
+
+## Further devolopment
+ - [ ] Use local temperature like [Agrawal](https://www.biorxiv.org/node/96892.full)
+ - [ ] Use proteins like [Cook and Marenduzzo](http://dx.doi.org/10.1093/nar/gkw135)