Class: Bio::Alignment::CIGAR

Inherits:
Object
  • Object
show all
Includes:
IteratePairs
Defined in:
lib/bio-sam-mutation/bio/alignment/cigar.rb

Overview

Parse a CIGAR string An example from Exonerate output. Ideally will also allow SAM file input to be used.

1 : CGGCTATGGGGTCGTGGGTCCCGCGTTG-CTCTGGGGCTCGGCACCCTGGGGCGGCACGGCCGT :  63
    | | || | ||||||||||||||||||| |||||||||||||||||||||||||||||||||||
1 : CAG-TA-GTGGTCGTGGGTCCCGCGTTGTCTCTGGGGCTCGGCACCCTGGGGCGGCACGGCCGT :  62

ref: CAGTAGTGGTCGTGGGTCCCGCGTTGTCTCTGG… cigar: SP-A12_D02_2015-01-16.seq 0 611 + SP-A3_ref 0 621 + 2514 M 3 I 1 M 2 I 1 M 21 D 1 M 306 D 9 M 89 I 1 M 126 D 1 M 24 D 1 M 8 I 1 M 6 D 1 M 5 D 1 M 17 I are not counted in reference Regexp (from SAM specification) - but in exonerate the number comes first: ([0-9][MIDNSHP])|*

Class Attribute Summary collapse

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(string, ref = nil, source = "") ⇒ CIGAR

Returns a new instance of CIGAR.



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
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 27

def initialize(string,ref=nil,source="")
	# strip out whitespace
	string.gsub!(/\s+/,"")

	# Auto-detect source if not supplied
	if !(Bio::Alignment::CIGAR.regexps.keys.include? source)
	  Bio::Alignment::CIGAR.regexps.each do |k,v|
			# Look for match at start of string
			if m = string.match(v)
			  source = k if m.offset(0)[0] == 0
		  end
		end
		if source == ""
			raise "Source (e.g. 'exonerate', 'sam') not given and failed to auto-detect."
		end
   end
	# Make an array of pairs of of cigar elements:
	@pairs = string.scan(Bio::Alignment::CIGAR.regexps[source])
	if source == "exonerate"
		@pairs.map!{|pair| [pair[0].to_s, pair[1].to_i]}
	else
		# Provision to have number and identifier the other way round
		@pairs.map!{|pair| [pair[1].to_s, pair[0].to_i]}
	end

	# Include reference sequence if provided
	@reference = ref
	# Check length of reference = sum(M+D)?
	#warn "Reference length is not equal to that implied by CIGAR string: #{@reference.length}, #{self.reference_length}." unless @reference.length == self.reference_length

end

Class Attribute Details

.query_operations ⇒ Object

Returns the value of attribute query_operations.



16
17
18
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 16

def query_operations
  @query_operations
end

.reference_operations ⇒ Object

Returns the value of attribute reference_operations.



16
17
18
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 16

def reference_operations
  @reference_operations
end

.regexps ⇒ Object

Returns the value of attribute regexps.



16
17
18
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 16

def regexps
  @regexps
end

.subexp ⇒ Object

Returns the value of attribute subexp.



16
17
18
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 16

def subexp
  @subexp
end

Instance Attribute Details

#pairs ⇒ Object

Returns the value of attribute pairs.



25
26
27
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 25

def pairs
  @pairs
end

#reference ⇒ Object

Returns the value of attribute reference.



25
26
27
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 25

def reference
  @reference
end

Instance Method Details

#combine_adjacent ⇒ Object

TODO combine adjacent operations of the same type into a single pair



198
199
200
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 198

def combine_adjacent

end

#deleted_length ⇒ Object 



96
97
98
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 96

def deleted_length
	count_type("D")
end

#hgnc(reference_pos = 0, insertions = [], type = "g", *subs) ⇒ Object

Output hgnc variant format given reference position. Only deletions can be accurately annotated from the cigar string; insertions or wild type seqeunces return nil NB mutation calling and annotation now implemented as extension to Bio::DB::Alignment (SAM)



148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 148

