[BioRuby-cvs] bioruby/lib/bio/util/restriction_enzyme/analysis calculated_cuts.rb, NONE, 1.1 cut_range.rb, NONE, 1.1 cut_ranges.rb, NONE, 1.1 fragment.rb, NONE, 1.1 fragments.rb, NONE, 1.1 horizontal_cut_range.rb, NONE, 1.1 permutation.rb, NONE, 1.1 sequence_range.rb, NONE, 1.1 tags.rb, NONE, 1.1 vertical_cut_range.rb, NONE, 1.1

Trevor Wennblom trevor at pub.open-bio.org
Wed Feb 1 07:27:27 UTC 2006


Update of /home/repository/bioruby/bioruby/lib/bio/util/restriction_enzyme/analysis
In directory pub.open-bio.org:/tmp/cvs-serv28844/analysis

Added Files:
	calculated_cuts.rb cut_range.rb cut_ranges.rb fragment.rb 
	fragments.rb horizontal_cut_range.rb permutation.rb 
	sequence_range.rb tags.rb vertical_cut_range.rb 
Log Message:
Bio::RestrictionEnzyme


--- NEW FILE: vertical_cut_range.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio/util/restriction_enzyme/analysis/cut_range'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/vertical_cut_range.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: vertical_cut_range.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/vertical_cut_range.rb -
=end
class VerticalCutRange < CutRange
  attr_reader :p_cut_left, :p_cut_right
  attr_reader :c_cut_left, :c_cut_right
  attr_reader :min, :max
  attr_reader :range

  def initialize( p_cut_left=nil, p_cut_right=nil, c_cut_left=nil, c_cut_right=nil )
    @p_cut_left = p_cut_left
    @p_cut_right = p_cut_right
    @c_cut_left = c_cut_left
    @c_cut_right = c_cut_right

    a = [@p_cut_left, @c_cut_left, @p_cut_right, @c_cut_right]
    a.delete(nil)
    a.sort!
    @min = a.first
    @max = a.last

    @range = nil
    @range = (@min.. at max) unless @min == nil or @max == nil
  end

  def include?(i)
    return false if @range == nil
    @range.include?(i)
  end

end

end
end

--- NEW FILE: tags.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

#require 'bio/util/restriction_enzyme/analysis/shared_information'
require 'bio'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/tags.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: tags.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/tags.rb -
=end
class Tags < Hash
end

end
end

--- NEW FILE: cut_range.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/cut_range.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: cut_range.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/cut_range.rb -
=end
class CutRange
end

end
end

