alignVariation

#!/usr/bin/env ruby

require 'optimist'
require 'bio'

ARGV.push("--help") if ARGV.empty?
opts = Optimist::options do
 banner File.basename($0)
 opt :input, "input file", :required=>true, :type=>:string
 opt :threshold, "variation fraction threshold", :default=>0.2
end

seqs = []
Bio::FlatFile.new(Bio::FastaFormat, File.new(opts.input)).each do |seq|
 seqs.push(seq.seq)
end

printf("Position\tVariation\n")
diffs = [0]*seqs.first.length
0.upto(diffs.length - 1) do |i|
 bases = []
 seqs.each do |seq|
 bases.push(seq[i])
 end
 counts = {}
 bases.uniq.each do |b|
 counts[b] = bases.count(b)
 end
 variation = (1000*(seqs.length-counts.values.max)/seqs.length)/1000.0
 if variation > opts.threshold
 printf("%3d\t%0.3f\n", i+1, variation)
 end
end