Module: FlAnalysis

Defined in:
lib/full_lengther_next/fl_analysis.rb

Instance Method Summary collapse

Instance Method Details

#analiza_orf_y_fl(seq, hit, options, db_name) ⇒ Object 



8
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
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 8

def analiza_orf_y_fl(seq, hit, options, db_name)
  query_fasta = seq.seq_fasta.upcase.dup # Upcase for prevents complications with masked sequences, dup for discard changes
  if hit.count > 1 # if the sequence has more than one hit, the frames are checked and fixed to get a single hit
      seq_unida = UneLosHit.new(hit, query_fasta)
      full_prot =    seq_unida.full_prot 
      query_fasta =  seq_unida.output_seq # repaired fasta
      final_hit =    seq_unida.final_hit   # single hit
      $global_warnings +=  seq_unida.msgs   # warning messages
  else
    query_fasta = reverse_seq(query_fasta, hit.first) if hit.first.q_frame < 0 # si la secuencia esta al reves le damos la vuelta
    final_hit = hit.first # single hit
  end
  query_fasta = exonerate_fix_frame_shift(query_fasta, hit) if options[:exonerate] 

  full_prot = query_fasta[final_hit.q_frame-1, query_fasta.length+1].translate
  original_query_coordinates = [final_hit.q_beg, final_hit.q_end] ## VERBOSE
  seq.show_alignment(final_hit, query_fasta, show_nts) if  $verbose > 2 ## VERBOSE
  atg_status, tmp_prot = set_start_codon(final_hit, options[:distance], full_prot, query_fasta)
  end_status, final_prot = find_end(final_hit, options[:distance], tmp_prot, query_fasta)

  puts "\n------------------- POST EXTENSION---------------------" if $verbose > 1 ## VERBOSE
  seq.show_alignment(final_hit, query_fasta, show_nts, original_query_coordinates) if  $verbose > 1 ## VERBOSE
  puts "ATG: #{atg_status}  STOP: #{end_status}" if  $verbose > 2 ## VERBOSE

  # decide the sequence status (Complete, Putative Complete, Internal, N-terminus, Putative N-terminus, C-terminus)
  type, status = determine_status(atg_status, end_status)
  status = compare_seq_length_with_subject(final_prot, options[:distance], final_hit, type, status)
  if final_prot.length >= 25 && final_prot.length.to_f/final_hit.full_subject_length >= options[:subject_coverage] # Prot length min of 25 aa and subject coverage by generated prot of 25%
    save_annotations(seq, final_hit, type, status, final_prot, query_fasta, db_name)
  end
end

#compare_seq_length_with_subject(final_prot, distance, final_hit, type, status) ⇒ Object 



183
184
185
186
187
188
189
190
191
192
193
194
195
196
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 183

def compare_seq_length_with_subject(final_prot, distance, final_hit, type, status)
  if final_prot.length - 2 * distance > final_hit.s_len
    $global_warnings << ['SeqLonger', final_prot.length, final_hit.s_len]
  elsif final_prot.length + 2 * distance < final_hit.s_len
    $global_warnings << ['SeqShorter', final_prot.length, final_hit.s_len]
    if final_prot.length + 100 < final_hit.s_len || final_prot.length*2 < final_hit.s_len       
      if type == COMPLETE
        status = false
        $global_warnings << 'VeryShorter'
      end
    end
  end
  return status
end

#determine_status(atg_status, end_status) ⇒ Object 



162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 162

def determine_status(atg_status, end_status)
  if atg_status != 'incomplete' && end_status != 'incomplete' # proteina completa
    type = COMPLETE
  elsif atg_status == 'incomplete' && end_status == 'incomplete' # region intermedia
    type = INTERNAL
  elsif atg_status != 'incomplete' && end_status == 'incomplete' # tenemos el principio de la proteina
    type = N_TERMINAL
  elsif atg_status == 'incomplete' && end_status != 'incomplete' # tenemos el final de la proteina
    type = C_TERMINAL
  end
  
  if atg_status == 'putative' || end_status == 'putative'
    status = false # Putative
  else
    status = true # Sure
  end

  return type, status