--- NEW FILE: permutation.rb ---
# = permutation.rb - Permutation class for Ruby
#
# == Author
#
# Florian Frank mailto:flori at ping.de
#
# == License
#
# This is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License Version 2 as published by the Free
# Software Foundation: www.gnu.org/copyleft/gpl.html
#
# == Download
#
# The latest version of <b>permutation</b> can be found at
#
# * http://rubyforge.org/frs/?group_id=291
#
# The homepage of this library is located at
#
# * http://permutation.rubyforge.org
#
# == Description
#
# This class has a dual purpose: It can be used to create permutations
# of a given size and to do some simple computations with/on
# permutations. The instances of this class don't require much memory
# because they don't include the permutation as a data structure. They
# only save the information necessary to create the permutation if asked
# to do so.
#
# To generate permutations the ranking/unranking method described in [SS97]
# is used. Because of Ruby's Bignum arithmetic it is useful also
# for permutations of very large size.
# 
# == Examples
#
# In this section some examples show what can be done with this class.
#
# Creating all permutations and project them on data:
#
#  perm = Permutation.new(3)
#  # => #<Permutation:0x57dc94 @last=5, @rank=0, @size=3>
#  perm.map { |p| p.value }
#  # => [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]
#  colors = [:r, :g, :b]
#  # => [:r, :g, :b]
#  perm.map { |p| p.project(colors) }
#  # => [[:r, :g, :b], [:r, :b, :g], [:g, :r, :b], [:g, :b, :r], [:b, :r, :g],
#  #    [:b, :g, :r]]
#  string = "abc"# => "abc"
#  perm.map { |p| p.project(string) }
#  # => ["abc", "acb", "bac", "bca", "cab", "cba"]
#
# Or perhaps more convenient to use:
#
#  perm = Permutation.for("abc")
#  perm.map { |p| p.project }
#  # => ["abc", "acb", "bac", "bca", "cab", "cba"]
#
# Finding the successor and predecessor of Permutations or a
# certain Permutation for a given rank:
#
#  perm = Permutation.new(7)
#  # => #<Permutation:0x8453c @rank=0, @size=7, @last=5039>
#  perm.succ!
#  # => #<Permutation:0x8453c @rank=1, @size=7, @last=5039>
#  perm.succ!
#  # => #<Permutation:0x8453c @rank=2, @size=7, @last=5039>
#  perm.succ!
#  # => #<Permutation:0x8453c @rank=3, @size=7, @last=5039>
#  perm.pred!
#  # => #<Permutation:0x8453c @rank=2, @size=7, @last=5039>
#  perm.rank = 3200
#  # => 3200
#  perm
#  # => #<Permutation:0x8453c @rank=3200, @size=7, @last=5039>
#  perm.value
#  # => [4, 2, 5, 1, 3, 0, 6]
#
# Generating random Permutations
#
#  perm = Permutation.new(10)
#  # => #<Permutation:0x59f4c0 @rank=0, @size=10, @last=3628799>
#  perm.random!.value
#  # => [6, 4, 9, 7, 3, 5, 8, 1, 2, 0]
#  perm.random!.value
#  # => [3, 7, 6, 1, 4, 8, 9, 2, 5, 0]
#  perm.random!.value
#  # => [2, 8, 4, 9, 3, 5, 6, 7, 0, 1]
#  perm.random!.project("ABCDEFGHIJ")
#  # => "DFJGAEBCIH"
#  perm.random!.project("ABCDEFGHIJ")
#  # => "BFADEGHJCI"
#
# Performing some mathematical operations on/with Permutations
#
#  p1 = Permutation.from_cycles([[1, 3, 2], [5, 7]], 10)
#  # => #<Permutation:0x593594 @rank=80694, @size=10, @last=3628799>
#  p2 = Permutation.from_value [3, 2, 0, 5, 6, 8, 9, 1, 4, 7]
#  # => #<Permutation:0x5897b0 @rank=1171050, @size=10, @last=3628799>
#  p3 = p1 * p2
#  # => #<Permutation:0x586a88 @rank=769410, @size=10, @last=3628799>
#  p3.value
#  # => [2, 1, 0, 7, 6, 8, 9, 3, 4, 5]
#  p3.cycles
#  # => [[0, 2], [3, 7], [4, 6, 9, 5, 8]]
#  p4 = p1 * -p2
#  # => #<Permutation:0x581a10 @rank=534725, @size=10, @last=3628799>
#  p4.value
#  # => [1, 5, 3, 0, 8, 2, 4, 9, 7, 6]
#  p4.cycles
#  # => [[0, 1, 5, 2, 3], [4, 8, 7, 9, 6]]
#  id = p1 * -p1
#  # => #<Permutation:0x583a7c @rank=0, @size=10, @last=3628799>
#
# == References
#
#  [SS97] The Algorithm Design Manual, Steven S. Skiena, Telos/Springer, 1997.
# 

