#!/usr/bin/perl -w # # compare_sequences.pl # A script to compare two sequences and output the position and base change at each position # # janet.higgins@bbsrc.ac.uk, 22/02/2011 # # Script used to compare the cured A and C versions of the naive reference sequences # The unigene identifier needs to be the same in both sequences # can also be used to compare a cured sequence to the naive reference sequence # Input the sequences to be compare as fasta files # Outputs a list of the unigene,position within the unigene,and the base in both versions of the sequence at this position use FileHandle; # Hashes my %A=(); my %C=(); # Check arguments unless(@ARGV == 3) { print "Usage: compare_sequences.pl \n"; exit; } # Open FHs open(A, "<$ARGV[0]") or die "Couldn't open inputA fasta file ($!)\n"; open(C, "<$ARGV[1]") or die "Couldn't open inputC fasta file ($!)\n"; open(OUT, ">$ARGV[2]") or die "Couldn't open output file ($!)\n"; #Read in genome A refseqs print "Reading in A genome cured refseq... "; while () { chomp; if (/^>(\S+)(.*)/) { ($id)= $1; } elsif (/^(\w+)$/) { $A{$id} .=$1; } } close A; print "done\n"; #Read in genome C refseqs print "Reading in C genome cured refseq... "; while () { chomp; if (/^>(\S+)(.*)/) { ($id)= $1; } elsif (/^(\w+)$/) { $C{$id} .=$1; } } close C; print "done\n"; foreach my $id (sort keys %A) { my @A = split //, $A{$id}; my @C = split //, $C{$id}; foreach my $pos (0 ... $#A) { if ( $A[$pos] ne $C[$pos]){ print OUT "$id\t"; print OUT "" . $pos+1 . "\t"; print OUT "$A[$pos]\t"; print OUT "$C[$pos]\t"; print OUT "\n"; } } } close OUT;