rlm: src/rlm/dna_melting.clj comparison

comparison src/rlm/dna_melting.clj @ 0:78a630e650d2

initial import

author	Robert McIntyre <rlm@mit.edu>
date	Tue, 18 Oct 2011 00:57:08 -0700
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:78a630e650d2
+(ns rlm.dna-melting
+(:refer-clojure :only [])
+(:require rlm.ns-rlm mobius.base))
+(rlm.ns-rlm/ns-clone mobius.base)
+(use :reload-all 'rlm.web-stuff)
+(def pre "/home/r/MIT/6/6.122/DNA/samples/")
+(defrecord dna-data [temp glow name temp-units glow-units]
+clojure.lang.IFn
+(invoke [this k] (get this k))
+(invoke [this k not-found] (get this k not-found))
+(toString [this] "whatever"))
+(defmethod print-method dna-data [d w]
+	   (.write w (.toString d)))
+(defmethod print-dup dna-data [d w]
+	   (print-method d w))
+(defn load-data [source]
+(let [data (read-dataset source :delim \tab)
+	temp (vec ($ :col0 data))
+	glow (vec ($ :col1 data))]
+(dna-data. temp glow
+	        (last (re-split #"/" source))
+		"Raw Temperature (Volts)"
+		"Raw Glow (Volts)")))
+(def ?-1 (load-data (str pre "sample4-1.txt")))
+;; I accidentally removed the cuvette while it was still recording!
+;; drop last 5 seconds of invalid data!
+(def ?-2* (load-data (str pre "sample4-2.txt")))
+(def ?-2 (assoc ?-2*
+	   :temp (drop-last 500 (:temp ?-2*))
+	   :glow (drop-last 500 (:glow ?-2*))))
+(def ?-3 (load-data (str pre "sample4-3.txt")))
+(def A-1 (load-data (str pre "sampleA-1.txt")))
+(def A-2 (load-data (str pre "sampleA-2.txt")))
+(def A-3 (load-data (str pre "sampleA-3.txt")))
+(def B-2 (load-data (str pre "sampleB-2.txt")))
+(def B-3 (load-data (str pre "sampleB-3.txt")))
+(def C-1 (load-data (str pre "sampleC-1.txt")))
+(def C-2 (load-data (str pre "sampleC-2.txt")))
+(def C-3 (load-data (str pre "sampleC-3.txt")))
+(def pure-bleaching (load-data (str pre "sampleB-1(failed-to-heat).txt")))
+(def double-cycle (load-data (str pre "sampleC-3(double).txt")))
+(def triple-cycle (load-data (str pre "sampleC-3(triple).txt")))
+(def all-samples
+[?-1 ?-3 A-1 A-2 A-3 B-2 B-3 C-1 C-2 C-3])
+(def everything
+(concat [pure-bleaching double-cycle triple-cycle] all-samples))
+(defn dna-plot [d]
+(view
+(->
+(visual (:temp d))
+(set-y-label (:temp-units d))
+(set-x-label "Time (deci-seconds)")
+(set-title  (:name d))))
+(view
+(->
+(visual (:glow d))
+(set-y-label (:glow-units d))
+(set-x-label "Time (deci-seconds)")
+(set-title  (:name d))))
+d)
+(defn dna-melting [d]
+(->
+(scatter-plot (:temp d) (:glow d))
+(set-y-label (:glow-units d))
+(set-x-label (:temp-units d))
+(set-title (:name d))))
+(defmethod visual clojure.lang.LazySeq [ls]
+	   (visual (vec ls)))
+(defmethod visual clojure.lang.PersistentVector [v]
+	   (set-y-range (xy-plot (vec (range (count v)))  v)
+			      (apply min v) (apply max v)))
+(defmethod visual dna-data [d]
+	   (dna-plot d))
+;;; Filter V_RTD and V_F to reduce noise
+;(let [a (int 1.8e4) b (int 2.2e4)]
+;  (visual (map (vec (temperature triple-cycle)) (range a b)))
+;  (visual (map (vec (dna-fraction triple-cycle)) (range a b))))
+;; Differentiate the resulting function with respect to temperature.
+(def test-temp-range (vec (range 260 380 0.1)))
+(defn trans-print [& args] (apply println args) args)
+;; (def search
+;;      #(for [x (range -1000 1000 100) y (range -1000 1000 100)]
+;; 	(do
+;; 	  (print [x y])
+;; 	  (try (trans-print (non-linear-model
+;; 			     (ideal-dna-fraction 30e-6)
+;; 			      test-dna-fraction test-temperature [x y]
+;; 			     :method :newton-raphson))
+;; 	       (catch Exception e (trans-print nil))))))
+;;(non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-184  -71e3])
+;;(time (def answer (non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-200  -90e3])))
+;;we're assuming that the SYBR Green 1 doesn't glow at all when there is only
+;;single stranded DNA, but it actually does. V_f only ever gets down to ~0.1 V
+;;even with everything melted.
+;;correct for thermal capacatance
+(defn find-max [v]
+(first (reduce (fn [a b] (if (> (last b) (last a)) b a)) (indexed v))))
+(defn find-min [v]
+(first (reduce (fn [a b] (if (< (last b) (last a)) b a)) (indexed v))))
+(defn plot-thermal-delay []
+(let [a (int 2.5e3) b (int 5e3)]
+(dna-plot (assoc A-1
+		 :temp (map (vec (:temp A-1)) (range a b))
+		 :glow (map (vec (:glow A-1)) (range a b))))))
+(defvar thermal-delay
+(let [a (int 2.5e3) b (int 5e3)]
+(-
+	(find-min (map (vec (:glow A-1)) (range a b)))
+	(find-max (map (vec (:temp A-1)) (range a b)))))
+"There are two sources of thermal delay: thermal capacatance
+of the cuvette holder and thermal capacatance from the time
+it takes DNA to unwind.
+Use A-1 to find the thermal delay. There was almost no delay
+between heating and cooling for A-1, and thus the extremum of
+temperature and glow should temporally coincide if there were
+no delays. It turns out that there is about a 30 second lag!" )
+(defn thermal-delay-correction*
+"There is a certain delay between the RTD reporting changes in
+temperature and the temperature of the DNA actually changing.
+This shifts the temperature values so that the delay is eliminated."
+[thermal-delay #^dna-data d]
+(assoc d
+:name (str (:name d) " |thermal-correction| ")
+:temp (drop-last thermal-delay (:temp d))
+:glow (drop thermal-delay (:glow d))))
+(def thermal-delay-correction (partial thermal-delay-correction* thermal-delay))
+;;correct for photo-bleaching
+(defn general-linear [[A B] x] (+ (* A x) B))
+(defvar bleaching-coef nil
+;; (let [[A B]
+;; 	 (:coefs
+;; 	  (let [a (int 3e4) b (length (:glow pure-bleaching))]
+;; 	    (non-linear-model
+;; 	     general-linear
+;; 	     (map (vec (:glow pure-bleaching)) (range a b))
+;; 	     (range a b)
+;; 	     [1 1])))]
+;;   A)
+"This function corrects V_f for photo-bleaching over a period of time.
+Photo-bleaching was found to be linear(!) over quite a long period of time.
+Data is taken from sample B-1, which I left in the holder for an hour
+with no heating")
+(defn bleaching-correction*
+[bleaching-coef #^dna-data d]
+;;the correction is clibrated to apply to the raw data only
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+;;we need to modulate this correction factor
+(let [baseline (first (drop 3e4 (:temp pure-bleaching)))]
+(assoc d
+:name (str (:name d) " |bleaching-correction| ")
+:glow (map (fn [[i n]] (- n (*
+				   (/ n baseline)
+				   ((partial general-linear [bleaching-coef 0]) i))))
+		 (indexed (:glow d))))))
+(def bleaching-correction (partial bleaching-correction* bleaching-coef))
+;; compensate for residual fluorescence when everything is melted
+(defn residual-fluorescence
+"the SYBR1 green dye doesn't completly stop glowing once everything
+is melted."
+[#^dna-data d]
+(let [low (apply min (:glow d))
+	high (apply max (:glow d))]
+(assoc d
+:name (str (:name d) " |residual-fluorescence| ")
+:glow (map (fn [f]
+		   (- f (* (- 1 (/ f high)) low))) (:glow d)))))
+;; extract the cooling phase
+(defn cooling-phase
+"first stretch where the temperature is going down for a while
+equivalent to where the f curve starts to recover"
+[#^dna-data d]
+(let [t (:temp d)
+	f (:glow d)
+	dividing-line (reduce (fn [a b] (if (< (a 1) (b 1)) a b)) (indexed f))
+	t-new (drop (dividing-line 0) t)
+	f-new (drop (dividing-line 0) f)]
+(assoc d
+:name (str (:name d) " |extract cooling| ")
+:temp (vec t-new) :glow (vec f-new))))
+;; filter noise
+;;(http://clojure-log.n01se.net/date/2009-10-13.html by ambient
+;;this is such a cool transform!!!
+(defn lo-pass
+[d coll]
+(reductions #(+ (* %2 d) (* %1 (- 1.0 d))) coll))
+;;)
+(def reasonable-low-pass (partial lo-pass 1e-2))
+(defn low-pass-correction
+"filter the raw data through a low pass filter to remove
+electronic noise"
+[#^dna-data d]
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+(assert (= (:temp-units d) "Raw Temperature (Volts)"))
+(assoc d
+:name (str (:name d) " |low-pass| ")
+:temp (reasonable-low-pass (:temp d))
+:glow (reasonable-low-pass (:glow d))))
+;; correct for units
+(defn RTD-voltage->temperature
+"transform RTD voltage to temperature.
+R_RTD = 1000 + 3.85 * (T - 273)
+T = ((R_RTD - 1000) / 3.85 ) + 273
+V_RTD = 15 * R_RTD / (R + R_RTD)
+R_RTD = R / (15/V_RTD - 1)"
+[coll]
+(let [R 14.97e3
+	V 15]
+(letfn [(T [R_RTD] (+ 273 (/ (- R_RTD 1000)  3.85))) ;; from lab notrd
+	    (R_RTD [V_RTD] (/ R (- (/ V V_RTD) 1)))]
+(map (comp T R_RTD) coll))))
+(defn normalize
+"Transform photodiode current to relative fluorescence."
+[coll]
+(let [low (apply min coll)
+	high (apply max coll)]
+(map #(/ (- % low) (- high low)) coll)))
+(defn unit-transform
+"turn raw voltage into temperature and dna-concentration"
+[#^dna-data d]
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+(assert (= (:temp-units d) "Raw Temperature (Volts)"))
+(assoc d
+:name (str (:name d) " |units| ")
+:temp-units "Temperature (Kelvin)"
+:glow-units "DNA-Fracton"
+:temp (RTD-voltage->temperature (:temp d))
+:glow (normalize (:glow d))))
+(defn temp-unit-transform
+"turn raw voltage into temperature and dna-concentration"
+[#^dna-data d]
+(assert (= (:temp-units d) "Raw Temperature (Volts)"))
+(assoc d
+:name (str (:name d) " |temp-units| ")
+:temp-units "Temperature (Kelvin)"
+:temp (RTD-voltage->temperature (:temp d))))
+(defn glow-unit-transform
+"turn raw voltage into temperature and dna-concentration"
+[#^dna-data d]
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+(assoc d
+:name (str (:name d) " |glow-units| ")
+:glow-units "DNA-Fracton"
+:glow (normalize (:glow d))))
+(defn thermal-efficiency
+[#^dna-data d]
+(assert (= (:temp-units d) "Temperature (Kelvin)"))
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+(assoc d
+:name (str (:name d) " |thermal-efficiency| ")
+:glow nil))
+(defn remove-first-values
+"remove the first minute or so of data so that the initial temperature of
+the cuvette itself doesn't matter"
+[#^dna-data d]
+(assoc d
+:name (str (:name d) " |remove-first| ")
+:temp (drop 600 (:temp d))
+:glow (drop 600 (:glow d))))
+(defn diff [coll]  (map - (rest coll) coll ))
+(defn temp-transform
+"my assumption will be that at the very end of the temperature data, the
+cuvette is exactly the same temperature as the RTD.  Then, I use newton's
+heat transfer and numerical integration to derive the heat of the cuvette
+for previous times"
+[k-value temperature-data]
+(let [T1 (vec (reverse temperature-data))
+;; want fast indexed access, and instead of doing the problem in reverse
+;; just reverse the vector at the start, treat it normally, and
+;; then reverse at the end
+	initial-cuvette-temp (first T1)
+	initial-RTD-temp     (first T1)
+	step (fn [prev-T2 current-T1]
+	  (let [new-T2 (+ prev-T2 (* k-value (- current-T1 prev-T2)))]
+	    new-T2))]
+(reverse (reductions step initial-cuvette-temp T1))))
+(def k-value 0.1)
+(defn temp-corr-2 [k-value T1*]
+(let [T1 (reverse T1*)]
+(reverse  (cons (first T1)  (map + (map (partial * k-value) (diff T1)) (rest T1))))))
+(defn temp-corr-3 [initial-cuvette-temp k-value T1]
+(let [step
+	(fn [prev-T2 current-T1]
+	  (+ prev-T2 (* k-value (- current-T1 prev-T2))))]
+(rest (reductions step initial-cuvette-temp T1))))
+(defn newton-temp-correction* [k-value #^dna-data d]
+(assoc d
+:name (str (:name d) " |newton-temp| ")
+:temp (temp-transform k-value (:temp d))))
+(defn newton-temp-correction2* [k-value #^dna-data d]
+(assoc d
+:name (str (:name d) " |newton-temp-reverse| ")
+:temp (temp-corr-2 k-value (:temp d))))
+(defn newton-temp-correction3* [init-temp k-value #^dna-data d]
+(assoc d
+:name (str (:name d) " |newton-temp-reverse| ")
+:temp (temp-corr-3 init-temp k-value (:temp d))))
+(defn turn-e [n] (format "%5.2e" n))
+(defn newton-temp-correction4*
+"now with AIR!!!"
+[k-block-sample k-sample-air initial-sample-temp room-temp #^dna-data d]
+(let [old-temp (:temp d)
+	step (fn [prev-T2 current-T1]
+	       (+ prev-T2
+		  (* k-block-sample (- current-T1 prev-T2))
+		  (* k-sample-air   (- room-temp  prev-T2))))
+	new-temp (rest (reductions step initial-sample-temp old-temp))]
+(assoc d
+:name (str (:name d)
+		 " |temp-adjusted {"
+		 "initial-sample:" (turn initial-sample-temp) " "
+		 "room-temp:" (turn room-temp) " "
+		 "k-block-sample:" (turn-e k-block-sample) " "
+		 "k-sample-air:" (turn-e k-sample-air) "}")
+:temp new-temp)))
+(defn plot-adjusted-temp
+"visualize the temperature correction"
+[adjuster #^dna-data d]
+(view
+(add-lines
+(visual (:temp d))
+(vec (range (count (:temp d))))
+(vec (:temp (adjuster d))))))
+;; process data
+(defn process [#^dna-data d]
+(cooling-phase
+(unit-transform
+(residual-fluorescence
+(bleaching-correction
+(low-pass-correction
+;(thermal-delay-correction
+	d))))))
+(defn first-pass
+"after this, the only thing stopping the melting curve from
+being a straight line is thermal capacatance"
+[#^dna-data data]
+(low-pass-correction
+;(remove-first-values
+(bleaching-correction
+data)))
+(defn explore [k-value]
+(dna-melting (newton-temp-correction* k-value (first-pass A-2))))
+(defn explore2 [k-value]
+(dna-melting (newton-temp-correction2* k-value (first-pass A-2))))
+(defn explore3 [init-temp k-value]
+(dna-melting (newton-temp-correction3* init-temp  k-value (unit-transform (first-pass A-2)))))
+;; get the three characteristic values we want.
+(defn ideal-dna-fraction* [concentration [delta-S delta-H] temperature]
+(let [R 8.31447215
+	  t temperature
+	  Ct concentration
+	  K_eq (exp (- (/ delta-S R) (/ delta-H (* R t))))
+	  f (/
+	     (+ 1 (* Ct K_eq)
+		(- (sqrt (+ 1 (* 2 Ct K_eq)))))
+	     (* Ct K_eq))]
+f))
+(def ideal-dna-fraction (partial ideal-dna-fraction* 30e-6))
+(defvar starting-values
+;;  [-184 -71e3]
+[-904 -32e4]
+"chosen from theoritical estimation")
+(defn analyze* [iterations #^dna-data d]
+(let [best-fit-results
+	(non-linear-model
+	 ideal-dna-fraction
+	 (:glow d)
+	 (:temp d)
+	 starting-values
+	:max-iter iterations)]
+best-fit-results))
+(def analyze (partial analyze* 200))
+(defn plot-results [d]
+(let [ans (analyze d)]
+(view (doto (xy-plot (:x ans) (:y ans))
+(set-title (:name d))
+(add-lines (:x ans) (:fitted ans))))
+ans))
+(defn zzz [d [A B]]
+(let [gar (vec (map (partial ideal-dna-fraction [A B]) (:temp d)))]
+(view (doto (xy-plot (:temp d) (:glow d))
+(set-title (:name d))
+(add-lines (:temp d) gar)))))
+(defn index-filter [pred coll]
+(when pred
+(for [[idx elt] (indexed coll) :when (pred elt)] idx)))
+;; we're going to estimate two correction factors,
+;; a constant capacative temperature delay,
+;; and a linear decrease in brightness.  We'll do this
+;; before unit-transforms.
+(defn zipmap-collect
+"Returns a map with the keys mapped to the corresponding vals.
+adds extra values to a list"
+[keys vals]
+(loop [map {}
+ks (seq keys)
+vs (seq vals)]
+(if (and ks vs)
+(recur (if (contains? map (first ks))
+		 (assoc map (first ks) (cons (first vs) (map (first ks))))
+		 (assoc map (first ks) (list (first vs))))
+(next ks)
+(next vs))
+map)))
+(defn little-process [thermal-cap bleaching #^dna-data d]
+((partial thermal-delay-correction thermal-cap)
+((partial bleaching-correction bleaching) d)))
+(def little-process (memoize little-process))
+;(defn lookup-map [#^dna-data d]
+(defn converge-function
+"this is the vertical discrepancy after applying the two transforms"
+[#^dna-data d [thermal-cap bleaching] x]
+(let [new-d (little-process thermal-cap bleaching d)
+	spreads
+	(map (vec (:glow new-d))
+	     (index-filter (partial = x) (vec (:temp new-d))))
+	spread (- (apply max spreads) (apply min spreads))]
+spread))
+(defn  blah  [k] (comp
+	   ( partial newton-temp-correction3* 280 k ) temp-unit-transform first-pass))
+;(def ggg (comp
+;	  dna-melting ( partial thermal-delay-correction* 120 ) (blah 0.0018)))
+;; Ensure that the melting function is uniquely valued by
+;;   combining samples with identical temperature values.
+(defn cooling-function
+"just put it in a map for this!"
+[source]
+(let [data (cooling-phase source)
+	lookup (apply sorted-map
+		    (interleave
+		     (:temperature data) (:dna-fraction data)))]
+{:temperature (vec (keys lookup))
+:dna-fraction (vec (map lookup (keys lookup)))}))
+(defn turn [n]
+(format "%5.2f" n))
+(defn bleach-correction-2 [bleaching-constant #^dna-data d]
+(assert (= (:glow-units d) "Raw Glow (Volts)"))
+(let [old-glow (vec (:glow d))
+	integral-old-glow (vec (reductions + old-glow))
+	correction (fn [t] (+
+			    (old-glow t)
+			    (* bleaching-constant (integral-old-glow t))))
+	new-glow (map + (map correction (range 0 (count old-glow))) old-glow)]
+;;glow is related to the number of particles,
+;;and the number destroyed is some fraction of
+;;the number present and glowing.
+(assoc d
+:name (str (:name d) " |bleaching-corrected constant = " (turn bleaching-constant) "| ")
+:glow new-glow)))
+(defn ultimate-correction
+([k-block-sample k-sample-air initial-sample-temp room-temp bleach-constant
+#^dna-data d
+]
+((comp
+(partial newton-temp-correction4*
+	       k-block-sample k-sample-air initial-sample-temp room-temp)
+temp-unit-transform
+(partial bleach-correction-2 bleach-constant)
+low-pass-correction) d))
+([{:keys [test-k-sample-air
+	    test-k-block-sample
+	    test-sample-temp
+	    test-room-temp
+	    test-bleach]}
+		#^dna-data d]
+(ultimate-correction
+(mean test-k-block-sample) (mean test-k-sample-air)
+(mean test-sample-temp) (mean test-room-temp) (mean test-bleach) d)))
+;;((uui 0.0020 0.0001 280 293) (bleach-correction-2 1.1e-5 ?-3))
+;;; this is for ?-3 after 35 iterations
+(def opt {:test-k-sample-air '(9.375630358472117E-4 9.375630359053611E-4), :test-k-block-sample '(0.0010343754438043107 0.0010343754455447197), :test-room-temp '(294.9804782861611 294.98047828674316), :test-sample-temp '(293.36120605294127 293.3612060546875), :test-bleach '(1.771081481128931E-5 1.7710814863094127E-5), :align 1.5796052678842228E-4})
+;;; C-1 after 20 iterations
+(def opt2
+{:test-k-sample-air '(9.936472587585448E-4 9.936491641998292E-4), :test-k-block-sample '(2.9390106201171873E-4 2.939580917358398E-4), :test-room-temp '(294.76531982421875 294.7653388977051), :test-sample-temp '(294.765625 294.765682220459), :test-bleach '(2.39045524597168E-5 2.3906250000000002E-5), :align 1.7562642331021098E-4})
+;;(time (def result-?-3 (doall (take 10 (iterate (partial improve ?-3) start-values-overestimates)))))
+;;(doall (map (comp dna-melting #(ultimate-correction % ?-3)) result-?-3))
+(defn subsample
+"subsample v every n points"
+[n v]
+(vec (map (vec v) (vec (range 0 (count (vec v)) n)))))
+(defn max-difference [v] (- (apply max v) (apply min v)))
+(defn alignment [#^dna-data d]
+(let [distribution (vec (vals
+			  (apply sorted-map
+				 (interleave (:temp d) (:glow d)))))]
+(mean (map variance (partition 10 distribution)))))
+(defn garigh [d]
+(let [distribution (vec (vals
+			   (apply sorted-map
+				  (interleave (:temp d) (:glow d)))))]
+(view (visual (vec  (map max-difference (partition 100 distribution)))))))
+(defn alignment-4 [#^dna-data d]
+(let [distribution (vec (vals
+			  (apply sorted-map
+				 (interleave (:temp d) (:glow d)))))]
+(reduce + (map sq (map max-difference (partition 100 distribution))))))
+(defn alignment-2 [#^dna-data d]
+(let [distribution
+	(vec
+	 (vals
+	  (apply sorted-map
+		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
+(reduce +     (map max-difference (partition 30 distribution)))))
+(defn alignment-3 [#^dna-data d]
+(let [distribution
+	(vec
+	 (vals
+	  (apply sorted-map
+		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
+(reduce max     (map max-difference (partition 30 distribution)))))
+(comment
+;; this takes a while
+(def opts (zipmap all-samples (map #(last (take 20 (iterate (partial improve %) start-values))) all-samples))))
+;(def opts (zipmap all-samples opts-vals))
+(defn mean-map [v s]
+(map mean (map (partial vector v) s)))
+(def start-values
+{:test-k-sample-air [0.000001 0.001]
+:test-k-block-sample   [0.001 0.03]
+:test-room-temp  [290 300]
+:test-sample-temp  [285 310]
+:test-bleach [0.00001 0.0009]})
+(def start-values-2
+{:test-k-sample-air [1e-15 0.0005]
+:test-k-block-sample   [0.001 0.002]
+:test-room-temp  [290 300]
+:test-sample-temp  [290 300]
+:test-bleach [1e-7 0.0005]})
+(def start-values-3
+{:test-k-block-sample   [0.0005  0.00125 0.002]
+:test-k-sample-air [0.0000001 0.000001 0.0000019]
+:test-sample-temp  [280 300]
+:test-room-temp  [300 306]
+:test-bleach [0.000001 0.00001]})
+(defn condense [[a b :as v]]
+(if (< (/ (abs (- (apply max v) (apply min v))) (apply max v)) 1e-2) (vector (mean v)) v))
+(defn mass-shift
+[value v]
+(let [frac 2/3]
+(cond (>= (count v) 2)
+	  (let [[a b] v]
+	    (condense (vec (map #(float (+ (* frac %) (* (- 1 frac) value))) v))))
+	  (= (count v) 1)
+	  v)))
+(defn improve* [#^dna-data d
+	       {:keys [test-k-sample-air
+		       test-k-block-sample
+		       test-room-temp
+		       test-sample-temp
+		       test-bleach] :as domain}]
+(do (print \newline)
+(reduce (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
+	      (for [k-sample-air (condense test-k-sample-air)
+		    k-block-sample (condense test-k-block-sample)
+		    room-temp (condense test-room-temp)
+		    initial-sample-temp (condense test-sample-temp)
+		    bleach-constant (condense test-bleach)]
+		(do
+		  (print ".")
+		  (let [align
+			(alignment-2
+			 (ultimate-correction
+			  k-block-sample k-sample-air initial-sample-temp room-temp
+			  bleach-constant d))
+			align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]
+		    {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
+		     :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
+		     :test-room-temp (mass-shift room-temp test-room-temp)
+		     :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
+		     :test-bleach (mass-shift bleach-constant test-bleach)
+		     :align align}))))))
+(def println-repl (bound-fn [& args] (apply println args)))
+(def print-repl (bound-fn [& args] (apply print args)))
+(defn improve [#^dna-data d
+	       {:keys [test-k-block-sample
+		       test-k-sample-air
+		       test-sample-temp
+		       test-room-temp
+		       test-bleach] :as domain}]
+(reduce
+(fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
+(let [return-hash
+	     (fn
+	       [[k-block-sample k-sample-air initial-sample-temp
+		room-temp bleach-constant]]
+	       (do
+		 (print-repl ".")
+		 (let [align
+		       (alignment-4
+			(ultimate-correction
+			 k-block-sample k-sample-air initial-sample-temp room-temp
+			 bleach-constant d))
+		       align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]
+		   {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
+		    :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
+		    :test-room-temp (mass-shift room-temp test-room-temp)
+		    :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
+		    :test-bleach (mass-shift bleach-constant test-bleach)
+		    :align align})))]
+	 (pmap return-hash
+	       (cartesian-product
+		 test-k-block-sample
+		 test-k-sample-air
+		 test-sample-temp
+		 test-room-temp
+		 test-bleach)))))
+(def improve (memoize improve))
+(defn iterate-test [n #^dna-data d]
+(last (take n (iterate (partial improve d) start-values-3))))
+(defn nested-vals [key coll]
+(for [x (tree-seq coll? seq coll) :when (contains? x key)]
+(get x key)))
+(defn quantum-efficiency* [q-constant #^dna-data d]
+;; 0.8 at room temperature
+;; quantum efficiency is somehow
+;; inversely proportional to temperature
+(let [correction
+	(fn [t g] (/ g   (* 0.8 (exp (* q-constant (- 292 t))))))]
+(assoc d
+	:name (str (:name d) " |quantum-correction| ")
+	:glow (map correction (:temp d) (:glow d)))))
+(def quantum-efficiency (partial quantum-efficiency* 1.2e-2))
+(defn exp-decay [[A] x]
+(* 0.8 (exp (* A (- 298 x)))))
+;(defn final-join [#^dna-data d]
+(def join-opts
+{
+A-1 {:test-k-sample-air [1.760184242508937E-6], :test-k-block-sample [0.0012930832259977858], :test-room-temp [304.41977945963544], :test-sample-temp [281.31689453125], :test-bleach [9.956068424799014E-6], :align 0.800815268639094}
+A-2 {:test-k-sample-air [3.017639376897326E-7], :test-k-block-sample [0.0011550507430608075], :test-room-temp [304.6666717529297], :test-sample-temp [282.6337890625], :test-bleach [9.5610075732111E-6], :align 0.8530570499185829}
+A-3 {:test-k-sample-air [2.6829059556281816E-6], :test-k-block-sample [0.002413294818252325], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6]}
+B-2 {:test-k-sample-air [1.6274943088016396E-6], :test-k-block-sample [0.0013714926317334175], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7631166194402621}
+B-3
+{:test-k-sample-air [1.455345833771086E-6], :test-k-block-sample [0.0014841814991086721], :test-room-temp [304.6666717529297], :test-sample-temp [287.9835662841797], :test-bleach [9.965317076421343E-6], :align 0.5358800697386289}
+C-1
+{:test-k-sample-air [1.6907095717518434E-6], :test-k-block-sample [0.001278722076676786], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7249279378017799}
+C-2
+{:test-k-sample-air [1.77524979487013E-6], :test-k-block-sample [0.001278722076676786], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 2.0474987460209793}
+C-3
+{:test-k-sample-air [1.010541116860016E-6], :test-k-block-sample [0.001601272417853276], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7191738524207253}
+?-1
+{:test-k-sample-air [1.5533593114014366E-6], :test-k-block-sample [0.0013761443551629782], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 1.3838506291264723}
+?-3
+{:test-k-sample-air [2.940294715851148E-7], :test-k-block-sample [0.0018434864080821474], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.48549623932134006}
+})
+(defn final-join [#^dna-data d]
+(let [
+	m (apply sorted-map (interleave (:temp d) (:glow d)))
+	distribution (vec (vals m))
+	ordinates (vec (keys m))
+			   ]
+(assoc d
+:name (str (:name d) " |joined| ")
+:temp (map (fn [region] (mean [(apply min region) (apply max region)]))
+		 (partition 30 ordinates))
+:glow (map (fn [region] (mean [(apply min region) (apply max region)]))
+		 (partition 30 distribution)))))
+(defn melting-point [#^dna-data d]
+(let [t (drop 150 (:temp d))
+	g (drop 150 (:glow d))]
+(nth t (find-min (map / (diff g) (diff t) )))))
+(def target-dir (file-str "/home/r/MIT/6/6.122/DNA/output"))
+(defn x-form-dna [#^dna-data d]
+(final-join (ultimate-correction (join-opts d) d)))
+(defn save-raw-dna-curve [#^dna-data d]
+(save (dna-melting d) (str target-dir File/separator (:name d) ".png")))
+(defn save-dna-curve-with-temp [#^dna-data d*]
+(let [
+	melt (melting-point d*)
+	d (assoc d* :name (str (:name d*) " |M.P. = " (turn melt) " |"))
+	]
+(save (add-lines (dna-melting d) (repeat (count (:glow d)) melt) (:glow d))
+	  (str target-dir File/separator (:name d) ".png"))))
+(defn process-all []
+(doall (pmap  save-dna-curve all-samples))
+(doall (pmap (comp save-dna-curve-with-temp x-form-dna) all-samples)))
+(def ?-temps (vec (map (comp melting-point x-form-dna) [?-1 ?-3])))
+(def A-temps (vec (map (comp melting-point x-form-dna) [A-1 A-2 A-3])))
+(def B-temps (vec (map (comp melting-point x-form-dna) [B-2 B-3])))
+(def C-temps (vec (map (comp melting-point x-form-dna) [C-1 C-2 C-3])))
+(defn report [temps]
+(sorted-map :data temps :std (sqrt (variance temps)) :mean (mean temps)))
+;;; now it's time to start data reporting
+;(defn process [#^dna-data d]
+;  (ultimate-correction (opts d) d)
+(report A-temps)
+{:data [348.8771481779248 349.3989434763391 348.76079920968573],
+:mean 349.0122969546499,
+:std 0.33986161910548834}
+(report B-temps)
+{:data [349.0655539116902 349.92155530431506],
+:mean 349.4935546080026,
+:std 0.6052843894485812}
+(report C-temps)
+{:data [349.92341967511976 349.68832687337056 348.2872189461977],
+:mean 349.299655164896,
+:std 0.884639745344559}
+(report ?-temps)
+{:data [350.4525305192239 349.87377093318975],
+:mean 350.1631507262068,
+:std 0.4092448279511269}

Mercurial > rlm

comparison src/rlm/dna_melting.clj @ 0:78a630e650d2