Skip to content

Instantly share code, notes, and snippets.

@floybix floybix/core.clj

Created Feb 14, 2017
Embed
What would you like to do?
GRN (gene regulatory network) experimental discrete formulation
(ns org.nfrac.gron.core
(:require [clojure.spec :as s]
[clojure.spec.gen :as gen]))
(s/def ::sep #{-1})
(s/def ::base (-> nat-int?
(s/with-gen #(s/gen (s/int-in 0 10)))))
(s/def ::seq-gene
(s/cat :sep1 ::sep
::promoter (s/* ::base)
:sep2 ::sep
::repressor (s/* ::base)
:sep3 ::sep
::produce (s/* ::base)
:sep4 ::sep
::consume (s/* ::base)))
(s/def ::genome
(s/+ ::seq-gene))
(s/def ::grn
(s/keys :req [::genome]))
(comment
(s/exercise ::seq-gene))
(s/def ::base-counts
(s/map-of ::base (-> nat-int? (s/with-gen #(s/gen (s/int-in 1 4))))))
(s/def ::cell ::base-counts)
(s/def ::promoter ::base-counts)
(s/def ::repressor ::base-counts)
(s/def ::produce ::base-counts)
(s/def ::consume ::base-counts)
(s/def ::gene (s/keys :req [::promoter
::repressor
::produce
::consume]))
(comment
(s/exercise ::gene))
(defn regulator-matched?
[regulator cell]
(and (seq regulator)
(every? (fn [[base n]] (>= (cell base 0) n)) regulator)))
(s/fdef regulator-matched?
:args (s/cat :regulator ::base-counts
:cell ::cell))
(defn gene-expressed?
[gene cell]
(and (regulator-matched? (::promoter gene) cell)
(not (regulator-matched? (::repressor gene) cell))))
(s/fdef gene-expressed?
:args (s/cat :gene ::gene
:cell ::cell))
(defn expressed-genes
[genes cell]
(loop [genes genes
bound-cell cell
expressed ()]
(if-let [gene (first genes)]
(if (gene-expressed? gene bound-cell)
(recur (rest genes)
(merge-with - bound-cell (::promoter gene) (::repressor gene))
(conj expressed gene))
;; gene not expressed
(recur (rest genes) bound-cell expressed))
;; checked all genes
expressed)))
(s/fdef expressed-genes
:args (s/cat :genes (s/every ::gene)
:cell ::cell)
:ret (s/every ::gene))
(defn gene-effects
[genes]
(let [pro (apply merge-with + (map ::produce genes))
con (apply merge-with + (map ::consume genes))]
(merge-with + pro (zipmap (keys con) (map - (vals con))))))
(s/fdef gene-effects
:args (s/cat :genes (s/every ::gene))
:ret (s/every-kv ::base int?))
(defn step
[genes cell inputs]
(let [icell (merge-with + cell inputs)
xg (expressed-genes genes icell)
ge (gene-effects xg)]
(reduce (fn [cell [base delta]]
(assoc cell base
(max 0 (+ (cell base 0) delta))))
cell ;icell
ge)))
(s/fdef step
:args (s/cat :genes (s/every ::gene)
:cell ::cell
:inputs ::base-counts)
:ret ::cell)
(defn genes-form
[seq-genome]
(->> seq-genome
(s/conform ::genome)
(map (fn [m]
(reduce (fn [m k]
(assoc m k (frequencies (get m k))))
m
[::promoter ::repressor ::produce ::consume])))))
(s/fdef genes-form
:args (s/cat :seq-genome ::genome)
:ret (s/every ::gene))
(comment
(require 'proto-repl-charts.charts)
(let [genes (genes-form (gen/generate (s/gen ::genome)))
cell-ts (reductions (fn [cell in]
(step genes cell in))
(gen/generate (s/gen ::base-counts))
(repeat 30 {}))
bases (range 10)
base-ts-map (reduce (fn [m b]
(assoc m b (map #(get % b 0) cell-ts)))
{}
bases)]
(println (count genes) "genes")
(proto-repl-charts.charts/line-chart
"Random GRNs"
base-ts-map
{:labels (range (count cell-ts))})))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.