class Permutation

    include Enumerable
    include Comparable

    # Creates a new Permutation instance of <code>size</code>
    # (and ranked with <code>rank</code>).
    def initialize(size, rank = 0)
        @size, @rank = size, rank
        @last = factorial(size) - 1
    end

    # Creates a new Permutation instance from the Array
    # <code>indices</code>, that should consist of a permutation of Fixnums
    # in the range of <code>0</code> and <code>indices.size - 1</code>. This is
    # for example the result of a call to the Permutation#value method.
    def self.from_value(indices)
        obj = new(indices.size)
        obj.instance_eval do
            self.rank = rank_indices(indices)
        end
        obj
    end

    # Creates a new Permutation instance from the Array of Arrays
    # <code>cycles</code>. This is for example the result of a
    # call to the Permutation#cycles method .
    def self.from_cycles(cycles, max = 0)
        indices = Array.new(max)
        cycles.each do |cycle|
            cycle.empty? and next
            for i in 0...cycle.size
                indices[ cycle[i - 1] ] = cycle[i]
            end
        end
        indices.each_with_index { |r, i| r or indices[i] = i }
        from_value(indices)
    end

  # A permutation instance of size collection.size is created with
  # collection as the default Permutation#project data object. A
  # collection should respond to size, [], and []=. The Permutation
  # instance will default to rank 0 if none is given.
  def self.for(collection, rank = 0)
    perm = new(collection.size, rank)
    perm.instance_variable_set(:@collection, collection)
    perm
  end

    # Returns the size of this permutation, a Fixnum.
    attr_reader :size

    # Returns the size of this permutation, a Fixnum in the range
    # of 0 and Permutation#last.
    attr_reader :rank

    # Returns the rank of the last ranked Permutation of size
    # Permutation#size .
    attr_reader :last

    # Assigns <code>m</code> to the rank attribute of this Permutation
    # instance. That implies that the indices produced by a call to the
    # Permutation#value method of this instance is the permutation ranked with
    # this new <code>rank</code>.
    def rank=(m)
        last = factorial(size) - 1
        while m > last do m -= last end
        while m < 0 do m += last end
        @rank = m
    end

    # Returns the indices in the range of 0    to Permutation#size - 1
    # of this permutation that is ranked with Permutation#rank.
    #
    # <b>Example:</b>
    #  perm = Permutation.new(6, 312)
    #  # => #<Permutation:0x6ae34 @last=719, @rank=312, @size=6>
    #  perm.value            
    #  # => [2, 4, 0, 1, 3, 5]
    def value
        unrank_indices(@rank)
    end

    # Returns the projection of this instance's Permutation#value
    # into the <code>data</code> object that should respond to
    # the #[] method. If this Permutation inbstance was created
    # with Permutation.for the collection used to create
    # it is used as a data object.
    #
    # <b>Example:</b>
    #  perm = Permutation.new(6, 312)
    #  # => #<Permutation:0x6ae34 @last=719, @rank=312, @size=6>
    #  perm.project("abcdef")
    #  # => "ceabdf"
    def project(data = @collection)
        data or raise ArgumentError.new("a collection is required to project")
        raise ArgumentError.new("data size is != #{size}!") if data.size != size
        projection = data.clone
        value.each_with_index { |i, j| projection[j] = data[i] }
        projection
    end

    # Switches this instances to the next ranked Permutation.
    # If this was the Permutation#last permutation it wraps around
    # the first (<code>rank == 0</code>) permutation.
    def next!
        @rank += 1
        last = factorial(size) - 1
        @rank = 0 if @rank > last
        self
    end

    alias succ! next!

    # Returns the next ranked Permutation instance.
    # If this instance is the Permutation#last permutation it returns the first
    # (<code>rank == 0</code>) permutation.
    def next
        clone.next!
    end

    alias succ next

    # Switches this instances to the previously ranked Permutation.
    # If this was the first permutation it returns the last (<code>rank ==
    # Permutation#last</code>) permutation.
    def pred!
        @rank -= 1
        last = factorial(size) - 1
        @rank = last if @rank < 0
        self
    end

    # Returns the previously ranked Permutation. If this was the first
    # permutation it returns the last (<code>rank == Permutation#last</code>)
    # permutation.
    def pred
        clone.pred!
    end

    # Switches this Permutation instance to random permutation
    # of size Permutation#size.
    def random!
        new_rank = rand(last).to_i
        self.rank = new_rank
        self
    end

    # Returns a random Permutation instance # of size Permutation#size.
    def random
        clone.random!
    end

    # Iterates over all permutations of size Permutation#size starting with the
    # first (<code>rank == 0</code>) ranked permutation and ending with the
    # last (<code>rank == Permutation#last</code>) ranked permutation while
    # yielding to a freshly created Permutation instance for every iteration
    # step.
    #
    # The mixed in methods from the Enumerable module rely on this method.
    def each # :yields: perm
        0.upto(last) do |r|
            klon = clone
            klon.rank = r
            yield klon
        end
    end

    # Does something similar to Permutation#each. It doesn't create new
    # instances (less overhead) for every iteration step, but yields to a
    # modified self instead. This is useful if one only wants to call a
    # method on the yielded value and work with the result of this call. It's
    # not a good idea to put the yielded values in a data structure because the
    # will all reference the same (this!) instance. If you want to do this
    # use Permutation#each.
    def each!
        old_rank = rank
        0.upto(last) do |r|
            self.rank = r
            yield self
        end
        self.rank = old_rank
    end

    # Compares to Permutation instances according to their Permutation#size
    # and the Permutation#rank.
    #
    # The mixed in methods from the Comparable module rely on this method.
    def <=>(other)
        size <=> other.size.zero? || rank <=> other.rank
    end

    # Returns true if this Permutation instance and the other have the same
    # value, that is both Permutation instances have the same Permutation#size
    # and the same Permutation#rank.
    def eql?(other)
        self.class == other.class && size == other.size && rank == other.rank
    end

    alias == eql?

    # Computes a unique hash value for this Permutation instance.
    def hash
        size.hash ^ rank.hash
    end

    # Switchtes this Permutation instance to the inverted permutation.
    # (See Permutation#compose for an example.)
    def invert!
        indices = unrank_indices(rank)
        inverted = Array.new(size)
        for i in 0...size
            inverted[indices[i]] = i
        end
        self.rank = rank_indices(inverted)
        self
    end

    # Returns the inverted Permutation of this Permutation instance.
    # (See Permutation#compose for an example.)
    def invert
        clone.invert!
    end

    alias -@ invert

    # Compose this Permutation instance and the other to
    # a new Permutation. Note that a permutation
    # composed with it's inverted permutation yields
    # the identity permutation, the permutation with rank 0.
    #
    # <b>Example:</b>
    #  p1 = Permutation.new(5, 42)
    #  # => #<Permutation:0x75370 @last=119, @rank=42, @size=5>
    #  p2 = p1.invert
    #  # => #<Permutation:0x653d0 @last=119, @rank=51, @size=5>
    #  p1.compose(p2)
    #  => #<Permutation:0x639a4 @last=119, @rank=0, @size=5>
    # Or a little nicer to look at:
    #  p1 * -p1
    #  # => #<Permutation:0x62004 @last=119, @rank=0, @size=5>
    def compose(other)
        size == other.size or raise ArgumentError.new(
            "permutations of unequal sizes cannot be composed!")
        indices = self.value
        composed = other.value.map { |i| indices[i] }
        klon = clone
        klon.rank = rank_indices(composed)
        klon
    end

    alias * compose

    # Returns the cycles representation of this Permutation instance.
    # The return value of this method can be used to create a
    # new Permutation instance with the Permutation.from_cycles method.
    #
    # <b>Example:</b>
    #  perm = Permutation.new(7, 23)
    #  # => #<Permutation:0x58541c @last=5039, @rank=23, @size=7>
    #  perm.cycles
    #  # => [[3, 6], [4, 5]]
    def cycles
        perm = value
        result = [[]]
        seen = {}
        current = nil
        until seen == perm.size
            current or current = perm.find { |x| !seen[x] }
            break unless current
            if seen[current]
                current = nil
                result << []
            else
                seen[current] = true
                result[-1] << current
                current = perm[current]
            end
        end
        result.pop
        result.select { |c| c.size > 1 }.map do |c|
            min_index = c.index(c.min)
            c[min_index..-1] + c[0...min_index]
        end
    end

    # Returns the signum of this Permutation instance.
    # It's -1 if this permutation is odd and 1 if it's
    # an even permutation.
    #
    # A permutation is odd if it can be represented by an odd number of
    # transpositions (cycles of length 2), or even if it can be represented of
    # an even number of transpositions.
    def signum
        s = 1
        cycles.each do |c|
            c.size % 2 == 0 and s *= -1
        end
        s
    end

    alias sgn signum

    # Returns true if this permutation is even, false otherwise.
    def even?
        signum == 1
    end

    # Returns true if this permutation is odd, false otherwise.
    def odd?
        signum == -1
    end

    private

    @@factorial_cache = {}

    def factorial(n)
        f = @@factorial_cache[n] and return f
        f = 1
        for i in 2..n do f *= i end
        @@factorial_cache[n] = f
    end

    def rank_indices(p)
        result = 0
        for i in 0...size
            result += p[i] * factorial(size - i - 1)
            for j in (i + 1)...size
                p[j] -= 1 if p[j] > p[i] 
            end
        end
        result
    end

    def unrank_indices(m)
        result = Array.new(size, 0)
        for i in 0...size
            f = factorial(i)
            x = m % (f * (i + 1))
            m -= x
            x /= f
            result[size - i - 1] = x
            x -= 1
            for j in (size - i)...size
                result[j] += 1 if result[j] > x
            end
        end
        result
    end

