Module: ChimericSeqs
- Defined in:
- lib/full_lengther_next/classes/chimeric_seqs.rb
Constant Summary collapse
- BEG =
0
- STOP =
1
- HIT =
2
Instance Method Summary collapse
- #cluster_query_hits(query) ⇒ Object
- #confirm_chimeras(homology_zones, db_path, ident_thresold) ⇒ Object
- #define_hit_limits(hits) ⇒ Object
- #define_homology_zones(query, options, query_fasta) ⇒ Object
- #delete_zones(overlapping_zones, zones) ⇒ Object
- #do_clustal(seq_fasta) ⇒ Object
- #duplicate_hits(query) ⇒ Object
- #fragment_chimera(query, seq, hit, hit_position, hit_limits, num_homology_zones, options, db_name, cut_positions) ⇒ Object
- #get_hits(query, ref_hit) ⇒ Object
- #get_limits(hit) ⇒ Object
- #get_limits_s(hit) ⇒ Object
- #hit_is_in?(h_beg, h_end, hit) ⇒ Boolean
- #hit_move_limits(hit, q_add, s_add) ⇒ Object
- #hit_set_q_len(hit, q_len) ⇒ Object
- #min_distance_between_homology_zones(homology_zones) ⇒ Object
- #move_limits(hit, q_add, s_add) ⇒ Object
- #overlapping_zones(zones) ⇒ Object
- #search_chimeras(seq, blast_query, options, db_name, db_path) ⇒ Object
- #search_ident_prots(homology_zones, distances, ident_thresold) ⇒ Object
- #set_cut_positions(homology_zones) ⇒ Object
- #set_limits(hit, q_beg, q_end, s_beg, s_end) ⇒ Object
Instance Method Details
#cluster_query_hits(query) ⇒ Object
211 212 213 214 215 216 217 218 219 220 221 222 223 224 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 211 def cluster_query_hits(query) hits = [] acc_hit = [] query.hits.each do |hit| ind = acc_hit.index(hit.acc) if ind.nil? acc_hit << hit.acc hits << [hit] else hits[ind] << hit end end return hits end |
#confirm_chimeras(homology_zones, db_path, ident_thresold) ⇒ Object
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 63 def confirm_chimeras(homology_zones, db_path, ident_thresold) acc_hit = homology_zones.map{|zone| zone[HIT].first.acc} seq_fasta = %x[blastdbcmd -db #{db_path} -entry #{acc_hit.join(',')}] seq_fasta << ">remove\nALGO\n" #Needed for clustal-omega display the dist-matrix, requires unless 3 sequences to do it clustal_matrix = do_clustal(seq_fasta) clustal_matrix.shift #Remove header clustal_matrix.pop #Remove false sequence clustal_hits = [] distances = [] clustal_matrix.each do |line| fields = line.split fields.pop #Remove data belong to false sequence fields.shift #Remove prot name distances << fields.map! {|field| field.to_f} end delete_positions = search_ident_prots(homology_zones, distances, ident_thresold) delete_zones(delete_positions, homology_zones) end |
#define_hit_limits(hits) ⇒ Object
187 188 189 190 191 192 193 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 187 def define_hit_limits(hits) limits=[] hits.each do |hit| limits << get_limits(hit) end return limits end |
#define_homology_zones(query, options, query_fasta) ⇒ Object
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 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 150 def define_homology_zones(query, , query_fasta) # Define hit limits #--------------------- hits = cluster_query_hits(query) #Hsp packages hits_limits = define_hit_limits(hits) # Define homology zones #------------------------ #First homology zone zones = [[hits_limits.first[BEG], hits_limits.first[STOP], hits.first]] ref_hit_beg = hits_limits.first[BEG] ref_hit_end = hits_limits.first[STOP] #Other homology zone hits_limits.each_with_index do |hit, i| coincidences = 0 zones.each do |zone| if hit_is_in?(zone[BEG], zone[STOP], hit) # Extender zona de homologia si coinciden en zona zone[BEG] = [zone[BEG],hit[BEG]].min zone[STOP] = [zone[STOP],hit[STOP]].max coincidences+=1 end end if coincidences == 0 zones << [hit[BEG], hit[STOP], hits[i]] end end zones.sort!{|e1,e2| e1[BEG] <=> e2[BEG]} # Delete overlapping homology zones #------------------------------------ overlapping_zones = overlapping_zones(zones) delete_zones(overlapping_zones, zones) return zones end |
#delete_zones(overlapping_zones, zones) ⇒ Object
226 227 228 229 230 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 226 def delete_zones(overlapping_zones, zones) overlapping_zones.each do |zone| zones.delete_at(zone) end end |
#do_clustal(seq_fasta) ⇒ Object
323 324 325 326 327 328 329 330 331 332 333 334 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 323 def do_clustal(seq_fasta) cmd='clustalo -i - -o /dev/null --percent-id --full --distmat-out=/dev/stdout --force' clustal_matrix = nil IO.popen(cmd,'w+') {|clustal| clustal.sync = TRUE clustal.write(seq_fasta) clustal.close_write clustal_matrix = clustal.readlines clustal.close_read } return clustal_matrix end |
#duplicate_hits(query) ⇒ Object
279 280 281 282 283 284 285 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 279 def duplicate_hits(query) dup_hits=[] query.hits.each do |hit| dup_hits << hit.dup end return dup_hits end |
#fragment_chimera(query, seq, hit, hit_position, hit_limits, num_homology_zones, options, db_name, cut_positions) ⇒ Object
102 103 104 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 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 102 def fragment_chimera(query, seq, hit, hit_position, hit_limits, num_homology_zones, , db_name, cut_positions) # Prepare new seq and query #---------------------------- query_bak = query.dup query_bak.hits = hit # Here, hit is an array of hsps query_bak.query_def += "_split_#{hit_position}" seq_bak = seq.dup seq_bak.reset_classification seq_bak.clean_warnings seq_bak.seq_name += "_split_#{hit_position}" seq_bak.clean_orfs seq_bak.save_fasta = TRUE seq_bak.ignore = FALSE # Cut sequence and move hit/hsps limits #---------------------------------------- if hit_position == 0 #First zone limit = 0 if hit.first.q_frame < 0 #Hit reversed hit.first.q_frame = -1 end else #Middle & last zone limit = cut_positions[BEG]#hit_limits[BEG] hit_move_limits(hit, -limit, 0) #Redefine hit limits on new sequence after cut if hit.first.q_frame >= 0 hit.first.q_frame=1 elsif hit_position < num_homology_zones-1 #Last zone keeps his original frame because it's composed by the hit and the terminal sequence (Here hit is reversed). hit.first.q_frame=-1 end end if hit_position == num_homology_zones-1 #Last zone seq_bak.seq_fasta = seq.seq_fasta[cut_positions[BEG]..seq.fasta_length-1]#[hit_limits[BEG]..seq.fasta_length-1] else # Beginning & Middle zone seq_bak.seq_fasta = seq.seq_fasta[limit..cut_positions[STOP]]#[limit..hit_limits[STOP]] end seq_length = seq_bak.seq_fasta.length query_bak.full_query_length = seq_length seq_bak.fasta_length = seq_length hit_set_q_len(hit, seq_length) # Full length analisys of fragment #---------------------------------------- analiza_orf_y_fl(seq_bak, query_bak.hits, , db_name) return seq_bak end |
#get_hits(query, ref_hit) ⇒ Object
255 256 257 258 259 260 261 262 263 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 255 def get_hits(query, ref_hit) all_hits=[] query.hits.each do |hit| if hit.acc == ref_hit.acc all_hits << hit end end return all_hits end |
#get_limits(hit) ⇒ Object
195 196 197 198 199 200 201 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 195 def get_limits(hit) coordenates=[] hit.map{|h| coordenates << h.q_beg; coordenates << h.q_end} # BEG END limits=[coordenates.min, coordenates.max] return limits end |
#get_limits_s(hit) ⇒ Object
203 204 205 206 207 208 209 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 203 def get_limits_s(hit) coordenates=[] hit.map{|h| coordenates << h.s_beg; coordenates << h.s_end} # BEG END limits=[coordenates.min, coordenates.max] return limits end |
#hit_is_in?(h_beg, h_end, hit) ⇒ Boolean
246 247 248 249 250 251 252 253 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 246 def hit_is_in?(h_beg, h_end, hit) is=FALSE # CONTIENE #OVERLAP if h_beg <= hit[BEG] && h_end > hit[BEG] || hit[BEG] <= h_beg && hit[STOP] > h_beg is=TRUE end return is end |
#hit_move_limits(hit, q_add, s_add) ⇒ Object
306 307 308 309 310 311 312 313 314 315 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 306 def hit_move_limits(hit, q_add, s_add) if hit.class.to_s == 'Array' hit.each do |hsp| move_limits(hsp, q_add, s_add) end elsif hit.class.to_s == 'Hit' #puts "\e[35m#{hit.acc}\t#{hit.q_beg}\t#{hit.q_end}\t#{hit.s_beg}\t#{hit.s_end}\t#{hit.reversed}\e[0m" move_limits(hit, q_add, s_add) end end |
#hit_set_q_len(hit, q_len) ⇒ Object
317 318 319 320 321 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 317 def hit_set_q_len(hit, q_len) hit.each do |hsp| hsp.q_len=q_len end end |
#min_distance_between_homology_zones(homology_zones) ⇒ Object
266 267 268 269 270 271 272 273 274 275 276 277 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 266 def min_distance_between_homology_zones(homology_zones) distance=nil homology_zones.each_with_index do |zone,i| if i > 0 local_distance=homology_zones[i][BEG] - homology_zones[i-1][STOP] if distance.nil? || distance > local_distance distance=local_distance end end end return distance end |
#move_limits(hit, q_add, s_add) ⇒ Object
294 295 296 297 298 299 300 301 302 303 304 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 294 def move_limits(hit, q_add, s_add) hit.q_beg+=q_add hit.q_end+=q_add hit.s_beg+=s_add hit.s_end+=s_add if hit.class.to_s == 'ExoBlastHit' && !hit.q_frameshift.empty? #There is frameshift hit.q_frameshift.map!{|fs| [fs.first + q_add, fs.last] } end end |
#overlapping_zones(zones) ⇒ Object
232 233 234 235 236 237 238 239 240 241 242 243 244 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 232 def overlapping_zones(zones) delete_zones=[] zones.each_with_index do |zone, i| if i>0 if zone[BEG]< zones[i-1][STOP] delete_zones << i delete_zones << i-1 end end end delete_zones.uniq! return delete_zones end |
#search_chimeras(seq, blast_query, options, db_name, db_path) ⇒ Object
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 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 9 def search_chimeras(seq, blast_query, , db_name, db_path) # DETECTION #---------------------- homology_zones = [] cut_positions = [] if blast_query.hits.length > 1 homology_zones = define_homology_zones(blast_query, , seq.seq_fasta) cut_positions = set_cut_positions(homology_zones) if homology_zones.length > 1 end # CONFIRMATION #---------------------- num_homology_zones = homology_zones.length if num_homology_zones > 1 && [:chimera].include?('r') confirm_chimeras(homology_zones, db_path, [:ident_thresold]) # Check if prots are differents or not num_homology_zones = homology_zones.length end # SPLICING #-------------------- new_seqs=[] if num_homology_zones > 1 #In this case the sequence is a chimera seq.format_chimera! homology_zones.each_with_index do |hom_zone, i| seq.hit << hom_zone[HIT].first.dup #Save hit before modified it for write output purposes hit_limits = get_limits(hom_zone[HIT])# Take beginning and end of hit on query, hit can be composed by unsorted or antisense hsps if [:chimera].include?('c') && hit_limits[STOP]-hit_limits[BEG]> [:min_nucleotides] new_seqs << fragment_chimera(blast_query, seq, hom_zone[HIT], i, hit_limits, num_homology_zones, , db_name, cut_positions[i]) seq.warnings('SOLVED') end end else new_seqs = nil #Sequence isn't chimera end return new_seqs end |
#search_ident_prots(homology_zones, distances, ident_thresold) ⇒ Object
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 85 def search_ident_prots(homology_zones, distances, ident_thresold) delete_positions = [] n_homology_zones = homology_zones.length n_homology_zones.times do |j| n_homology_zones.times do |i| next if i == j if distances[j][i] >= ident_thresold delete_positions << j delete_positions << i end end end delete_positions.uniq! return delete_positions end |
#set_cut_positions(homology_zones) ⇒ Object
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 46 def set_cut_positions(homology_zones) cut_positions = [] last_cut = -1 homology_zones.each_with_index do |hom_zone, i| if i > 0 positions = [] positions << last_cut + 1 # Start of fragment cut_position = homology_zones[i-1][STOP] + (hom_zone[BEG] - homology_zones[i-1][STOP])/2 positions << cut_position # End of fragment last_cut = cut_position cut_positions << positions end end cut_positions << [last_cut, homology_zones.last[HIT].first.q_len-1] return cut_positions end |
#set_limits(hit, q_beg, q_end, s_beg, s_end) ⇒ Object
287 288 289 290 291 292 |
# File 'lib/full_lengther_next/classes/chimeric_seqs.rb', line 287 def set_limits(hit, q_beg, q_end, s_beg, s_end) hit.q_beg = q_beg hit.q_end = q_end hit.s_beg = s_beg hit.s_end = s_end end |