end

#exonerate_fix_frame_shift(query_fasta, hit) ⇒ Object 



239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 239

def exonerate_fix_frame_shift(query_fasta, hit)
  frame_shifts = []
  added_nts = 0
  hit.each_with_index do |hsp, num|
    if hsp.class.to_s == 'ExoBlastHit' #Only this type of class of BlastHit has frameshift attributes
      if !hsp.q_frameshift.empty? #There is frameshift
        hsp.q_frameshift.each do |position, num_nts|
          local_add = 3 - num_nts
          fs_final_position = position + num_nts 
          $global_warnings << ['ExFrameS', fs_final_position]
          frame_shifts << [fs_final_position, local_add]
          added_nts += local_add
        end
      end
    end
    hsp.q_beg += added_nts if num > 0
    hsp.q_end += added_nts
  end
  add = 0
  frame_shifts.each do |position, num_nts|
    query_fasta = query_fasta.insert(position+add, 'n'*num_nts)
    add += num_nts
  end
  return query_fasta
end

#find_end(final_hit, max_distance, tmp_prot, query_fasta) ⇒ Object 



105
106
107
108
109
110
111
112
113
114
115
116
117
118
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 105

def find_end(final_hit, max_distance, tmp_prot, query_fasta)
  frame_shift = check_frame_shift(final_hit)
  beg_end_string =(final_hit.q_end-final_hit.q_beg)/3 - max_distance # Begin of terminal region (Coordinate) in tmp_prot
  atg_substring = tmp_prot[0..beg_end_string] # prot without terminal region
  end_substring = tmp_prot[beg_end_string + 1 ..tmp_prot.length-1] #Take 3' of unigen
  #puts "\e[32m\nfinal_hit.q_end-final_hit.q_beg: #{final_hit.q_end-final_hit.q_beg} /3  - max_distance: #{max_distance}\e[0m"
  #puts "\e[33mbeg_end_string: #{beg_end_string}\e[0m"
  #puts "\e[35mtmp_prot.length: #{tmp_prot.length}\e[0m"
  if beg_end_string < 0 || end_substring.nil? #Sequences whose homology is at end of it and dont't exits the 3' part of unigene
    atg_substring = tmp_prot
    end_substring = ''
    end_status = 'incomplete'
  else
    end_status = 'putative'
    putative_end = end_substring.index('*')
    end_substring = end_substring[0 .. putative_end] if putative_end
    
    s_end_resto = final_hit.s_len - (final_hit.s_end + 1) # en el subject, numero de aas que necesito cubrir
    q_end_resto = (query_fasta.length - final_hit.q_end)/3 # en el query, numero de aas que tengo 
    sq_end_distance = q_end_resto - s_end_resto # La diferencia se hace a partir del final del hit para que el calculo no quede sesgado en caso de que la secuecia este truncada por 5'
    
    if (final_hit.align_len == final_hit.s_len && putative_end)||(sq_end_distance.abs  <= max_distance && putative_end && putative_end <= max_distance*2) #Stop in a Full-length. max_distance *2 is set by de margin of +-15aa at the end of aligment 
      end_status = 'complete'
    elsif sq_end_distance  < max_distance # si no tenemos suficiente secuencia para tener el stop (nos faltan 15 aas o mas)
      end_status = 'incomplete'
      if putative_end
        $global_warnings << ['UnexpSTOP3pDist', sq_end_distance.abs]
      else
        $global_warnings << ['DistSubj', sq_end_distance.abs]
      end
    else # tenemos suficiente secuencia
      if putative_end # tenemos un stop
        #beg_end_string indica en que punto del unigen se encuentra el area de busqueda del codon stop
        stop_q_s = beg_end_string + putative_end - final_hit.s_len # Space between query's stop and subject's stop
        if stop_q_s.abs <= max_distance #Stop codon is in search region
          end_status = 'complete'
        elsif stop_q_s < 0
          $global_warnings << 'UnexpSTOP3p'
        elsif stop_q_s > 0
          end_status = 'complete'
          $global_warnings << 'QueryTooLong'
        end
      else # no tenemos codon de parada pero tenemos suficiente secuencia
        end_status = 'incomplete'
        $global_warnings << 'ProtFusion'
      end
    end
  end
  final_prot = atg_substring + end_substring
  end_status = 'complete' if final_prot.length == final_hit.s_len+1 && final_prot[final_prot.length-1] == '*'
  new_q_end = final_hit.q_beg-1 + final_prot.length * 3 + frame_shift
  modify_3p_align(new_q_end, final_hit, query_fasta, final_prot) if  $verbose > 1
  final_hit.q_end = new_q_end  
  return end_status, final_prot