end

--- NEW FILE: fragment.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio/util/restriction_enzyme/analysis/tags'
require 'bio/util/restriction_enzyme/analysis/cut_ranges'
require 'bio/util/restriction_enzyme/analysis/horizontal_cut_range'
require 'bio'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/fragment.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: fragment.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/fragment.rb -
=end
class Fragment

  attr_reader :size
#  attr_reader :tags

  def initialize( primary_bin, complement_bin )
#    @tags = []
    @primary_bin = primary_bin
    @complement_bin = complement_bin
  end

#  def add_tag( index, info=nil )
#    @tags[index] = info
#  end

  DisplayFragment = Struct.new(:primary, :complement)

  def for_display(p_str=nil, c_str=nil)
    df = DisplayFragment.new
    df.primary = ''
    df.complement = ''

    both_bins = (@primary_bin + @complement_bin).sort.uniq
    both_bins.each do |item|
      @primary_bin.include?(item) ? df.primary << p_str[item] : df.primary << ' '
      @complement_bin.include?(item) ? df.complement << c_str[item] : df.complement << ' '
    end

    df
  end



end

end
end

--- NEW FILE: horizontal_cut_range.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio/util/restriction_enzyme/analysis/cut_range'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/horizontal_cut_range.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: horizontal_cut_range.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/horizontal_cut_range.rb -
=end
class HorizontalCutRange < CutRange
  attr_reader :p_cut_left, :p_cut_right
  attr_reader :c_cut_left, :c_cut_right
  attr_reader :min, :max
  attr_reader :hcuts

  def initialize( left, right=left )
    raise "left > right" if left > right

    # The 'range' here is actually off by one on the left
    # side in relation to a normal CutRange, so using the normal
    # variables from CutRange would result in unpredictable
    # behavior.

    @p_cut_left = nil
    @p_cut_right = nil
    @c_cut_left = nil
    @c_cut_right = nil
    @min = nil
    @max = nil
    @range = nil

    @hcuts = (left..right)
  end

  def include?(i)
    @range.include?(i)
  end

