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
39
40
41
42
43
44
45
46
47
48
|
# File 'bin/acnfp-trans-counter', line 8
def analyse(filename)
alignments = Bio::FlatFile.open(filename)
alignments.to_a.combination(2) do |a,b|
transition_count = 0
transversion_count = 0
indel_count = 0
informative_count = 0
alignment = Bio::Alignment::OriginalAlignment.new([a,b])
alignment.each_site do |site|
next if site.has_gap?
if site[0] == site[1] or site.include? "N"
informative_count += 1
next
end
case site.sort.map{|char| char.upcase}
when %w{C T}, %w{A G}
transition_count += 1
when %w{A C}, %w{A T}, %w{C G}, %w{G T}
transversion_count += 1
else
raise "Unknown pairing: #{site}"
end
end
mutation_count = transversion_count + transition_count.to_f
= "%s vs %s" % [a.definition, b.definition]
puts
puts "=" * .length
puts "Transversions: %d (%.2f%%)" % [transversion_count, 100 * transversion_count / mutation_count]
puts "Transitions: %d (%.2f%%)" % [transition_count, 100 * transition_count / mutation_count]
puts "Transversions / Transitions: %.2f" % [transversion_count.to_f / transition_count]
puts "Non-indel sites: %d" % informative_count
puts
end
end
|