Class: GeneValidator::AlignmentValidation
- Inherits:
-
ValidationTest
- Object
- ValidationTest
- GeneValidator::AlignmentValidation
- Extended by:
- Forwardable
- Defined in:
- lib/genevalidator/validation_alignment.rb
Overview
This class contains the methods necessary for validations based on multiple alignment
Instance Attribute Summary collapse
-
#index_file_name ⇒ Object
readonly
Returns the value of attribute index_file_name.
-
#multiple_alignment ⇒ Object
readonly
Returns the value of attribute multiple_alignment.
-
#raw_seq_file ⇒ Object
readonly
Returns the value of attribute raw_seq_file.
-
#raw_seq_file_load ⇒ Object
readonly
Returns the value of attribute raw_seq_file_load.
Attributes inherited from ValidationTest
#cli_name, #description, #header, #hits, #prediction, #run_time, #short_header, #type, #validation_report
Instance Method Summary collapse
-
#array_to_ranges(ar) ⇒ Object
converts an array of integers into array of ranges.
-
#consensus_validation(prediction_raw, consensus) ⇒ Object
Returns the percentage of consesnsus residues from the ma that are in the prediction Params:
prediction_raw
:String
corresponding to the prediction sequenceconsensus
:String
corresponding to the statistical model Output:Fixnum
with the score. -
#extra_sequence_validation(prediction_raw, sm) ⇒ Object
Returns the percentage of extra sequences in the prediction with respect to the statistical model Params:
prediction
:String
corresponding to the prediction sequencesm
:String
corresponding to the statistical model Output:Fixnum
with the score. -
#gap_validation(prediction_raw, sm) ⇒ Object
Returns the percentage of gaps in the prediction with respect to the statistical model Params:
prediction
:String
corresponding to the prediction sequencesm
:String
corresponding to the statistical model Output:Fixnum
with the score. -
#get_consensus(ma = @multiple_alignment) ⇒ Object
Returns the consensus regions among a set of multiple aligned sequences i.e positions where there is the same element in all sequences Params:
ma
: array of Strings, corresponding to the multiple aligned sequences Output:String
with the consensus regions. -
#get_sm_pssm(ma = @multiple_alignment, threshold = 0.7) ⇒ Object
Builds a statistical model from a set of multiple aligned sequences based on PSSM (Position Specific Matrix) Params:
ma
: array of Strings, corresponding to the multiple aligned sequencesthreshold
: percentage of genes that are considered in statistical model Output:String
representing the statistical modelArray
with the maximum frequeny of the majoritary residue for each position. -
#initialize(prediction, hits) ⇒ AlignmentValidation
constructor
Initilizes the object Params:
prediction
: aSequence
object representing the blast queryhits
: a vector ofSequence
objects (representing blast hits)plot_path
: name of the fasta file. -
#isalpha(str) ⇒ Object
Returns true if the string contains only letters and false otherwise.
-
#multiple_align_mafft(prediction = @prediction, hits = @hits) ⇒ Object
Builds the multiple alignment between all the hits and the prediction using MAFFT tool Also creates a fasta file with the alignment Params:
prediction
: aSequence
object representing the blast queryhits
: a vector ofSequience
objects (usually representing blast hits)path
: path of mafft installation Output: Array of Strings, corresponding to the multiple aligned sequences the prediction is the last sequence in the vector. -
#plot_alignment(freq, ma = @multiple_alignment) ⇒ Object
Generates a json file cotaining data used for plotting lines for multiple hits alignment, prediction and statistical model Params:
freq
:String
residue frequency from the statistical modeloutput
: plot_path of the json filema
:String
array with the multiple alignmened hits and prediction. -
#remove_isolated_residues(seq, len = 2) ⇒ Object
Remove isolated residues inside long gaps from a given sequence Params:
seq
:String
: sequence of residueslen
:Fixnum
: number of isolated residues to be removed Output:String
: the new sequence. -
#run(n = 10) ⇒ Object
Find gaps/extra regions based on the multiple alignment of the first n hits Output:
AlignmentValidationOutput
object.
Constructor Details
#initialize(prediction, hits) ⇒ AlignmentValidation
Initilizes the object Params: prediction
: a Sequence
object representing the blast query hits
: a vector of Sequence
objects (representing blast hits) plot_path
: name of the fasta file
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 |
# File 'lib/genevalidator/validation_alignment.rb', line 102 def initialize(prediction, hits) super @short_header = 'MissingExtraSequences' @cli_name = 'align' @header = 'Missing/Extra Sequences' @description = 'Finds missing and extra sequences in the' \ ' prediction, based on the multiple alignment of' \ ' the best hits. Also counts the percentage of' \ ' the conserved regions that appear in the' \ ' prediction.' @raw_seq_file = opt[:raw_sequences] @index_file_name = config[:raw_seq_file_index] @raw_seq_file_load = config[:raw_seq_file_load] @db = opt[:db] @multiple_alignment = [] @num_threads = opt[:num_threads] @type = config[:type] end |
Instance Attribute Details
#index_file_name ⇒ Object (readonly)
Returns the value of attribute index_file_name.
93 94 95 |
# File 'lib/genevalidator/validation_alignment.rb', line 93 def index_file_name @index_file_name end |
#multiple_alignment ⇒ Object (readonly)
Returns the value of attribute multiple_alignment.
91 92 93 |
# File 'lib/genevalidator/validation_alignment.rb', line 91 def multiple_alignment @multiple_alignment end |
#raw_seq_file ⇒ Object (readonly)
Returns the value of attribute raw_seq_file.
92 93 94 |
# File 'lib/genevalidator/validation_alignment.rb', line 92 def raw_seq_file @raw_seq_file end |
#raw_seq_file_load ⇒ Object (readonly)
Returns the value of attribute raw_seq_file_load.
94 95 96 |
# File 'lib/genevalidator/validation_alignment.rb', line 94 def raw_seq_file_load @raw_seq_file_load end |
Instance Method Details
#array_to_ranges(ar) ⇒ Object
converts an array of integers into array of ranges
390 391 392 393 394 395 396 397 398 399 |
# File 'lib/genevalidator/validation_alignment.rb', line 390 def array_to_ranges(ar) prev = ar[0] ranges = ar.slice_before { |e| prev, prev2 = e, prev prev2 + 1 != e }.map { |a| a[0]..a[-1] } ranges end |
#consensus_validation(prediction_raw, consensus) ⇒ Object
Returns the percentage of consesnsus residues from the ma that are in the prediction Params: prediction_raw
: String
corresponding to the prediction sequence consensus
: String
corresponding to the statistical model Output: Fixnum
with the score
306 307 308 309 310 311 312 313 314 315 316 317 |
# File 'lib/genevalidator/validation_alignment.rb', line 306 def consensus_validation(prediction_raw, consensus) return 1 if prediction_raw.length != consensus.length # no of conserved residues among the hits no_conserved_residues = consensus.length - consensus.scan(/[\?-]/).length return 1 if no_conserved_residues == 0 # no of conserved residues from the hita that appear in the prediction no_conserved_pred = consensus.split(//).each_index.count { |j| consensus[j] != '-' && consensus[j] != '?' && consensus[j] == prediction_raw[j] } no_conserved_pred / (no_conserved_residues + 0.0) end |
#extra_sequence_validation(prediction_raw, sm) ⇒ Object
Returns the percentage of extra sequences in the prediction with respect to the statistical model Params: prediction
: String
corresponding to the prediction sequence sm
: String
corresponding to the statistical model Output: Fixnum
with the score
287 288 289 290 291 292 293 294 295 296 |
# File 'lib/genevalidator/validation_alignment.rb', line 287 def extra_sequence_validation(prediction_raw, sm) return 1 if prediction_raw.length != sm.length # find residues that are in the prediction # but not in the statistical model no_insertions = 0 (0..sm.length - 1).each do |i| no_insertions += 1 if prediction_raw[i] != '-' && sm[i] == '-' end no_insertions / (sm.length + 0.0) end |
#gap_validation(prediction_raw, sm) ⇒ Object
Returns the percentage of gaps in the prediction with respect to the statistical model Params: prediction
: String
corresponding to the prediction sequence sm
: String
corresponding to the statistical model Output: Fixnum
with the score
268 269 270 271 272 273 274 275 276 277 |
# File 'lib/genevalidator/validation_alignment.rb', line 268 def gap_validation(prediction_raw, sm) return 1 if prediction_raw.length != sm.length # find gaps in the prediction and # not in the statistical model no_gaps = 0 (0..sm.length - 1).each do |i| no_gaps += 1 if prediction_raw[i] == '-' && sm[i] != '-' end no_gaps / (sm.length + 0.0) end |
#get_consensus(ma = @multiple_alignment) ⇒ Object
Returns the consensus regions among a set of multiple aligned sequences i.e positions where there is the same element in all sequences Params: ma
: array of Strings, corresponding to the multiple aligned sequences Output: String
with the consensus regions
255 256 257 258 |
# File 'lib/genevalidator/validation_alignment.rb', line 255 def get_consensus(ma = @multiple_alignment) align = Bio::Alignment.new(ma) align.consensus end |
#get_sm_pssm(ma = @multiple_alignment, threshold = 0.7) ⇒ Object
Builds a statistical model from a set of multiple aligned sequences based on PSSM (Position Specific Matrix) Params: ma
: array of Strings, corresponding to the multiple aligned sequences threshold
: percentage of genes that are considered in statistical model Output: String
representing the statistical model Array
with the maximum frequeny of the majoritary residue for each position
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 |
# File 'lib/genevalidator/validation_alignment.rb', line 330 def get_sm_pssm(ma = @multiple_alignment, threshold = 0.7) sm = '' freq = [] (0..ma[0].length - 1).each do |i| freqs = Hash.new(0) ma.map { |seq| seq[i] }.each { |res| freqs[res] += 1 } # get the residue with the highest frequency max_freq = freqs.map { |_res, n| n }.max residue = (freqs.map { |res, n| n == max_freq ? res : [] }.flatten)[0] if residue == '-' freq.push(0) else freq.push(max_freq / (ma.length + 0.0)) end if max_freq / (ma.length + 0.0) >= threshold sm << residue else sm << '?' end end [sm, freq] end |
#isalpha(str) ⇒ Object
Returns true if the string contains only letters and false otherwise
384 385 386 |
# File 'lib/genevalidator/validation_alignment.rb', line 384 def isalpha(str) !str.match(/[^A-Za-z]/) end |
#multiple_align_mafft(prediction = @prediction, hits = @hits) ⇒ Object
Builds the multiple alignment between all the hits and the prediction using MAFFT tool Also creates a fasta file with the alignment Params: prediction
: a Sequence
object representing the blast query hits
: a vector of Sequience
objects (usually representing blast hits) path
: path of mafft installation Output: Array of Strings, corresponding to the multiple aligned sequences the prediction is the last sequence in the vector
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 |
# File 'lib/genevalidator/validation_alignment.rb', line 224 def multiple_align_mafft(prediction = @prediction, hits = @hits) fail unless prediction.is_a?(Query) && hits[0].is_a?(Query) = ['--maxiterate', '1000', '--localpair', '--anysymbol', '--quiet', '--thread', "#{@num_threads}"] mafft = Bio::MAFFT.new('mafft', ) sequences = hits.map(&:raw_sequence) sequences.push(prediction.protein_translation) report = mafft.query_align(sequences) # Accesses the actual alignment. align = report.alignment align.each_with_index do |s, _i| @multiple_alignment.push(s.to_s) end @multiple_alignment rescue raise NoMafftInstallationError end |
#plot_alignment(freq, ma = @multiple_alignment) ⇒ Object
Generates a json file cotaining data used for plotting lines for multiple hits alignment, prediction and statistical model Params: freq
: String
residue frequency from the statistical model output
: plot_path of the json file ma
: String
array with the multiple alignmened hits and prediction
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 |
# File 'lib/genevalidator/validation_alignment.rb', line 407 def plot_alignment(freq, ma = @multiple_alignment) # get indeces of consensus in the multiple alignment consensus = get_consensus(@multiple_alignment[0..@multiple_alignment.length - 2]) consensus_idxs = consensus.split(//).each_index.select { |j| isalpha(consensus[j]) } consensus_ranges = array_to_ranges(consensus_idxs) consensus_all = get_consensus(@multiple_alignment) consensus_all_idxs = consensus_all.split(//).each_index.select { |j| isalpha(consensus_all[j]) } consensus_all_ranges = array_to_ranges(consensus_all_idxs) match_alignment = ma[0..ma.length - 2].each_with_index.map { |seq, _j| seq.split(//).each_index.select { |j| isalpha(seq[j]) } } match_alignment_ranges = [] match_alignment.each { |arr| match_alignment_ranges << array_to_ranges(arr) } query_alignment = ma[ma.length - 1].split(//).each_index.select { |j| isalpha(ma[ma.length - 1][j]) } query_alignment_ranges = array_to_ranges(query_alignment) len = ma[0].length # plot statistical model data = freq.each_with_index.map { |h, j| { 'y' => ma.length, 'start' => j, 'stop' => j + 1, 'color' => 'orange', 'height' => h } } + # hits match_alignment_ranges.each_with_index.map { |ranges, j| ranges.map { |range| { 'y' => ma.length - j - 1, 'start' => range.first, 'stop' => range.last, 'color' => 'red', 'height' => -1 } } }.flatten + ma[0..ma.length - 2].each_with_index.map { |_seq, j| consensus_ranges.map { |range| { 'y' => j + 1, 'start' => range.first, 'stop' => range.last, 'color' => 'yellow', 'height' => -1 } } }.flatten + # plot prediction [{ 'y' => 0, 'start' => 0, 'stop' => len, 'color' => 'gray', 'height' => -1 }] + query_alignment_ranges.map { |range| { 'y' => 0, 'start' => range.first, 'stop' => range.last, 'color' => 'red', 'height' => -1 } }.flatten + # plot consensus consensus_all_ranges.map { |range| { 'y' => 0, 'start' => range.first, 'stop' => range.last, 'color' => 'yellow', 'height' => -1 } }.flatten y_axis_values = 'Prediction' (1..ma.length - 1).each { |i| y_axis_values << ", hit #{i}" } y_axis_values << ', Statistical Model' Plot.new(data, :align, 'Missing/Extra sequences Validation: Multiple Align. &' \ 'Statistical model of hits', 'Conserved Region, Yellow', 'Offset in the Alignment', '', ma.length + 1, y_axis_values) end |
#remove_isolated_residues(seq, len = 2) ⇒ Object
Remove isolated residues inside long gaps from a given sequence Params: seq
:String
: sequence of residues len
:Fixnum
: number of isolated residues to be removed Output: String
: the new sequence
363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 |
# File 'lib/genevalidator/validation_alignment.rb', line 363 def remove_isolated_residues(seq, len = 2) gap_starts = seq.to_enum(:scan, /(-\w{1,#{len}}-)/i).map { |_m| $`.size + 1 } # remove isolated residues gap_starts.each do |i| (i..i + len - 1).each do |j| seq[j] = '-' if isalpha(seq[j]) end end # remove isolated gaps res_starts = seq.to_enum(:scan, /([?\w]-{1,2}[?\w])/i).map { |_m| $`.size + 1 } res_starts.each do |i| (i..i + len - 1).each do |j| seq[j] = '?' if seq[j] == '-' end end seq end |
#run(n = 10) ⇒ Object
Find gaps/extra regions based on the multiple alignment of the first n hits Output: AlignmentValidationOutput
object
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 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 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 |
# File 'lib/genevalidator/validation_alignment.rb', line 126 def run(n = 10) n = 50 if n > 50 fail NotEnoughHitsError unless hits.length >= n fail unless prediction.is_a?(Query) && hits[0].is_a?(Query) start = Time.new # get the first n hits less_hits = @hits[0..[n - 1, @hits.length].min] useless_hits = [] # get raw sequences for less_hits less_hits.map do |hit| next unless hit.raw_sequence.nil? hit.raw_sequence = FetchRawSequences.run(hit.identifier, hit.accession_no) useless_hits.push(hit) if hit.raw_sequence.nil? end useless_hits.each { |hit| less_hits.delete(hit) } fail NoInternetError if less_hits.length == 0 # in case of nucleotide prediction sequence translate into protein # translate with the reading frame of all hits considered for alignment reading_frames = less_hits.map(&:reading_frame).uniq fail ReadingFrameError if reading_frames.length != 1 if @type == :nucleotide s = Bio::Sequence::NA.new(prediction.raw_sequence) prediction.protein_translation = s.translate(reading_frames[0]) end # multiple align sequences from less_hits with the prediction # the prediction is the last sequence in the vector multiple_align_mafft(prediction, less_hits) out = get_sm_pssm(@multiple_alignment[0..@multiple_alignment.length - 2]) sm = out[0] freq = out[1] # remove isolated residues from the predicted sequence index = @multiple_alignment.length - 1 prediction_raw = remove_isolated_residues(@multiple_alignment[index]) # remove isolated residues from the statistical model sm = remove_isolated_residues(sm) a1 = get_consensus(@multiple_alignment[0..@multiple_alignment.length - 2]) plot1 = plot_alignment(freq) gaps = gap_validation(prediction_raw, sm) extra_seq = extra_sequence_validation(prediction_raw, sm) consensus = consensus_validation(prediction_raw, a1) @validation_report = AlignmentValidationOutput.new(@short_header, @header, @description, gaps, extra_seq, consensus) @validation_report.plot_files.push(plot1) @validation_report.run_time = Time.now - start @validation_report rescue NotEnoughHitsError @validation_report = ValidationReport.new('Not enough evidence', :warning, @short_header, @header, @description) rescue NoMafftInstallationError @validation_report = ValidationReport.new('Mafft error', :error, @short_header, @header, @description) @validation_report.errors.push NoMafftInstallationError rescue NoInternetError @validation_report = ValidationReport.new('Internet error', :error, @short_header, @header, @description) @validation_report.errors.push NoInternetError rescue ReadingFrameError @validation_report = ValidationReport.new('Multiple reading frames', :error, @short_header, @header, @description) @validation_report.errors.push 'Multiple reading frames Error' rescue @validation_report = ValidationReport.new('Unexpected error', :error, @short_header, @header, @description) @validation_report.errors.push 'Unexpected Error' end |