end

end
end

--- NEW FILE: cut_ranges.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/cut_ranges.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: cut_ranges.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/cut_ranges.rb -
=end
class CutRanges < Array
  def min; self.collect{|a| a.min}.flatten.sort.first; end
  def max; self.collect{|a| a.max}.flatten.sort.last; end
  def include?(i); self.collect{|a| a.include?(i)}.include?(true); end
end

end
end

--- NEW FILE: fragments.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/fragments.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: fragments.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/fragments.rb -
=end
class Fragments < Array
  
  attr_accessor :primary
  attr_accessor :complement

  def initialize(primary, complement)
    @primary = primary
    @complement = complement
  end

  DisplayFragment = Struct.new(:primary, :complement)

  def for_display(p_str=nil, c_str=nil)
    p_str ||= @primary
    c_str ||= @complement
    pretty_fragments = []
    self.each { |fragment| pretty_fragments << fragment.for_display(p_str, c_str) }
    pretty_fragments
  end

end

end
end

--- NEW FILE: sequence_range.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio/util/restriction_enzyme/analysis/tags'
require 'bio/util/restriction_enzyme/analysis/cut_ranges'
require 'bio/util/restriction_enzyme/analysis/horizontal_cut_range'
require 'bio/util/restriction_enzyme/analysis/vertical_cut_range'
require 'bio/util/restriction_enzyme/analysis/calculated_cuts'
require 'bio/util/restriction_enzyme/analysis/fragments'
require 'bio/util/restriction_enzyme/analysis/fragment'
require 'bio'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/sequence_range.rb -
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: sequence_range.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/sequence_range.rb -
=end
class SequenceRange

  attr_reader :p_left, :p_right
  attr_reader :c_left, :c_right

  attr_reader :left, :right
  attr_reader :size
  attr_reader :tags
  attr_reader :cut_ranges

  def initialize( p_left = nil, p_right = nil, c_left = nil, c_right = nil )
    @__fragments_current = false
    raise ArgumentError if p_left == nil and c_left == nil
    raise ArgumentError if p_right == nil and c_right == nil
    (raise ArgumentError unless p_left <= p_right) unless p_left == nil or p_right == nil
    (raise ArgumentError unless c_left <= c_right) unless c_left == nil or c_right == nil

    @p_left  = p_left
    @p_right = p_right
    @c_left  = c_left
    @c_right = c_right

    tmp = [p_left, c_left]
    tmp.delete(nil)
    @left = tmp.sort.first

    tmp = [p_right, c_right]
    tmp.delete(nil)
    @right = tmp.sort.last

    @size = (@right - @left) + 1 unless @left == nil or @right == nil

    @tags = Tags.new
    @cut_ranges = CutRanges.new
  end

