code_ruby/lib/maasha/seq.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 require 'maasha/digest'
  26 require 'maasha/patternmatcher'
  27 #require 'maasha/patscan'
  28
  29 # Residue alphabets
  30 DNA     = %w[a t c g]
  31 RNA     = %w[a u c g]
  32 PROTEIN = %w[f l s y c w p h q r i m t n k v a d e g]
  33 INDELS  = %w[. - _ ~]
  34
  35 # Quality scores bases
  36 SCORE_PHRED    = 33
  37 SCORE_ILLUMINA = 64
  38
  39 # Error class for all exceptions to do with Seq.
  40 class SeqError < StandardError; end
  41
  42 class Seq
  43   #include Patscan
  44   include PatternMatcher
  45
  46   attr_accessor :seq_name, :seq, :type, :qual
  47
  48   # Method that generates all possible oligos of a specifed length and type.
  49   def self.generate_oligos(length, type)
  50     raise SeqError, "Cannot generate negative oligo length: #{length}" if length <= 0
  51
  52     case type.downcase
  53     when /dna/     then alph = DNA
  54     when /rna/     then alph = RNA
  55     when /protein/ then alph = PROTEIN
  56     else
  57       raise SeqError, "Unknown sequence type: #{type}"
  58     end
  59
  60     oligos = [""]
  61
  62     (1 .. length).each do
  63       list = []
  64
  65       oligos.each do |oligo|
  66         alph.each do |char|
  67           list << oligo + char
  68         end
  69       end
  70
  71       oligos = list
  72     end
  73
  74     oligos
  75   end
  76
  77   # Initialize a sequence object with the following arguments:
  78   # - seq_name: Name of the sequence.
  79   # - seq: The sequence.
  80   # - type: The sequence type - DNA, RNA, or protein
  81   # - qual: An Illumina type quality scores string.
  82   def initialize(seq_name = nil, seq = nil, type = nil, qual = nil)
  83     @seq_name = seq_name
  84     @seq      = seq
  85     @type     = type
  86     @qual     = qual
  87   end
  88
  89   # Method that guesses and returns the sequence type
  90   # by inspecting the first 100 residues.
  91   def type_guess
  92     raise SeqError, "Guess failed: sequence is nil" if self.seq.nil?
  93
  94     case self.seq[0 ... 100].downcase
  95     when /[flpqie]/ then return "protein"
  96     when /[u]/      then return "rna"
  97     else                 return "dna"
  98     end
  99   end
 100
 101   # Method that guesses and sets the sequence type
 102   # by inspecting the first 100 residues.
 103   def type_guess!
 104     self.type = self.type_guess
 105   end
 106
 107   # Returns the length of a sequence.
 108   def length
 109     self.seq.nil? ? 0 : self.seq.length
 110   end
 111
 112   alias :len :length
 113
 114   # Return the number indels in a sequence.
 115   def indels
 116     regex = Regexp.new(/[#{Regexp.escape(INDELS.join(""))}]/)
 117     self.seq.scan(regex).size
 118   end
 119
 120   # Method that returns true is a given sequence type is DNA.
 121   def is_dna?
 122     self.type == 'dna'
 123   end
 124
 125   # Method that returns true is a given sequence type is RNA.
 126   def is_rna?
 127     self.type == 'rna'
 128   end
 129
 130   # Method that returns true is a given sequence type is protein.
 131   def is_protein?
 132     self.type == 'protein'
 133   end
 134
 135   # Method to transcribe DNA to RNA.
 136   def to_rna
 137     raise SeqError, "Cannot transcribe 0 length sequence" if self.length == 0
 138     raise SeqError, "Cannot transcribe sequence type: #{self.type}" unless self.is_dna?
 139     self.type = 'rna'
 140     self.seq.tr!('Tt','Uu')
 141   end
 142
 143   # Method to reverse-transcribe RNA to DNA.
 144   def to_dna
 145     raise SeqError, "Cannot reverse-transcribe 0 length sequence" if self.length == 0
 146     raise SeqError, "Cannot reverse-transcribe sequence type: #{self.type}" unless self.is_rna?
 147
 148     self.type = 'dna'
 149     self.seq.tr!('Uu','Tt')
 150   end
 151
 152   # Method that given a Seq entry returns a Biopieces record (a hash).
 153   def to_bp
 154     raise SeqError, "Missing seq_name" if self.seq_name.nil?
 155     raise SeqError, "Missing seq"      if self.seq.nil?
 156
 157     record             = {}
 158     record[:SEQ_NAME] = self.seq_name
 159     record[:SEQ]      = self.seq
 160     record[:SEQ_LEN]  = self.length
 161     record[:SCORES]   = self.qual if self.qual
 162     record
 163   end
 164
 165   # Method that given a Seq entry returns a FASTA entry (a string).
 166   def to_fasta(wrap = nil)
 167     raise SeqError, "Missing seq_name" if self.seq_name.nil?
 168     raise SeqError, "Missing seq"      if self.seq.nil?
 169
 170     seq_name = self.seq_name
 171     seq      = self.seq
 172
 173     unless wrap.nil?
 174       seq.gsub!(/(.{#{wrap}})/) do |match|
 175         match << "\n"
 176       end
 177
 178       seq.chomp!
 179     end
 180
 181     ">#{seq_name}\n#{seq}\n"
 182   end
 183
 184   # Method that generates a unique key for a
 185   # DNA sequence and return this key as a Fixnum.
 186   def to_key
 187     key = 0
 188
 189     self.seq.upcase.each_char do |char|
 190       key <<= 2
 191
 192       case char
 193       when 'A' then key |= 0
 194       when 'C' then key |= 1
 195       when 'G' then key |= 2
 196       when 'T' then key |= 3
 197       else raise SeqError, "Bad residue: #{char}"
 198       end
 199     end
 200
 201     key
 202   end
 203
 204   # Method to reverse complement sequence.
 205   def reverse_complement
 206     self.reverse
 207     self.complement
 208   end
 209
 210   alias :revcomp :reverse_complement
 211
 212   # Method to reverse the sequence.
 213   def reverse
 214     self.seq.reverse!
 215   end
 216
 217   # Method that complements sequence including ambiguity codes.
 218   def complement
 219     raise SeqError, "Cannot complement 0 length sequence" if self.length == 0
 220
 221     if self.is_dna?
 222       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn' )
 223     elsif self.is_rna?
 224       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn' )
 225     else
 226       raise SeqError, "Cannot complement sequence type: #{self.type}"
 227     end
 228   end
 229
 230   # Method that generates a random sequence of a given length and type.
 231   def generate(length, type)
 232     raise SeqError, "Cannot generate sequence length < 1: #{length}" if length <= 0
 233
 234     case type.downcase
 235     when "dna"
 236       alph = DNA
 237     when "rna"
 238       alph = RNA
 239     when "protein"
 240       alph = PROTEIN
 241     else
 242       raise SeqError, "Unknown sequence type: #{type}"
 243     end
 244
 245     seq_new   = Array.new(length) { alph[rand(alph.size)] }.join("")
 246     self.seq  = seq_new
 247     self.type = type.downcase
 248     seq_new
 249   end
 250
 251   # Method that returns a subsequence of from a given start position
 252   # and of a given length.
 253   def subseq(start, length = self.length - start)
 254     raise SeqError, "subsequence start: #{start} < 0"                                                if start  < 0
 255     raise SeqError, "subsequence length: #{length} < 1"                                              if length <= 0
 256     raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length
 257
 258     stop = start + length - 1
 259
 260     seq  = self.seq[start .. stop]
 261     qual = self.qual[start .. stop] unless self.qual.nil?
 262
 263     Seq.new(self.seq_name, seq, self.type, qual)
 264   end
 265
 266   # Method that replaces a sequence with a subsequence from a given start position
 267   # and of a given length.
 268   def subseq!(start, length = self.length - start)
 269     raise SeqError, "subsequence start: #{start} < 0"                                                if start  < 0
 270     raise SeqError, "subsequence length: #{length} < 1"                                              if length <= 0
 271     raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length
 272
 273     stop = start + length - 1
 274
 275     self.seq  = self.seq[start .. stop]
 276     self.qual = self.qual[start .. stop] unless self.qual.nil?
 277   end
 278
 279   # Method that returns a subsequence of a given length
 280   # beginning at a random position.
 281   def subseq_rand(length)
 282     if self.length - length + 1 == 0
 283       start = 0
 284     else
 285       start = rand(self.length - length + 1)
 286     end
 287
 288     self.subseq(start, length)
 289   end
 290
 291   # Method that returns the residue compositions of a sequence in
 292   # a hash where the key is the residue and the value is the residue
 293   # count.
 294   def composition
 295     comp = Hash.new(0);
 296
 297     self.seq.upcase.each_char do |char|
 298       comp[char] += 1
 299     end
 300
 301     comp
 302   end
 303
 304   # Method that returns the length of the longest homopolymeric stretch
 305   # found in a sequence.
 306   def homopol_max(min = 1)
 307     return 0 if self.seq.nil? or self.seq.empty?
 308
 309     found = false
 310
 311     self.seq.upcase.scan(/A{#{min},}|T{#{min},}|G{#{min},}|C{#{min},}|N{#{min},}/) do |match|
 312       found = true
 313       min   = match.size > min ? match.size : min
 314     end
 315
 316     return 0 unless found
 317
 318     min
 319   end
 320
 321   # Method that returns the percentage of hard masked residues
 322   # or N's in a sequence.
 323   def hard_mask
 324     ((self.seq.upcase.scan("N").size.to_f / (self.len - self.indels).to_f) * 100).round(2)
 325   end
 326
 327   # Method that returns the percentage of soft masked residues
 328   # or lower cased residues in a sequence.
 329   def soft_mask
 330     ((self.seq.scan(/[a-z]/).size.to_f / (self.len - self.indels).to_f) * 100).round(2)
 331   end
 332
 333   # Method to convert the quality scores from a specified base
 334   # to another base.
 335   def convert_phred2illumina!
 336     self.qual.gsub!(/./) do |score|
 337       score_phred  = score.ord - SCORE_PHRED
 338       raise SeqError, "Bad Phred score: #{score} (#{score_phred})" unless (0 .. 40).include? score_phred
 339       score_illumina = score_phred + SCORE_ILLUMINA
 340       score          = score_illumina.chr
 341     end
 342   end
 343
 344   # Method to convert the quality scores from Solexa odd/ratio to
 345   # Illumina format.
 346   def convert_solexa2illumina!
 347     self.qual.gsub!(/./) do |score|
 348       score = solexa_char2illumina_char(score)
 349     end
 350   end
 351
 352   private
 353
 354   # Method to convert a Solexa score (odd ratio) to
 355   # a phred (probability) integer score.
 356   def solexa2phred(score)
 357     (10.0 * Math.log(10.0 ** (score / 10.0) + 1.0, 10)).to_i
 358   end
 359
 360   # Method to convert a Solexa score encoded using base
 361   # 64 ASCII to a Phred score encoded using base 64 ASCII.
 362   def solexa_char2illumina_char(char)
 363     score_solexa = char.ord - 64
 364     score_phred  = solexa2phred(score_solexa)
 365     (score_phred + 64).chr
 366   end
 367 end
 368
 369 __END__