end

#find_start(subject_start, amine_seq, distance) ⇒ Object 



68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 68

def find_start(subject_start, amine_seq, distance)    
  atg_status = 'putative' # complete, incomplete or putative
  stop_codon = amine_seq.rindex('*')
  if !stop_codon.nil? # tenemos un codon de parada en el amine_seq 5 prima
    _5prime_UTR = amine_seq.length - 10 - subject_start # marcamos la distancia al s_beg desde el principio del amine_seq
    amine_seq = amine_seq[stop_codon + 1 .. amine_seq.length - 1]
    first_m = amine_seq.index('M')
    if stop_codon <= _5prime_UTR # Ver si stop está en zona 5 prima UTR
      if first_m # tenemos M 
        amine_seq = amine_seq[first_m .. amine_seq.length - 1]          
        atg_status = 'complete'
      else # con STOP pero sin M 
        $global_warnings << 'noM1'
      end
    else # esta despues, un cambio de fase impide analizar el principio
      $global_warnings << 'UnexpSTOP5p'
        amine_seq = amine_seq[first_m .. amine_seq.length - 1] if first_m # tenemos M
    end
  else # no hay stop codon
    first_m = amine_seq.index('M')
    if first_m # tenemos M
      amine_seq = amine_seq[first_m .. amine_seq.length - 1]
      m_distance = (subject_start - amine_seq.length).abs - 10
      if m_distance.abs > distance*2 # con atg pero muy lejos del inicio que marca el subject
        $global_warnings << 'NoStopMfar'
        atg_status = 'incomplete'
      else # Tenemos M
        atg_status = 'complete'
      end
    else # sin M
      $global_warnings << 'noM2'
    end
  end
  return amine_seq, atg_status
end

#mark_nt_seqs(final_hit, query_fasta) ⇒ Object 



221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 221

def mark_nt_seqs(final_hit, query_fasta)
  atg = query_fasta[final_hit.q_beg..final_hit.q_beg + 2]
  mark_atg = nil
  if atg == 'ATG'
    mark_atg = '_-_'  
  end
  stop = query_fasta[final_hit.q_end - 2..final_hit.q_end]
  mark_stop = nil
  if stop == 'TAG' || stop == 'TGA' || stop == 'TAA'
    mark_stop = '___'
  end
  seq5p = query_fasta[0..final_hit.q_beg-1]
  orf = query_fasta[final_hit.q_beg..final_hit.q_end]
  seq3p = query_fasta[final_hit.q_end..query_fasta.length]
  nt_seq = "#{seq5p}#{mark_atg}#{orf}#{mark_stop}#{seq3p}"
  return nt_seq
end

#modify_3p_align(new_q_end, final_hit, query_fasta, final_prot) ⇒ Object

For visual report



274
275
276
277
278
279
280
281
282
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 274

def modify_3p_align(new_q_end, final_hit, query_fasta, final_prot) ## For visual report
  if new_q_end > final_hit.q_end #There is an align extension
    extend_align = query_fasta[final_hit.q_end+1 .. new_q_end].translate
    final_hit.q_seq = final_hit.q_seq + extend_align
  elsif new_q_end < final_hit.q_end #The align is cutted
    upper_limit = final_prot.length - 1 + final_hit.q_seq.count('-')
    final_hit.q_seq = final_hit.q_seq[0 .. upper_limit]
  end