=begin
Special Case: Horizontal cuts at beginning or end of strand
=end

  Bin = Struct.new(:c, :p)

  def fragments
    return @__fragments if @__fragments_current == true
    @__fragments_current = true

    cc = CalculatedCuts.new(@size)
    cc.add_cuts_from_cut_ranges(@cut_ranges)
    cc.remove_incomplete_cuts

    p_cut = cc.vc_primary
    c_cut = cc.vc_complement
    h = cc.hc_between_strands 

    if @circular
    # NOTE
    # if it's circular we should start at the beginning of a cut for orientation
    # scan for it, hack off the first set of hcuts and move them to the back
    else
#      last_index = @size - 1
      p_cut.unshift(-1) unless p_cut.include?(-1)
#      p_cut.push(last_index) unless p_cut.include?(last_index)
      c_cut.unshift(-1) unless c_cut.include?(-1)
#      c_cut.push(last_index) unless c_cut.include?(last_index)
    end


    if @circular
      largest_bin = 0
    else
      largest_bin = -1
    end
    p_bin = largest_bin
    c_bin = largest_bin
    bins = { largest_bin => Bin.new }  # bin_id, bin
    bins[ largest_bin ].p = []
    bins[ largest_bin ].c = []

    x = lambda do |bin_id|
      largest_bin += 1
      bins[ bin_id ] = Bin.new
      bins[ bin_id ].p = []
      bins[ bin_id ].c = []
    end

    -1.upto(@size-1) do |idx|

      # if bins are out of sync but the strands are attached
      if p_bin != c_bin and h.include?(idx) == false
        bins.delete( [p_bin, c_bin].sort.last )
        p_bin = c_bin = [p_bin, c_bin].sort.first
        largest_bin -= 1
      end

      bins[ p_bin ].p << idx
      bins[ c_bin ].c << idx

      if p_cut.include? idx
        p_bin = largest_bin + 1
        x.call(p_bin)
      end

      if c_cut.include? idx
        c_bin = largest_bin + 1
        x.call(c_bin)
      end

    end

    # Easy way to indicate the start of a strand just in case
    # there is a horizontal cut at position 0
    bins.delete(-1) unless @circular

#    require 'pp'
#    pp bins

#NOTE
    str1 = nil
    str2 = nil

    num_txt_repeat = lambda { num_txt = '0123456789'; (num_txt * ( @size / num_txt.size.to_f ).ceil)[0.. at size-1] }
    (str1 == nil) ? a = num_txt_repeat.call : a = str1.dup
    (str2 == nil) ? b = num_txt_repeat.call : b = str2.dup

    fragments = Fragments.new(a,b)

    bins.sort.each do |k, bin|
      fragment = Fragment.new( bin.p, bin.c )
      @tags.each { |k,v| fragment.add_tag(k,v) if (ts.left..ts.right).include?(k) }
      fragments << fragment
    end