def hgnc(reference_pos=0,insertions=[],type="g",*subs)
	if insertions
		if insertions.is_a? String
			insertions = [insertions]
		end
	end
	first_match = true
	total = 0
	hgnc_format = []
	@pairs.each do |pair|
		case pair[0]
			when "M"
				#break if first_match == false
				reference_pos += pair[1]
				total += pair[1]
				first_match = false
			when "D"
				deleted_bases = @reference[total,pair[1]].upcase
				if (pair[1] == 1)
					string = (reference_pos + 1).to_s
				else
					string = (reference_pos + 1).to_s + "_" + (reference_pos + pair[1]).to_s
				end
				string = string + "del" + deleted_bases
				hgnc_format << string
				total += pair[1]
			when "I"
				inserted_bases = (insertions.length == 0) ? "N" : insertions.shift

				hgnc_format << (reference_pos).to_s + "_" + (reference_pos + 1).to_s + "ins" + inserted_bases.upcase
		end
	end
	# Use for substitutions, but could also pass any other annotation to include in here, as an array of strings
	subs = subs.first # >1 arguments discarded
	if subs
		if (subs.length > 0 && (subs.is_a? Array))
			hgnc_format = hgnc_format + subs
		end
	end
	if hgnc_format.length == 0
		nil
	elsif hgnc_format.length == 1
		type.to_s + "." + hgnc_format[0]
	else
		type.to_s + "." + "[" + hgnc_format.join(";") + "]"
	end

end

#inserted_length ⇒ Object 



100
101
102
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 100

def inserted_length
	count_type("I")
end

#masked_length ⇒ Object 



104
105
106
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 104

def masked_length
	count_type(/[SH]/)
end

#matched_length ⇒ Object 



92
93
94
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 92

def matched_length
	count_type("M")
end

#positions(type) ⇒ Object

Returns a hash (keyed by operation type) of three element arrays: the start positions on the reference of operations of the given type(s) and the length of the operation, followed by query position (for e.g. retrieving inserted bases from SAM). A regexp can be used to specify multiple types e.g. /[ID]/.



204
205
206
207
208
209
210
211
212
213
214
215
216
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 204

def positions(type)
	total = 0
	qtotal = 0
	hash = Hash.new{|h,k| h[k] = []}
	@pairs.each do |pair|
		if pair[0].match(type)
			hash[$&] << [total, pair[1], qtotal]
		end
		total += pair[1] if pair[0].match Bio::Alignment::CIGAR.reference_operations
		qtotal += pair[1] if pair[0].match Bio::Alignment::CIGAR.query_operations
	end
	hash
end

#query(insertions = nil) ⇒ Object

Output a representation of the query: replace deleted portions with “-”, flag insertions with “*” or sim. Optionally provide the sequence (or symbols to use) of insertions, in order of appearence. Should be able to accept an array TODO: Add support for substitution highlighting (e.g lowercasing)



119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 119

def query(insertions=nil)
	if (insertions && (insertions.is_a? String))
		insertions = [insertions]
	end
	sequence = []
	total = 0
	@pairs.each do |pair|
		if pair[0].match("M")
			sequence << @reference[total..total+pair[1]-1].upcase
			total += pair[1]
		end
		if pair[0].match("I")
			if (insertions)
				insertion = insertions.shift.to_s
			else
				insertion = '['+pair[1].to_s+']'
			end
			sequence << insertion
		end
		if pair[0].match("D")
			pair[1].times{ sequence << "-" }
			total += pair[1]
		end
	end
	sequence.join("")
end

#query_length ⇒ Object 



112
113
114
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 112

def query_length
	count_type(Bio::Alignment::CIGAR.query_operations)
end

#reference_length ⇒ Object 



108
109
110
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 108

def reference_length
	count_type(Bio::Alignment::CIGAR.reference_operations)
end

#remove_empty! ⇒ Object 



87
88
89
90
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 87

def remove_empty!
	self.pairs.keep_if{|pair| pair[1] != 0 }
	self
end

#remove_small!(threshold = 1) ⇒ Object 



80
81
82
83
84
85
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 80

def remove_small!(threshold=1)
	# Deletions convert to matches, insertions just remove
	deletions_to_matches(threshold)
	remove_small_nonmatches(threshold)
	self
end

#subalignment(offset, length, regexp = Bio::Alignment::CIGAR.reference_operations) ⇒ Object Also known as: slice

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



60
61
62
63
64
65
66
67
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 60

def subalignment(offset,length,regexp=Bio::Alignment::CIGAR.reference_operations)
	new_array = iterate_pairs(@pairs,offset,length,regexp)
	# Return a CIGAR instance with just the new alignment
	new_string = new_array.join(" ")
	# -1 from offset as ruby string starts at zero
	new_cigar = Bio::Alignment::CIGAR.new(new_string,@reference[offset-1,length])
	new_cigar.remove_empty!
end

#subcigar(offset, length) ⇒ Object

Given a CIGAR-based [not reference - use subalignment] offset and length, return a subregion



71
72
73
74
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 71

def subcigar(offset,length)
	# No regexp - includes everything
	self.subalignment(offset,length,//)
end

#unmasked(offset, length) ⇒ Object 



76
77
78
 
# File 'lib/bio-sam-mutation/bio/alignment/cigar.rb', line 76

def unmasked(offset,length)
	self.subalignment(offset,length,/[MDI]/)
end