Point estimate

bin/VERSION

@@ -1 +1 @@
-v0.11
+v0.11.1

bin/whippet-delta.jl

@@ -45,6 +45,13 @@ function parse_cmd()
       help     = "Calculate max probability of |deltaPsi| greater than this value (default is 0.0, it is not advisable to change this)."
       arg_type = Float64
       default  = 0.0 =#
+    "--pseudo-adjust"
+      help     = "When a single replicate is given, pseudo-adjust a \"second\" for mle fitting by this value."
+      arg_type = Float64
+      default  = 0.001
+    "--use-depth"
+      help     = "Sample each replicate's PSI according to read-depth of the event. (default off)"
+      action   = :store_true
     "--emperical-size", "-e"
       help     = "Emperical distribution size to sample from."
       arg_type = Int64
@@ -92,7 +99,9 @@ function main()
                              min_samp=args["min-samples"], 
                              min_reads=args["min-reads"],
                              amt=0.0,
-                             size=args["emperical-size"] ) 
+                             size=args["emperical-size"],
+                             point_est=!args["use-depth"], 
+                             pseudo_adj=args["pseudo-adjust"]) 
    println(STDERR, "Whippet $ver done." )
 end
 

src/diff.jl

@@ -21,11 +21,24 @@ end
 # combine multiple beta posteriors for replicates into one
 # if unpaired, we don't need to maintain sample1_rep1 to sample2_rep1
 # so we can re-sample from the new joint beta
-function PosteriorPsi( set::Vector{PosteriorPsi}; paired=false, size=1000 )
-   emperical = vcat( map( x->x.emperical, set )... )
+function PosteriorPsi( set::Vector{PosteriorPsi}; 
+                       paired::Bool=false, size::Int=1000,
+                       point_est::Bool=true, pseudo_adj::Float64=0.001)
+   if point_est
+      psiarr = [x.psi for x in set]
+      if length(set) > 1 && var(psiarr) > 0.0
+         emperical = [x.psi for x in set]
+      else
+         one = psiarr[1]
+         two = one < (1.0 - pseudo_adj) ? one + pseudo_adj : one - pseudo_adj
+         emperical = vcat(psiarr, two)
+      end
+   else
+      emperical = vcat( map( x->x.emperical, set )... )
+   end
    beta = fit(Beta, emperical)
    psi  = mean(beta)
-   if !paired
+   if !paired || point_est
       emperical = rand(beta, size)
    end
    PosteriorPsi( beta, emperical, psi )
@@ -84,7 +97,8 @@ end
 
 parse_complexity( c::S ) where {S <: AbstractString} = split( c, COMPLEX_CHAR, keep=false )[1] |> x->parse(Int,x)
 
-function process_psi_line( streams::Vector{BufferedStreams.BufferedInputStream}; min_reads=5, size=1000 )
+function process_psi_line( streams::Vector{BufferedStreams.BufferedInputStream}; 
+                           min_reads=5, size=1000 )
    postvec = Vector{PosteriorPsi}()
    event   = split( "", "" )
    complex = 0
@@ -110,7 +124,9 @@ end
 
 function process_psi_files( outfile, a::Vector{BufferedStreams.BufferedInputStream}, 
                                      b::Vector{BufferedStreams.BufferedInputStream}; 
-                                     min_samp=1, min_reads=5, amt=0.0, size=1000 )
+                                     min_samp::Int=1, min_reads::Int=5, 
+                                     amt::Float64=0.0, size::Int=1000,
+                                     point_est::Bool=true, pseudo_adj::Float64=0.001 )
    io = open( outfile, "w" )
    stream = ZlibDeflateOutputStream( io )
    output_diff_header( stream )
@@ -124,8 +140,8 @@ function process_psi_files( outfile, a::Vector{BufferedStreams.BufferedInputStre
       entropy = a_entropy > b_entropy ? a_entropy : b_entropy
       @assert( a_event == b_event, "Incorrect events matched!!" )
       if length(a_post) >= min_samp && length(b_post) >= min_samp
-         a_post   = PosteriorPsi( a_post ) # fit new posterior
-         b_post   = PosteriorPsi( b_post ) 
+         a_post   = PosteriorPsi( a_post, point_est=point_est, pseudo_adj=pseudo_adj ) # fit new posterior
+         b_post   = PosteriorPsi( b_post, point_est=point_est, pseudo_adj=pseudo_adj ) 
          fwdprob  = probability( a_post, b_post, amt=amt )
          revprob  = probability( b_post, a_post, amt=amt )
          prob     = fwdprob > revprob ? fwdprob : revprob

test/runtests.jl

@@ -613,5 +613,9 @@ IIIIIIIIIIII
          @test cnt == 128482 # correct number of reads in file
          #run(`julia ../bin/whippet-quant.jl --ebi SRR1199010`)
       end
+
+      @testset "Diff and Delta" begin
+
+      end
    end
 end