#pp fragments.for_display
#    pp fragments
#    exit

    @__fragments = fragments
    return fragments
  end
  
  def add_tag( index, info=nil )
    @__fragments_current = false

    raise IndexError unless index >= @left and index <= @right
    @tags[index] = info
  end

# Cut occurs immediately after the index supplied.
# For example, a cut at '0' would mean a cut occurs between 0 and 1.
  def add_cut_range( p_cut_left=nil, p_cut_right=nil, c_cut_left=nil, c_cut_right=nil )
    @__fragments_current = false

    if p_cut_left.kind_of? CutRange
      @cut_ranges << p_cut_left
    else
      (raise IndexError unless p_cut_left >= @left and p_cut_left <= @right) unless p_cut_left == nil
      (raise IndexError unless p_cut_right >= @left and p_cut_right <= @right) unless p_cut_right == nil
      (raise IndexError unless c_cut_left >= @left and c_cut_left <= @right) unless c_cut_left == nil
      (raise IndexError unless c_cut_right >= @left and c_cut_right <= @right) unless c_cut_right == nil

      @cut_ranges << VerticalCutRange.new( p_cut_left, p_cut_right, c_cut_left, c_cut_right )
    end
  end

  def add_cut_ranges(*cut_ranges)
    cut_ranges.flatten!
    cut_ranges.each do |cut_range|
      raise TypeError, "Not of type CutRange" unless cut_range.kind_of? CutRange
      self.add_cut_range( cut_range )
    end
  end

  def add_horizontal_cut_range( left, right=left )
    @__fragments_current = false
    @cut_ranges << HorizontalCutRange.new( left, right )
  end


end

end
end

--- NEW FILE: calculated_cuts.rb ---
require 'pathname'
libpath = Pathname.new(File.join(File.dirname(__FILE__), ['..'] * 5, 'lib')).cleanpath.to_s
$:.unshift(libpath) unless $:.include?(libpath)

require 'bio/util/restriction_enzyme/cut_symbol'
require 'bio/util/restriction_enzyme/string_formatting'

module Bio; end
class Bio::RestrictionEnzyme
class Analysis

#
# bio/util/restriction_enzyme/analysis/calculated_cuts.rb - 
#
# Copyright::  Copyright (C) 2006 Trevor Wennblom <trevor at corevx.com>
# License::    LGPL
#
#  $Id: calculated_cuts.rb,v 1.1 2006/02/01 07:34:11 trevor Exp $
#
#
#--
#
#  This library is free software; you can redistribute it and/or
#  modify it under the terms of the GNU Lesser General Public
#  License as published by the Free Software Foundation; either
#  version 2 of the License, or (at your option) any later version.
#
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with this library; if not, write to the Free Software
#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
#
#++
#
#

=begin rdoc
bio/util/restriction_enzyme/analysis/calculated_cuts.rb - 


   1 2 3 4 5 6 7
   G A|T T A C A
      +-----+
   C T A A T|G T
   1 2 3 4 5 6 7