end

#modify_5p_align(new_q_beg, final_hit, query_fasta) ⇒ Object

For visual report



285
286
287
288
289
290
291
292
293
294
295
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 285

def modify_5p_align(new_q_beg, final_hit, query_fasta) ## For visual report
  if new_q_beg < final_hit.q_beg #There is an align extension
    extend_align = query_fasta[new_q_beg .. final_hit.q_beg-1].translate
    final_hit.q_seq = extend_align + final_hit.q_seq
  elsif new_q_beg > final_hit.q_beg #The align is cut
    seq_cut = (new_q_beg - final_hit.q_beg)/3
    gaps = final_hit.q_seq[0..seq_cut].count('-')
    seq_cut += gaps
    final_hit.q_seq = final_hit.q_seq[seq_cut .. final_hit.q_seq.length-1]
  end
end

#save_annotations(seq, final_hit, type, status, final_prot, query_fasta, db_name) ⇒ Object 



199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 199

def save_annotations(seq, final_hit, type, status, final_prot, query_fasta, db_name)
  # if the sequence is Complete or it hasn't previous info will be saved
  if seq.type == UNKNOWN || (type == COMPLETE && seq.type != COMPLETE)
    seq.type = type
    seq.status = status
    seq.db_name = db_name
    seq.seq_fasta = query_fasta
    seq.seq_aa = final_prot
    seq.hit = final_hit
    seq.warnings($global_warnings)
    $global_warnings = [] # Clean all warnings for current sequence
    seq.seq_nt = mark_nt_seqs(final_hit, query_fasta)
    if type == COMPLETE
      seq.ignore = true
    end
  end
  if  $verbose > 2
    puts "\e[1mStruct annot: #{seq.prot_annot_calification}\e[0m"
  end
end

#set_start_codon(final_hit, distance, full_prot, query_fasta) ⇒ Object 



41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 41

def set_start_codon(final_hit, distance, full_prot, query_fasta)
  q_index_start = contenidos_en_prot(final_hit.q_seq, full_prot) 
  atg_status = nil
  _5prima = q_index_start + distance

  if  final_hit.s_beg == 0 && final_hit.q_seq[0] == 'M' && final_hit.s_seq[0] == 'M' #there is M in query and subject at first position of alignment and subject's M is in first position
    atg_status = 'complete'
    tmp_prot = full_prot[q_index_start..full_prot.length]
  elsif _5prima >= final_hit.s_beg
    amine_seq = full_prot[0, _5prima] #Contiene parte amino de la proteina
    carboxile_seq = full_prot[_5prima, full_prot.length - _5prima] #Contiene parte carboxilo de la proteina hasta el fin de la secuencia
    length_before_cut = amine_seq.length
    amine_seq, atg_status = find_start(final_hit.s_beg, amine_seq, distance) # to look for the beginning of the protein
    tmp_prot = "#{amine_seq}#{carboxile_seq}" # merge seqs in prot
    new_q_beg = final_hit.q_frame-1 + (length_before_cut - amine_seq.length) * 3
    modify_5p_align(new_q_beg, final_hit, query_fasta)  if  $verbose > 1 ## VERBOSE, Modify query align
    final_hit.q_beg = new_q_beg # to get the value of the start_ORF index
  else
    $global_warnings << 'UnexpStopBegSeq' if full_prot[0, q_index_start].rindex('*')
    atg_status = 'incomplete'
    tmp_prot = full_prot
  end

  return atg_status, tmp_prot
end

#show_nts ⇒ Object

VERBOSE METHODS



267
268
269
270
271
 
# File 'lib/full_lengther_next/fl_analysis.rb', line 267

def show_nts
  show = false 
  show = true if $verbose && $verbose > 3
  return show
end