diff --git a/bin/VERSION b/bin/VERSION index 5416288..4f7638f 100644 --- a/bin/VERSION +++ b/bin/VERSION @@ -1 +1 @@ -v0.11 +v0.11.1 diff --git a/bin/whippet-delta.jl b/bin/whippet-delta.jl index 6675e75..1414484 100755 --- a/bin/whippet-delta.jl +++ b/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 diff --git a/src/diff.jl b/src/diff.jl index 65cc12d..4601b1c 100644 --- a/src/diff.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index f5d1605..47fa95e 100644 --- a/test/runtests.jl +++ b/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