Primary cut = 2
Complement cut = 5
Horizontal cuts = 3, 4, 5
=end
class CalculatedCuts
  include CutSymbol
  include StringFormatting

  # Vertical cuts on the primary strand
  attr_reader :vc_primary
  
  # Vertical cuts on the complement strand
  attr_reader :vc_complement

  # Horizontal cuts
  attr_reader :hc_between_strands

  # Set to +true+ if the fragment CalculatedCuts is working on is circular
  attr_accessor :circular

  def initialize(size=nil, circular=false)
    @size = size
    @circular = circular
    @vc_primary = []
    @vc_complement = []
    @hc_between_strands = []
  end

  def add_cuts_from_cut_ranges(cut_ranges)
    @strands_for_display_current = false

    cut_ranges.each do |cut_range|
      @vc_primary += [cut_range.p_cut_left, cut_range.p_cut_right]
      @vc_complement += [cut_range.c_cut_left, cut_range.c_cut_right]

      if cut_range.class == VerticalCutRange
        ( cut_range.min + 1 ).upto( cut_range.max ){|i| @hc_between_strands << i} if cut_range.min < cut_range.max
      elsif cut_range.class == HorizontalCutRange
        ( cut_range.hcuts.first ).upto( cut_range.hcuts.last ){|i| @hc_between_strands << i}
      end
    end
    clean_all
  end

  def remove_incomplete_cuts(size=nil)
    @strands_for_display_current = false
    @size = size if size
    raise IndexError, "Size of the strand must be provided here or during initalization." if !@size.kind_of?(Fixnum) and not @circular

    vcuts = (@vc_primary + @vc_complement).uniq.sort
    hcuts = @hc_between_strands
    last_index = @size - 1
    good_hcuts = []
    potential_hcuts = []

    if @circular
    # NOTE
    # if it's circular we should start at the beginning of a cut for orientation
    # scan for it, hack off the first set of hcuts and move them to the back
    else
      vcuts.unshift(-1) unless vcuts.include?(-1)
      vcuts.push(last_index) unless vcuts.include?(last_index)
    end

    hcuts.each do |hcut|
      raise IndexError if hcut < -1 or hcut > last_index
      # skipped a nucleotide
      potential_hcuts.clear if !potential_hcuts.empty? and (hcut - potential_hcuts.last).abs > 1

      if potential_hcuts.empty?
        if vcuts.include?( hcut ) and vcuts.include?( hcut - 1 )
          good_hcuts += [hcut]
        elsif vcuts.include?( hcut - 1 )
          potential_hcuts << hcut
        end
      else
        if vcuts.include?( hcut )
          good_hcuts += potential_hcuts + [hcut]
          potential_hcuts.clear
        else
          potential_hcuts << hcut
        end
      end
    end

    check_vc = lambda do |vertical_cuts, opposing_vcuts|
      # opposing_vcuts is here only to check for blunt cuts, so there shouldn't
      # be any out-of-order problems with this
      good_vc = []
      vertical_cuts.each { |vc| good_vc << vc if good_hcuts.include?( vc ) or good_hcuts.include?( vc + 1 ) or opposing_vcuts.include?( vc ) }
      good_vc
    end

    @vc_primary = check_vc.call(@vc_primary, @vc_complement)
    @vc_complement = check_vc.call(@vc_complement, @vc_primary)
    @hc_between_strands = good_hcuts

    clean_all
  end

  def strands_for_display(str1 = nil, str2 = nil, vcp=nil, vcc=nil, hc=nil)
    return @strands_for_display if @strands_for_display_current
    vcs = '|'
    hcs = '-'
    vhcs = '+'
      
    num_txt_repeat = lambda { num_txt = '0123456789'; (num_txt * ( @size / num_txt.size.to_f ).ceil)[0.. at size-1] }
    (str1 == nil) ? a = num_txt_repeat.call : a = str1.dup
    (str2 == nil) ? b = num_txt_repeat.call : b = str2.dup

    vcp = @vc_primary if vcp==nil
    vcc = @vc_complement if vcc==nil
    hc = @hc_between_strands if hc==nil

    vcuts = (vcp + vcc).uniq.sort

    vcp.reverse.each { |c| a.insert(c+1, vcs) }
    vcc.reverse.each { |c| b.insert(c+1, vcs) }

    between = ' ' * @size
    hc.each {|hcut| between[hcut,1] = hcs }

    s_a = add_spacing(a, vcs)
    s_b = add_spacing(b, vcs)
    s_bet = add_spacing(between)

    # NOTE watch this for circular
    i = 0
    0.upto( s_a.size-1 ) do
      if (s_a[i,1] == vcs) or (s_b[i,1] == vcs)
        s_bet[i] = vhcs 
      elsif i != 0 and s_bet[i-1,1] == hcs and s_bet[i+1,1] == hcs
        s_bet[i] = hcs 
      end
      i+=1
    end

    @strands_for_display_current = true
    @strands_for_display = [s_a, s_bet, s_b]
  end

=begin
  def vc_primary_add(c)
    @vc_primary << c
    @current = false
  end

  def vc_complement_add(c)
    @vc_complement << c
    @current = false
  end

  def hc_between_strands_add(c)
    @hc_between_strands << c
    @current = false
  end
=end
  #########
  protected
  #########

  def clean_all
    [@vc_primary, @vc_complement, @hc_between_strands].collect { |a| a.delete(nil); a.uniq!; a.sort! }
  end


end

end
end




More information about the bioruby-cvs mailing list