Class: Sequence

Inherits:
String
  • Object
show all
Defined in:
lib/parse_fasta/sequence.rb

Overview

Provide some methods for dealing with common tasks regarding nucleotide sequences.

Instance Method Summary collapse

Constructor Details

#initialize(str) ⇒ Sequence

Strips whitespace from the str argument before calling super

Examples:

Removes whitespace

Sequence.new "AA CC TT" #=> "AACCTT"

Raises:



40
41
42
43
44
45
46
 
# File 'lib/parse_fasta/sequence.rb', line 40

def initialize(str)
  if str.match(/>/)
    raise ParseFasta::SequenceFormatError
  end

  super(str.gsub(/ +/, ""))
end

Instance Method Details

#base_counts(count_ambiguous_bases = nil) ⇒ Hash

Returns a map of base counts

This method will check if the sequence is DNA or RNA and return a count map appropriate for each. If a truthy argument is given, the count of ambiguous bases will be returned as well.

If a sequence has both T and U present, will warn the user and keep going. Will return a map with counts of both, however.

Examples:

Get base counts of DNA sequence without ambiguous bases

Sequence.new('AcTGn').base_counts
#=> { a: 1, c: 1, t: 1, g: 1 }

Get base counts of DNA sequence with ambiguous bases

Sequence.new('AcTGn').base_counts(true)
#=> { a: 1, c: 1, t: 1, g: 1, n: 1 }

Get base counts of RNA sequence without ambiguous bases

Sequence.new('AcUGn').base_counts
#=> { a: 1, c: 1, u: 1, g: 1 }

Get base counts of DNA sequence with ambiguous bases

Sequence.new('AcUGn').base_counts(true)
#=> { a: 1, c: 1, u: 1, g: 1, n: 1 }

Returns:

  • (Hash) —

    A hash with base as key, count as value



101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
 
# File 'lib/parse_fasta/sequence.rb', line 101

def base_counts(count_ambiguous_bases=nil)
  s = self.downcase
  t = s.count('t')
  u = s.count('u')
  counts = { a: s.count('a'), c: s.count('c'), g: s.count('g') }

  if t > 0 && u == 0
    counts[:t] = t
  elsif t == 0 && u > 0
    counts[:u] = u
  elsif t > 0 && u > 0
    warn('ERROR: A sequence contains both T and U')
    counts[:t], counts[:u] = t, u
  end

  counts[:n] = s.count('n') if count_ambiguous_bases

  counts
end

#base_frequencies(count_ambiguous_bases = nil) ⇒ Hash

Returns a map of base frequencies

Counts bases with the base_counts method, then divides each count by the total bases counted to give frequency for each base. If a truthy argument is given, ambiguous bases will be included in the total and their frequency reported. Can discern between DNA and RNA.

If default or falsy argument is given, ambiguous bases will not be counted in the total base count and their frequency will not be given.

Examples:

Get base frequencies of DNA sequence without ambiguous bases

Sequence.new('AcTGn').base_counts
#=> { a: 0.25, c: 0.25, t: 0.25, g: 0.25 }

Get base counts of DNA sequence with ambiguous bases

Sequence.new('AcTGn').base_counts(true)
#=> { a: 0.2, c: 0.2, t: 0.2, g: 0.2, n: 0.2 }

Returns:

  • (Hash) —

    A hash with base as key, frequency as value



141
142
143
144
145
146
147
 
# File 'lib/parse_fasta/sequence.rb', line 141

def base_frequencies(count_ambiguous_bases=nil)
  base_counts = self.base_counts(count_ambiguous_bases)
  total_bases = base_counts.values.reduce(:+).to_f
  base_freqs =
    base_counts.map { |base, count| [base, count/total_bases] }.flatten
  Hash[*base_freqs]
end

#gc ⇒ 0, Float

Calculates GC content

Calculates GC content by dividing count of G + C divided by count of G + C + T + A + U. If there are both T's and U's in the Sequence, things will get weird, but then again, that wouldn't happen, now would it! Ambiguous bases are ignored similar to BioRuby.

Examples:

Get GC of a Sequence

Sequence.new('ACTg').gc #=> 0.5

Using with FastaFile#each_record

FastaFile.open('reads.fna', 'r').each_record do |header, sequence|
  puts [header, sequence.gc].join("\t")
end

Returns:

  • (0) —

    if the Sequence is empty or there are no A, C, T, G or U present

  • (Float) —

    if the GC content is defined for the Sequence



66
67
68
69
70
71
72
73
74
75
76
 
# File 'lib/parse_fasta/sequence.rb', line 66

def gc
  s = self.downcase
  c = s.count('c')
  g = s.count('g')
  t = s.count('t')
  a = s.count('a')
  u = s.count('u')

  return 0 if c + g + t + a + u == 0
  return (c + g) / (c + g + t + a + u).to_f
end

#rev_comp ⇒ Sequence

Note:

If Sequence contains non-IUPAC characters, these are not complemented

Returns a reverse complement of self

Examples:

Hanldes any IUPAC character and capitalization properly

Sequence.new("gARKbdctymvhu").rev_comp #=> "adbkraghvMYTc"

Leaves non IUPAC characters

Sequence.new("cccc--CCCcccga").rev_comp #=> "tcgggGGG--gggg""

Returns:

  • (Sequence) —

    a Sequence that is the reverse complement of self



162
163
164
165
166
167
168
169
170
171
172
173
 
# File 'lib/parse_fasta/sequence.rb', line 162

def rev_comp
  # if self.match(/T/i) && self.match(/U/i)
  #   raise Sequence::AmbiguousSequenceError
  # end

  # if self.match(/[^ATUGCYRSWKMBDHVN]/i)
  #   warn "WARNING: Sequence contains non IUPAC characters"
  # end

  self.reverse.tr("ATUGCYRSWKMBDHVNatugcyrswkmbdhvn",
                  "TAACGRYSWMKVHDBNtaacgryswmkvhdbn")
end