Class: Bio::DB::Tag::MD

Inherits:
Object
  • Object
show all
Includes:
Alignment::IteratePairs
Defined in:
lib/bio-sam-mutation/bio/db/tag/md.rb

Constant Summary collapse

@@regexp = 
/MD:Z:([\w^]+)/
@@format = 
/[\w^]+/
@@splitter = 
/(?<match>\d+)|(?<substitution>[GATCN]+)|\^(?<deletion>[GATCN]+)/
@@reference =

Operations that consume reference seqeunce:

/[msd]/

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(data) ⇒ MD

Returns a new instance of MD.



9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 9

def initialize(data)
		if data.is_a? String
			if data.match(@@regexp)
				@tag = $~[1]
			elsif data.match(@@format)
				#Assume tag given without MD:Z: leader
				@tag = data
			else
				raise "Tag not of expected format."
			end
		elsif data.is_a? Bio::DB::Tag
			@tag = data.value
			warn "Not an MD tag" if data.tag == "MD"
		else
			raise "Tag not of expected format."
		end

	# Splits the string into operations using the splitter regexp class variable, returns array of two-element arrays describing operations
	spl = @tag.scan(@@splitter)
	# Returns an array of matches [match,substition,deletion]
	# Although regexp captures are named, these don't get included automatically with scan as it doesn't return MatchData objects.
	spl.map! do |a|
		array = [["m", a[0]],["s", a[1]],["d", a[2]]]
		# Only one of these will be non-nil
		array.keep_if{|i| i[1]}
		array.map!{|i| if i[0] == "m" then i[1] = i[1].to_i end; i}
		array[0]
	end
	@pairs = spl

	@cumulative = []
	cumulative_length = 0
	read_length = 0
	@pairs.each do |q|
		p = q.dup
		case p[0]
			when "m"
				len = p[1]
				rlen = p[1]
			when "s"
				len = p[1].length
				rlen = p[1].length
			when "d"
				len = p[1].length
				# Deleted bases don't appear in the read, so don't count to the length
				rlen = 0
		end
		# third element in each array will be the total preceding length on the reference, i.e. the position of the operation.
		# fourth element is similar for the read.
		@cumulative << p.dup.push(cumulative_length).push(read_length)
		cumulative_length += len
		read_length += rlen
	end
end

Instance Attribute Details

#cumulative ⇒ Object

Returns the value of attribute cumulative.



3
4
5
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 3

def cumulative
  @cumulative
end

#pairs ⇒ Object

Returns the value of attribute pairs.



3
4
5
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 3

def pairs
  @pairs
end

#tag ⇒ Object

Returns the value of attribute tag.



3
4
5
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 3

def tag
  @tag
end

Instance Method Details

#deletions ⇒ Object 



64
65
66
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 64

def deletions
	report(/d/)
end

#reconstruct_tag(array = @pairs) ⇒ Object

Reconstruct a MD:Z tag from the pairs array



84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 84

def reconstruct_tag(array=@pairs)
	new_tag = []
	array.each do |p|
		case p[0]
			when "m"
				string = p[1].to_s
			when "s"
				string = p[1]
			when "d"
				string = "^"+p[1]
		end
		new_tag << string
	end
	new_tag.join("")
end

#ref_length ⇒ Object

Sums the total length of the reference sequence represented by the MD:Z tag (or part of)



103
104
105
106
107
108
109
110
111
112
113
114
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 103

def ref_length
	#Need the sum of all "movement" operations (i.e. numbers) as well as any substituted bases (count 1 each)
	if @tag =~ /^\d+$/
		@tag.to_i
	else
		temp_tag = @tag.dup
		temp_tag.gsub!(/\^/,"")  # Deletions need to be counted - sub the caret character out and count the remaining base characters
		movements = temp_tag.split(/[GATCN]+/).map(&:to_i).reduce(:+) # Sum numbers
		deletions = temp_tag.split(/\d+/).map(&:length).reduce(:+) # Sum number of base chars
		movements + deletions
	end
end

#report(regexp = /[sd]/) ⇒ Object

Report the positions of given events



73
74
75
76
77
78
79
80
81
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 73

def report(regexp=/[sd]/)
	to_return = []
	@cumulative.each do |p|
		if p[0] =~ regexp
			to_return << p
		end
	end
	to_return
end

#slice(offset, length) ⇒ Object

Given an offset in reference sequence and length, return an object corresponding to that subregion of the alignment



116
117
118
119
120
121
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 116

def slice(offset,length)
	new_array = iterate_pairs(@pairs,offset,length,@@reference)
	# Return a MDZ instance with just the new alignment
	new_tag = reconstruct_tag(new_array)
	Bio::DB::Tag::MD.new(new_tag)
end

#substitutions ⇒ Object 



68
69
70
 
# File 'lib/bio-sam-mutation/bio/db/tag/md.rb', line 68

def substitutions
	report(/s/)
end