User:WikiFew/BE131

=hw1.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program reads a FASTA file with DNA sequences.  It will then output the reverse complement of this sequence in either DNA or RNA, as specified by the user of the program.

open FILE, $ARGV[0] or die "This is not a valid file; try again\n";

while ($line = ) {   chomp($line); if ($line =~ />/) {       if (defined $len) { # i.e. if there was a previous sequence print "$name, $len bp\n"; $reverse = reverse(lc($origDNA)); if ($ARGV[1] eq "dna") {               $reverse =~ tr/acgt/tgca/; }           elsif ($ARGV[1] eq "rna") {               $reverse =~ tr/acgt/ugca/; }           print "$reverse\n"; }       # initialize variables $name = $line; $len = 0; $origDNA = ""; }    else {       $origDNA .= $line; $len += length($line); } }

print "$name, $len bp\n"; $reverse = reverse(lc($origDNA)); if ($ARGV[1] eq "dna") {   $reverse =~ tr/acgt/tgca/; } elsif ($ARGV[1] eq "rna") {   $reverse =~ tr/acgt/ugca/; } print "$reverse\n";
 * 1) This section of code is used for the final sequence in the file

close FILE;

=hw2.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program compares two FASTA files, printing sequence names that are present in one file, but not the other as well as sequence names that are present in both but refer to different sequences.
 * 2) Syntax for use: perl -w hw2.pl filename1 filename2

open(FILE1, $ARGV[0]) or die("Please input a functional file.\n"); open(FILE2, $ARGV[1]) or die("Please input a functional file.\n");
 * 1) Error handling- if file does not exist or if file has duplicate sequence names

$filename1 = $ARGV[0]; $filename2 = $ARGV[1];

open(FILE1, $filename1); open(FILE2, $filename2);

if ($filename1 eq $filename2) {   exit; } else {   %sequences1 = make_hash($filename1); %sequences2 = make_hash($filename2); @names1 = keys(%sequences1); @names2 = keys(%sequences2); # Comparing sequence names @same_names = ""; @in1not2 = ""; @in2not1 = ""; $length1 = @names1 + 0; $length2 = @names2 + 0; # Determining what is in file 1, but not in file 2 for ($j = 0; $j < $length2; ++$j) {       for ($i = 0; $i < $length1; ++$i) {           if ($names1[$i] =~ /$names2[$j]/g) {               push @same_names, $names1[i]; }           else {               push @in1not2, $names1[i]; }       }    }    # Printing what is in file 1, but not in file 2 print "File 1 but not file 2: \n"; for ($k = 0; $k < @in1not2 + 0; ++$k) {       print "$in1not2[$i] \n"; }   # Determining what is in file 2, but not in file 1 for ($j = 0; $j < $length1; ++$j) {       for ($i = 0; $i < $length2; ++$i) {           if (!($names2[$i] =~ /$names1[$j]/g)) {               push @in2not1, $names2[i]; }       }    }    # Printing what is in file 2, but not in file 1 print "File 2 but not file 1: \n"; for ($k = 0; $k < @in1not2 + 0; ++$k) {       print "$in2not1[$i] \n"; }   # If both files contain duplicate sequence names if (@in1not2 + 0 == 1 && @in2not1 + 0 == 1) {       die("Files contain duplicat sequence names.\n"); }   # If two sequences names are the same, but have different sequences @diff_sequences = diff_seq(@same_names); print "Both files, different sequences:" for ($k = 0; $k < @in1not2 + 0; ++$k) {       print "$diff_sequences[$i] \n"; }
 * 1) If files are equivalent, exit without printing
 * 1) Else, print sequence names that are present in one, but not the other or present in both but with different sequences

}

sub make_hash {   my ($filename) = @_; my (%name2seq, $name, $seq); open FILE, $filename; while () {       chomp; if (/>/) {           s/>//; if (defined $name) {               $name2seq{$name} = $seq; }           $name = $_; $seq = ""; }        else {           $seq .= $_; }   }    $name2seq{$name} = $seq; close FILE; return %name2seq; }
 * 1) This subroutine will put all sequence names into a hash

sub diff_seq {   my (@names) = @_; @different = ""; for ($j = 0; $j < @names + 0; ++$j) {       if (!($sequences1{$j} =~ /$sequences2{$j}/gi)) {           push @different, @names[j]; }   }    return @different; }
 * 1) This subroutine will determine if two sequence names contain different sequences

close (FILE1); close (FILE2);

=hw3a.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program reads a FASTA file supplied by the user and tests three conditions. It will output one line per sequence containing the sequence name, the length, and a description of whether or not it passed the tests.
 * 2) Syntax for use: perl -w hw3a.pl filename

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

while($line = ) {   chomp $line; if ($line =~ /^>gi\|[0-9]+\|(.+)/) {       if (defined($name)) {           $status = testing($sequence); $length = length($sequence); print "$name $length bp $status\n"; }       $name = $line; $name =~ s/>//; $sequence = ""; }   else {       $sequence .= $line; } } $status = testing($sequence); $length = length($sequence); print "$name $length bp $status\n";

sub testing {   my($length) = length($sequence); if ($sequence =~ /[bd-fh-sv-z0-9]/i) {       # if a foreign character is present $status = "FAIL"; }   elsif ($length%3 != 0) {       # if the sequence length is not divisible by 3 $status = "FAIL"; }   else {       # searching for a stop codon while ($sequence =~ /(TAG|TGA|TAA)/gi) {           $position = pos($sequence) - 2; if ($position % 3 == 1) # to ensure it is an actual codon {               $stop_positions .= " $position"; }       }        if (defined($stop_positions)) {           $status = "STOP $stop_positions"; }       else {           $status = "PASS"; }   }    return $status; }
 * 1) This subroutine tests the sequence for three conditions.  It determines if all characters are nucleotides, the length is divisible by three, and if any stop codons are present.

=hw3b.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program reads a FASTA file supplied by the user and determines the ratio of CpG divided by the product of the probability of a C and the probability of a G.
 * 2) Syntax for use: perl -w hw3b.pl filename


 * 1) This ratio is useful because two different genes will contain different CG content.  However, as two species begin to diverge, their DNA content will be modified slightly, thus changing their CG content.  Therefore, we can use this probability to determine the homology of two genes or to distinguish between bacterial and mammalian sequences.

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

while($line = ) {   chomp $line; if ($line =~ /^>gi\|[0-9]+\|(.+)/) {       if (defined($name)) {           $ratio = probability($sequence); print "$name\n", "$ratio\n"; }       $name = $line; $name =~ s/>//; $sequence = ""; }   else {       $sequence .= $line; } } $ratio = probability($sequence); print "$name\n", "$ratio\n";

sub probability {   $length = length($sequence); #initializing variables $C_counter = 0; $G_counter = 0; $CG_counter = 0;
 * 1) This subroutine calculates the probability of a nucleotide being a C or a G and the probability of a dinucleotide being a CG.  It then computes a ratio of those probabilities.

while ($sequence =~ /C/gi) {       $C_counter++; }   while ($sequence =~ /G/gi) {       $G_counter++; }   while ($sequence =~ /CG/gi) {       $CG_counter++; }   $prob_C = $C_counter / $length; $prob_G = $G_counter / $length; #for x nucleotides, there are (x-1) possible dinucleotides $prob_CG = $CG_counter / ($length - 1); $ratio = $prob_CG/($prob_C * $prob_G); return $ratio; }

=hw4.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program uses the Needleman-Wunsch algorithm to determine the best alignment of two sequences given a substitution matrix and the gap score.  The output will be the Score Matrix as well as the best alignment.
 * 2) Syntax for use: perl -w hw4.pl matrixfilename sequencesfilename

open(FILE1, $ARGV[0]) or die("Please input a functional file.\n"); open(FILE2, $ARGV[1]) or die("Please input a functional file.\n");

$row = 0; # Initializing a row counter while ($line = ) {   if ($line =~ /[a-z]/i) {       $matrix = uc($line); $matrix =~ s/ //g; chomp($matrix); @matrix = split //, $matrix; }   elsif ($row < length($matrix)) {       $col = 0; # Initializing a column counter while ($line =~ /([-]*\d+)/g) {           $A = $matrix[$row]; $B = $matrix[$col]; $key = $A.$B; $hash{$key} = $1; $col++; }       $row++; }   else {       $gap = $line; chomp($gap); }           }
 * 1) Reading the matrix file

$i = 1; # Initializing a control variable that will stop after 2 sequences while($line = ) {
 * 1) Reading the file of sequences

last if ($i==0); chomp $line; if ($line =~ /^>/) {       if (!defined($name1)) {           $name1 = $line; $name1 =~ s/>//; $sequence1 = ""; }       elsif (!defined($name2)) {           $name2 = $line; $name2 =~ s/>//; $sequence2 = ""; }       else {           $i = 0; }   }    else {       if (defined($name2)) {           $sequence2 .= uc($line); }       elsif (defined($name1)) {           $sequence1 .= uc($line); }   } } @sequence1 = split //, $sequence1; @sequence2 = split //, $sequence2;

if ($sequence1 =~ /[^$matrix]i*/) {   die "Sequence 1 is not valid.\n"; } if ($sequence2 =~ /[^$matrix]i*/) {   die "Sequence 2 is not valid.\n"; }
 * 1) Making sure the sequences contain only characters present in the matrix file

$ScoreMatrix[0][0] = "*"; $ScoreMatrix[0][1] = "*"; for($i = 0; $i < length($sequence1); $i++) {   $ScoreMatrix[0][2+$i] = $sequence1[$i]; }
 * 1) Making the first row of the Score Matrix using Sequence 1

$ScoreMatrix[1][0] = "*"; for($i = 0; $i < length($sequence2); $i++) {   $ScoreMatrix[2+$i][0] = $sequence2[$i]; }
 * 1) Making the first column of the Score Matrix using Sequence 2

$ScoreMatrix[1][1] = 0; for($i = 0; $i < length($sequence1); $i++) {   $ScoreMatrix[1][2+$i] = $ScoreMatrix[1][1+$i] + $gap; } for($i = 0; $i < length($sequence2); $i++) {   $ScoreMatrix[2+$i][1] = $ScoreMatrix[1+$i][1] + $gap; }
 * 1) Establishing the Boundary Conditions for the Score Matrix

$most_digits = 1; #variable that marks the longest score to help build a square Score Matrix for($r = 2; $r <= length($sequence2)+1; $r++) {   for ($c = 2; $c <= length($sequence1)+1; $c++) {       $above = $ScoreMatrix[$r-1][$c] + $gap; $left = $ScoreMatrix[$r][$c-1] + $gap; $diagonal = $ScoreMatrix[$r-1][$c-1] + $hash{$sequence1[$c-2].$sequence2[$r-2]}; $score = find_max; $ScoreMatrix[$r][$c] = $score; if (length($score) > length($most_digits)) #comparing the longest score {           $most_digits = $score; }   } }
 * 1) Building the Score Matrix

for($r = 0; $r <= length($sequence2)+1; $r++) {   print "|"; for ($c = 0; $c <= length($sequence1)+1; $c++) {       $spaces = length($most_digits) - length($ScoreMatrix[$r][$c]); for ($i = 0; $i < $spaces; $i++) {           print " "; }       print "$ScoreMatrix[$r][$c]|"; }   print "\n"; }
 * 1) Printing the Score Matrix

$flag = 1; $c = length($sequence1); $r = length($sequence2); $aligned1 = ""; $aligned2 = "";
 * 1) Determining the alignment

while ($flag == 1) {   if ($r == 0) {       if ($c == 0) {           $flag = 0; }       else {           $aligned1 = $sequence1[$c-1].$aligned1; $aligned2 = "-".$aligned2; $c--; }   }    elsif ($c == 0) {       if ($r == 0) {           $flag = 0; }       else {           $aligned1 = "-".$aligned1; $aligned2 = $sequence2[$r-1].$aligned2; $r--; }   }    else {       if ($ScoreMatrix[$r+1][$c+1] == $ScoreMatrix[$r][$c+1] + $gap) {           $aligned1 = "-".$aligned1; $aligned2 = $sequence2[$r-1].$aligned2; $r--; }       elsif ($ScoreMatrix[$r+1][$c+1] == $ScoreMatrix[$r+1][$c] + $gap) {           $aligned1 = $sequence1[$c-1].$aligned1; $aligned2 = "-".$aligned2; $c--; }       else {           $aligned1 = $sequence1[$c-1].$aligned1; $aligned2 = $sequence2[$r-1].$aligned2; $r--; $c--; }   } }

print "$name1 aligned to $name2\n"; for ($j = 0; $j < length($aligned1)/80; $j++) {   print substr($aligned1, 80 * $j, 80), "\n"; print substr($aligned2, 80 * $j, 80), "\n\n"; }
 * 1) Printing the final alignment

sub find_max {   if ($above >= $left) {       if ($above >= $diagonal) {           $max = $above; }   }    if ($left >= $above) {       if ($left >= $diagonal) {           $max = $left; }   }    if ($diagonal >= $above) {       if ($diagonal >= $left) {           $max = $diagonal; }   }    return $max; }
 * 1) This subroutine will determine the maximum value of three numbers

close (FILE1); close (FILE2);

matrix.txt
A R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X 4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -1 5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4

proteing_seq.pl
>First protein Sequence AERTWPPRPTIL

>Second protein Sequence AQRTWPGRPTIL

>Third protein sequence AQRTWPRPTIL

references.pl

 * 1) ! /usr/bin/perl -w

@AoA = ([1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]); $AoA_ref = \@AoA;

print "$AoA[0][1]\n"; print "$AoA_ref->[2][1]\n";

$ArrayOfArrays[0][0] = 5; print "$ArrayOfArrays[0][0]\n";

$sequence1 = "ACGTT"; @sequence1 = split //, $sequence1;

$ScoreMatrix[0][0] = "*"; $ScoreMatrix[0][1] = "*"; for($i = 0; $i < length($sequence1); $i++) {   $ScoreMatrix[0][2+$i] = $sequence1[$i]; }

print "$ScoreMatrix[0][6]\n";

=hw5.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program calculates entropy and mutual information of a multiple sequence alignment.
 * 2) Syntax for use: perl -w hw5.pl Stockholmfilename

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

$number_of_rows = 0;
 * 1) Initializing variables

while ($line = ) {   chomp $line; if (!($line =~ /^\#/)) {       if (!($line =~ /^[\/]{2}/)) {           if ($line =~ /[\S]+[\s]+(.+)/) {               @temp_sequence = split //, $1; for ($number_of_cols = 0; $number_of_cols < length($1); $number_of_cols++) {                   $SequenceMatrix[$number_of_rows][$number_of_cols] = $temp_sequence[$number_of_cols];
 * 1) Reading and working with the file.

}               $number_of_rows++ }       }    } } $number_of_rows--; $number_of_cols--;


 * 1) Calculating the probability distribution for a column and the entropy of that distribution.

for ($col = 0; $col < $number_of_cols; $col++) {   $A_counter = 0; $C_counter = 0; $G_counter = 0; $U_counter = 0; $gap_counter = 0; for ($row = 0; $row < $number_of_rows; $row++) {       if ($SequenceMatrix[$row][$col] eq "A") {           $A_counter++; }       elsif ($SequenceMatrix[$row][$col] eq "C") {           $C_counter++; }       elsif ($SequenceMatrix[$row][$col] eq "G") {           $G_counter++; }       elsif ($SequenceMatrix[$row][$col] eq "U") {           $U_counter++; }       elsif ($SequenceMatrix[$row][$col] eq ".") {           $gap_counter++; }   }    for ($row = 0; $row < $number_of_rows; $row++) {       $probability[$row][$col]{"A"} = $A_counter/($number_of_rows); $probability[$row][$col]{"C"} = $C_counter/($number_of_rows); $probability[$row][$col]{"G"} = $G_counter/($number_of_rows); $probability[$row][$col]{"U"} = $U_counter/($number_of_rows); $probability[$row][$col]{"gap"} = $gap_counter/($number_of_rows); $entropy[$row][$col] = entropy; } }

$i = 0; for ($row = 0; $row < $number_of_rows; $row++) {   for ($col = 0; $col < $number_of_cols; $col++) {       $all_entropy[$i] = $entropy[$row][$col]; $i++; } } @low_entropy = sort {$a <=> $b} @all_entropy;
 * 1) Sorting the entropy

$i = 0; print "The 'bottom 10' columns with lowest entropy:\n"; for ($r = 0; $r < $number_of_rows; $r++) {   for ($c = 0; $c < $number_of_cols; $c++) {       if ($low_entropy[$i] == $entropy[$r][$c]) {           if ($i < 10) {               print "Column $c with entropy $entropy[$r][$c]\n"; $i++; }       }    } }
 * 1) Matching each bottom entropy with the corresponding column

for ($col1 = 0; $col1 < $number_of_cols; $col1++) {   for ($row = 0; $row < $number_of_rows; $row++) {       $A_counter = 0; $C_counter = 0; $G_counter = 0; $U_counter = 0; $gap_counter = 0; $A_counter2 = 0; $C_counter2 = 0; $G_counter2 = 0; $U_counter2 = 0; $gap_counter2 = 0; $AA_counter = 0; $AC_counter = 0; $AG_counter = 0; $AU_counter = 0; $Agap_counter = 0; $CA_counter = 0; $CC_counter = 0; $CG_counter = 0; $CU_counter = 0; $Cgap_counter = 0; $GA_counter = 0; $GC_counter = 0; $GG_counter = 0; $GU_counter = 0; $Ggap_counter = 0; $UA_counter = 0; $UC_counter = 0; $UG_counter = 0; $UU_counter = 0; $Ugap_counter = 0; $gapA_counter = 0; $gapC_counter = 0; $gapG_counter = 0; $gapU_counter = 0; $gapgap_counter = 0; if ($SequenceMatrix[$row][$col1] eq "A") {           $A_counter++; for ($col2 = $col1 + 1; $col2 < $number_of_cols; $col2++) {               if ($SequenceMatrix[$row][$col2] eq "A") {                   $A_counter2++; $AA_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "C") {                   $C_counter2++; $AC_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "G") {                   $G_counter2++; $AG_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "U") {                   $U_counter2++; $AU_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq ".") {                   $gap_counter2++; $Agap_counter++; }           }        }        elsif ($SequenceMatrix[$row][$col1] eq "C") {           $C_counter++; for ($col2 = $col1 + 1; $col2 < $number_of_cols; $col2++) {               if ($SequenceMatrix[$row][$col2] eq "A") {                   $A_counter2++; $CA_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "C") {                   $C_counter2++; $CC_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "G") {                   $G_counter2++; $CG_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "U") {                   $U_counter2++; $CU_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq ".") {                   $gap_counter2++; $Cgap_counter++; }           }        }        elsif ($SequenceMatrix[$row][$col1] eq "G") {           $G_counter++; for ($col2 = $col1 + 1; $col2 < $number_of_cols; $col2++) {               if ($SequenceMatrix[$row][$col2] eq "G") {                   $A_counter2++; $GA_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "C") {                   $C_counter2++; $GC_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "G") {                   $G_counter2++; $GG_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "U") {                   $U_counter2++; $GU_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq ".") {                   $gap_counter2++; $Ggap_counter++; }           }        }        elsif ($SequenceMatrix[$row][$col1] eq "U") {           $U_counter++; for ($col2 = $col1 + 1; $col2 < $number_of_cols; $col2++) {               if ($SequenceMatrix[$row][$col2] eq "A") {                   $A_counter2++; $UA_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "C") {                   $C_counter2++; $UC_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "G") {                   $G_counter2++; $UG_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "U") {                   $U_counter2++; $UU_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq ".") {                   $gap_counter2++; $Ugap_counter++; }           }        }        elsif ($SequenceMatrix[$row][$col1] eq ".") {           $gap_counter++; for ($col2 = $col1 + 1; $col2 < $number_of_cols; $col2++) {               if ($SequenceMatrix[$row][$col2] eq "A") {                   $A_counter2++; $gapA_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "C") {                   $C_counter2++; $gapC_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "G") {                   $G_counter2++; $gapG_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq "U") {                   $U_counter2++; $gapU_counter++; }               elsif ($SequenceMatrix[$row][$col2] eq ".") {                   $gap_counter2++; $gapgap_counter++; }           }        }        $mutual[$row][$col1] = mutual; $mutual_hash{"Columns $col1 and $col2"} = $mutual[$row][$col1]; } }
 * 1) Determining the probability for a pair and the mutual information

@sorted = sort { $mutual_hash{$a} cmp $mutual_hash{$b} } keys %mutual_hash; print "The 'top 20' columns with highest mutual information:\n"; $i = 0; while ($i < 20) {   $key = $sorted[$i]; print "$key\n"; $i++; }
 * 1) Sorting the mutual information

print "The 'top 50' columns with highest mutual information:\n"; $i = 0; while ($i < 50) {   $key = $sorted[$i]; print "$key\n"; $i++; }

sub entropy {   if ($probability[$row][$col]{"A"} != 0) {       $S_A = -$probability[$row][$col]{"A"} * (log($probability[$row][$col]{"A"})/log(2)); }   else {       $S_A = 0; }   if ($probability[$row][$col]{"C"} != 0) {       $S_C = -$probability[$row][$col]{"C"} * (log($probability[$row][$col]{"C"})/log(2)); }   else {       $S_C = 0; }   if ($probability[$row][$col]{"G"} != 0) {       $S_G = -$probability[$row][$col]{"G"} * (log($probability[$row][$col]{"G"})/log(2)); }   else {       $S_G = 0; }   if ($probability[$row][$col]{"U"} != 0) {       $S_U = -$probability[$row][$col]{"U"} * (log($probability[$row][$col]{"U"})/log(2)); }   else {       $S_U = 0; }   if ($probability[$row][$col]{"gap"} != 0) {       $S_gap = -$probability[$row][$col]{"gap"} * (log($probability[$row][$col]{"gap"})/log(2)); }   else {       $S_gap = 0; }   $S = $S_A + $S_C + $S_G + $S_U + $S_gap; return $S; }
 * 1) This subroutine calculates the entropy of the probability distribution

sub mutual {   $joint_AA = $AA_counter/$number_of_cols; $prob_A = $A_counter/$number_of_cols; $prob_A = $A_counter2/$number_of_cols; if ($joint_AA == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   else {       $M1 = $joint_AA * log($joint_AA/($prob_A*$prob_A)); }   $joint_AC = $AC_counter/$number_of_cols; $prob_A = $A_counter/$number_of_cols; $prob_C = $C_counter2/$number_of_cols; if ($joint_AC == 0) {       $M2 = 0; }   elsif ($prob_A == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   else {       $M2 = $joint_AC * log($joint_AC/($prob_A*$prob_C)); }   $joint_AG = $AG_counter/$number_of_cols; $prob_A = $A_counter/$number_of_cols; $prob_G = $G_counter2/$number_of_cols; if ($joint_AG == 0) {       $M3 = 0; }   elsif ($prob_A == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   else {       $M3 = $joint_AG * log($joint_AG/($prob_A*$prob_G)); }   $joint_AU = $AU_counter/$number_of_cols; $prob_A = $A_counter/$number_of_cols; $prob_U = $U_counter2/$number_of_cols; if ($joint_AU == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   else {       $M4 = $joint_AU * log($joint_AU/($prob_A*$prob_U)); }   $joint_Agap = $Agap_counter/$number_of_cols; $prob_A = $A_counter/$number_of_cols; $prob_gap = $gap_counter2/$number_of_cols; if ($joint_Agap == 0) {       $M5 = 0; }   elsif ($prob_A == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   else {       $M5 = $joint_Agap * log($joint_Agap/($prob_A*$prob_gap)); }   $M_A = $M1 + $M2 + $M3 + $M4 + $M5;
 * 1) This subroutine calculates mutual information

$joint_CA = $CA_counter/$number_of_cols; $prob_C = $C_counter/$number_of_cols; $prob_A = $A_counter2/$number_of_cols; if ($joint_CA == 0) {       $M1 = 0; }   elsif ($prob_C == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   else {       $M1 = $joint_CA * log($joint_CA/($prob_C*$prob_A)); }   $joint_CC = $CC_counter/$number_of_cols; $prob_C = $C_counter/$number_of_cols; $prob_C = $C_counter2/$number_of_cols; if ($joint_CC == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   else {       $M2 = $joint_CC * log($joint_CC/($prob_C*$prob_C)); }   $joint_CG = $CG_counter/$number_of_cols; $prob_C = $C_counter/$number_of_cols; $prob_G = $G_counter2/$number_of_cols; if ($joint_CG == 0) {       $M3 = 0; }   elsif ($prob_C == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   else {       $M3 = $joint_CG * log($joint_CG/($prob_C*$prob_G)); }   $joint_CU = $CU_counter/$number_of_cols; $prob_C = $C_counter/$number_of_cols; $prob_U = $U_counter2/$number_of_cols; if ($joint_CU == 0) {       $M4 = 0; }   elsif ($prob_C == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   else {       $M4 = $joint_CU * log($joint_CU/($prob_C*$prob_U)); }   $joint_Cgap = $Cgap_counter/$number_of_cols; $prob_C = $C_counter/$number_of_cols; $prob_gap = $gap_counter2/$number_of_cols; if ($joint_Cgap == 0) {       $M5 = 0; }   elsif ($prob_C == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   else {       $M5 = $joint_Cgap * log($joint_Cgap/($prob_C*$prob_gap)); }   $M_C = $M1 + $M2 + $M3 + $M4 + $M5;

$joint_GA = $GA_counter/$number_of_cols; $prob_G = $G_counter/$number_of_cols; $prob_A = $A_counter2/$number_of_cols; if ($joint_GA == 0) {       $M1 = 0; }   elsif ($prob_G == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   else {       $M1 = $joint_GA * log($joint_GA/($prob_G*$prob_A)); }   $joint_GC = $GC_counter/$number_of_cols; $prob_G = $G_counter/$number_of_cols; $prob_C = $C_counter2/$number_of_cols; if ($joint_GC == 0) {       $M2 = 0; }   elsif ($prob_G == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   else {       $M2 = $joint_GC * log($joint_GC/($prob_G*$prob_C)); }   $joint_GG = $GG_counter/$number_of_cols; $prob_G = $G_counter/$number_of_cols; $prob_G = $G_counter2/$number_of_cols; if ($joint_GG == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   else {       $M3 = $joint_GG * log($joint_GG/($prob_G*$prob_G)); }   $joint_GU = $GU_counter/$number_of_cols; $prob_G = $G_counter/$number_of_cols; $prob_U = $U_counter2/$number_of_cols; if ($joint_GU == 0) {       $M4 = 0; }   elsif ($prob_G == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   else {       $M4 = $joint_GU * log($joint_GU/($prob_G*$prob_U)); }   $joint_Ggap = $Ggap_counter/$number_of_cols; $prob_G = $G_counter/$number_of_cols; $prob_gap = $gap_counter2/$number_of_cols; if ($joint_Ggap == 0) {       $M5 = 0; }   elsif ($prob_G == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   else {       $M5 = $joint_Ggap * log($joint_Ggap/($prob_G*$prob_gap)); }   $M_G = $M1 + $M2 + $M3 + $M4 + $M5;

$joint_UA = $UA_counter/$number_of_cols; $prob_U = $U_counter/$number_of_cols; $prob_A = $A_counter2/$number_of_cols; if ($joint_UA == 0) {       $M1 = 0; }   elsif ($prob_U == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   else {       $M1 = $joint_UA * log($joint_UA/($prob_U*$prob_A)); }   $joint_UC = $UC_counter/$number_of_cols; $prob_U = $U_counter/$number_of_cols; $prob_C = $C_counter2/$number_of_cols; if ($joint_UC == 0) {       $M2 = 0; }   elsif ($prob_U == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   else {       $M2 = $joint_UC * log($joint_UC/($prob_U*$prob_C)); }   $joint_UG = $UG_counter/$number_of_cols; $prob_U = $U_counter/$number_of_cols; $prob_G = $G_counter2/$number_of_cols; if ($joint_UG == 0) {       $M3 = 0; }   elsif ($prob_U == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   else {       $M3 = $joint_UG * log($joint_UG/($prob_U*$prob_G)); }   $joint_UU = $UU_counter/$number_of_cols; $prob_U = $U_counter/$number_of_cols; $prob_U = $U_counter2/$number_of_cols; if ($joint_UU == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   else {       $M4 = $joint_UU * log($joint_UU/($prob_U*$prob_U)); }   $joint_Ugap = $Ugap_counter/$number_of_cols; $prob_U = $U_counter/$number_of_cols; $prob_gap = $gap_counter2/$number_of_cols; if ($joint_Ugap == 0) {       $M5 = 0; }   elsif ($prob_U == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   else {       $M5 = $joint_Ugap * log($joint_Ugap/($prob_U*$prob_gap)); }   $M_U = $M1 + $M2 + $M3 + $M4 + $M5;

$joint_gapA = $gapA_counter/$number_of_cols; $prob_gap = $gap_counter/$number_of_cols; $prob_A = $A_counter2/$number_of_cols; if ($joint_gapA == 0) {       $M1 = 0; }   elsif ($prob_gap == 0) {       $M1 = 0; }   elsif ($prob_A == 0) {       $M1 = 0; }   else {       $M1 = $joint_gapA * log($joint_gapA/($prob_gap*$prob_A)); }   $joint_gapC = $gapC_counter/$number_of_cols; $prob_gap = $gap_counter/$number_of_cols; $prob_C = $C_counter2/$number_of_cols; if ($joint_gapC == 0) {       $M2 = 0; }   elsif ($prob_gap == 0) {       $M2 = 0; }   elsif ($prob_C == 0) {       $M2 = 0; }   else {       $M2 = $joint_gapC * log($joint_gapC/($prob_gap*$prob_C)); }   $joint_gapG = $gapG_counter/$number_of_cols; $prob_gap = $gap_counter/$number_of_cols; $prob_G = $G_counter2/$number_of_cols; if ($joint_gapG == 0) {       $M3 = 0; }   elsif ($prob_gap == 0) {       $M3 = 0; }   elsif ($prob_G == 0) {       $M3 = 0; }   else {       $M3 = $joint_gapG * log($joint_gapG/($prob_gap*$prob_G)); }   $joint_gapU = $gapU_counter/$number_of_cols; $prob_gap = $gap_counter/$number_of_cols; $prob_U = $U_counter2/$number_of_cols; if ($joint_gapU == 0) {       $M4 = 0; }   elsif ($prob_gap == 0) {       $M4 = 0; }   elsif ($prob_U == 0) {       $M4 = 0; }   else {       $M4 = $joint_gapU * log($joint_gapU/($prob_gap*$prob_U)); }   $joint_gapgap = $gapgap_counter/$number_of_cols; $prob_gap = $gap_counter/$number_of_cols; $prob_gap = $gap_counter2/$number_of_cols; if ($joint_gapgap == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   elsif ($prob_gap == 0) {       $M5 = 0; }   else {       $M5 = $joint_gapgap * log($joint_gapgap/($prob_gap*$prob_gap)); }   $M_gap = $M1 + $M2 + $M3 + $M4 + $M5;

$M = $M_A +$M_C +$M_G +$M_U +$M_gap; return $M; } close(FILE);

=hw6a.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program calculates the probability that a method of randomly selecting a sequence was used given a sequence.  The two possible methods are 1) selecting nucleotides at random, and 2) selecting nucleotides at random, but rechoosing if a stop codon is selected.
 * 2) Syntax for use: perl -w hw6a.pl filename


 * 1) The log-ratio tells you a relation between the two posterior probabilities that just the posterior probability P(G=1|S) does not tell you.  If you wanted to compare the two probabilities, the log-ratio would be more useful.

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

while ($line = ) {   chomp $line; if ($line =~ /^>/) {       if (defined($name)) {           probability; }       $name = $line; $name =~ s/>//; $A_counter = 0; $C_counter = 0; $G_counter = 0; $U_counter = 0; $total_length = 0; $sequence = ""; $flag = 0; }   elsif ($line =~ /\S/) {       $sequence .= $line; @sequence = split //, $sequence; $total_length += length($sequence); while ($line =~ /[A]/gi) {           $A_counter++; }       while ($line =~ /[C]/gi) {           $C_counter++; }       while ($line =~ /[G]/gi) {           $G_counter++; }       while ($line =~ /[U]/gi) {           $U_counter++; }       if ($line =~ /UAA|UAG|UGA/g) {           $flag = 1; }   } } probability;
 * 1) Reading the file

sub probability {   $prob = 1; if ($flag == 0) {       foreach $x (@sequence) {           if ($x eq "A") {               $prob *= ($A_counter/$total_length); }           elsif ($x eq "C") {               $prob *= ($C_counter/$total_length); }           elsif ($x eq "G") {               $prob *= ($G_counter/$total_length); }           elsif ($x eq "U") {               $prob *= ($U_counter/$total_length); }       }
 * 1) This subroutine calculates the posterior probability and log-ratio and prints the results.

$probG1 = (0.5*$prob)/(0.5*$prob + 0.5*$prob); # P(G=1|S) #                  P(G=1)P(S|G=1) #P(G=1|S) = --- #           P(G=1)P(S|G=1)+P(G=0)P(S|G=0)

$probG0 = (0.5*$prob)/(0.5*$prob + 0.5*$prob); # P(G=0|S) #                  P(G=0)P(S|G=0) #P(G=0|S) = --- #           P(G=0)P(S|G=0)+P(G=1)P(S|G=1) $log = log($probG1/$probG0)/log(2);

print "$name:\n"; print "posterior probability: P(G=1|S) = $probG1\n"; print "log-ratio = $log\n\n"; }   else {       $probG1 = 0; # P(G = 1|S) $probG0 = 1; # P(G = 0|S) print "$name:\n"; print "posterior probability: P(G=1|S) = $probG1\n"; print "log-ratio is undefined\n\n"; } }

close(FILE);

hw6a.pl again

 * 1) ! /usr/bin/perl -w


 * 1) This program calculates the probability that a method of randomly selecting a sequence was used given a sequence.  The two possible methods are 1) selecting nucleotides at random, and 2) selecting nucleotides at random, but rechoosing if a stop codon is selected.
 * 2) Syntax for use: perl -w hw6a.pl filename


 * 1) The log-ratio does not tell more because it can be derived from the posterior probability, however it is a useful way to compare and see the information.

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

while ($line = ) {   chomp $line; if ($line =~ /^>/) {       $name = $line; $name =~ s/>//; }   else {       $line = uc($line); $line =~ tr/U/T/; $sequences{$name} .= $line; } }
 * 1) Reading the file
 * 2) Basic file-handling; separating a sequence name from the sequence itself.  I converted the entire sequence to upper case to avoid issues with conflicting case.  In case an RNA sequence is input, all U's are converted to T's so the corresponding DNA strand is analyzed.

@nucleotides = ('A', 'T', 'G', 'C');

foreach $sequence (values %sequences) {   $total_length += length($sequence); for ($i = 0; $i < length($sequence); $i++) {       $counts{substr($sequence, $i, 1)}++; # counts is a hash that stores the number of occurrences of each nucleotide # e.g. counts{'A'} = 4 if A appeared 4 times } } foreach $nuc (@nucleotides) {   $prob{$nuc} = $counts{$nuc}/$total_length; # %prob stores the probability distribution of the nucleotides }
 * 1) Calculating p(x) for the nucleotides
 * 2) Sum the number of times a nucleotide appears in a sequence, then divide that number by the length.

$number_of_nucleotides = scalar(@nucleotides); $probability_of_stop = $prob{'T'} * $prob{'A'} * $prob{'A'}; $probability_of_stop += $prob{'T'} * $prob{'G'} * $prob{'A'}; $probability_of_stop += $prob{'T'} * $prob{'A'} * $prob{'G'};
 * 1) Calculating the probability of a stop codon

for ($i = 0; $i < $number_of_nucleotides; $i++) {   for ($j = 0; $j < $number_of_nucleotides; $j++) {       for ($k = 0; $k < $number_of_nucleotides; $k++) {           # If G = 0, it doesn't matter if a stop codon is drawn. # Therefore it is a straightforward calculation. # If G = 1, you have to generate a new codon if a stop codon is created. # Therefore you have to divide the probability of making some nonstop codon # by the probability of making a nonstop codon, (1 - probability of a stop). $probG0{$nucleotides[$i].$nucleotides[$j].$nucleotides[$k]} = $prob{$nucleotides[$i]} * $prob{$nucleotides[$j]} * $prob{$nucleotides[$k]}; $probG1{$nucleotides[$i].$nucleotides[$j].$nucleotides[$k]} = $prob{$nucleotides[$i]} * $prob{$nucleotides[$j]} * $prob{$nucleotides[$k]}; $probG1{$nucleotides[$i].$nucleotides[$j].$nucleotides[$k]} = $probG1{$nucleotides[$i].$nucleotides[$j].$nucleotides[$k]} / (1 - $probability_of_stop); }   } }

$probG1{'TAA'} = $probG1{'TGA'} = $probG1{'TAG'} = 0;
 * 1) The probability of a stop codon if G = 1 is 0 because, by definition, if a stop codon is generated it will be thrown out

foreach $key (keys %sequences) {   print "$key:\n"; $sequence = $sequences{$key}; $S_G1 = calculate($sequence, \%probG1); #P(S|G = 1) $S_G0 = calculate($sequence, \%probG0); #P(S|G = 0)
 * 1) Calculating posterior probabilities

$G1_S = $S_G1/($S_G1 + $S_G0); #P(G = 1|S) #                   P(G=1)P(S|G=1) # P(G=1|S) = --- #            P(G=1)P(S|G=1)+P(G=0)P(S|G=0) # Factor out P(G=0)=P(G=1) to simplify the calculation.

$G0_S = 1 - $G1_S; #P(G = 0|S) # There are only two ways to generate a sequence. # Therefore the P(G = 0) = 1 - P(G = 1).

print "  P(G = 1|S) = $G1_S\n";

# Calculating log ratio if ($G1_S > 0) {       $log_ratio = log($G1_S/$G0_S)/log(2); }   else {       $log_ratio = "undefined"; }   print "   log-ratio = $log_ratio\n"; }

sub calculate {   $seq = shift(@_); $prob_ref = shift(@_); $answer = 1; for($i = 0; $i < length($seq)/3; $i++) {       $answer *= $prob_ref->{substr($seq, $i*3, 3)}; }   return $answer; }
 * 1) This subroutine determines the probability a sequence was generated given a table with probabilities of codons.
 * 2) It takes a sequence with length divisible by 3 and analyzes each codon (i.e. each substring 3 characters long).  It takes each codon and uses that as a key to access a hash with the probability of making that codon (%probG0 or %probG1).  It takes the product of the probabilities of making each codon and returns the final product as the answer.

close(FILE);

=hw6b.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program will determine the probability a neuron is activated during a third experiment given data for an activated and a resting neuron,


 * 1) x = 0: resting
 * 2) x = 1: activated
 * 3) y = 0: non-spike
 * 4) y = 1: spike


 * 1) Syntax for use: perl -w hw6b.pl E1 E2 E3

open(FILE1, $ARGV[0]) or die("Please input a functional 'E1' file.\n"); open(FILE2, $ARGV[1]) or die("Please input a functional 'E2' file.\n"); open(FILE3, $ARGV[2]) or die("Please input a functional 'E3' file.\n");

$E1spikes = 0; while ($line = ) {   chomp $line; $E1spikes++; $E1_final_time = $line; } $y1x0 = $E1spikes/$E1_final_time; # P(y = 1|x = 0)
 * 1) Reading E1

$E2spikes = 0; while ($line = ) {   chomp $line; $E2spikes++; $E2_final_time = $line; } $y1x1 = $E2spikes/$E2_final_time; # P(y = 1|x = 1)
 * 1) Reading E2

$x1 = 0.01; #P(x=1)

close(FILE1); close(FILE2); close(FILE3);

hw6b.pl again

 * 1) ! /usr/bin/perl -w


 * 1) This program will determine the probability a neuron is activated during a third experiment given data for an activated and a resting neuron.


 * 1) x = 0: resting
 * 2) x = 1: activated
 * y: data for the third trial


 * 1) Syntax for use: perl -w hw6b.pl E1 E2 E3

$file1 = $ARGV[0]; $file2 = $ARGV[1]; $file3 = $ARGV[2];

read_file($file1, \@train1); read_file($file2, \@train2); read_file($file3, \@train3);
 * 1) Reading input files

$lambda_resting = (scalar(@train1) - 1)/$train1[-1]; $lambda_activated = (scalar(@train2) - 1)/$train2[-1];
 * 1) Calculating lambda, the rate of a spike (spikes/second)
 * 2) Use the quantity (scalar(@train) - 1) because the last entry is not a spike.
 * 3) Use the variable $train[-1] because the last entry is the total time for which the experiment ran.

$total_spikes = scalar(@train3) - 1; $total_time = $train3[-1];
 * 1) Calculating data for the third experiment
 * 2) Use the quantity (scalar(@train3) - 1) because the last entry is not a spike.
 * 3) Use the variable $train3[-1] because the last entry is the total time for which the experiment ran.

$y_x1 = &poisson($total_spikes, $lambda_activated, $total_time); # P(y|x = 1) $y_x0 = &poisson($total_spikes, $lambda_resting, $total_time); # P(y|x = 0)
 * 1) Calculating probabilities

$x1_y = (.01*$y_x1)/(.01*$y_x1 + .99*$y_x0);
 * 1)                          P(x = 1)*P(y|x = 1)
 * 2) P(x = 1|y) = ---
 * 3)               P(x = 1)*P(y|x = 1) + P(x = 0)*P(y|x = 0)
 * 4) where P(x = 1) = 0.01
 * 5)       P(x = 0 = 1 - 0.01 = 0.99

print "The posterior probability the neuron is activated during the third experiment is $x1_y\n";

sub poisson {   $num_events = shift(@_); #takes $total_spikes $lambda = shift(@_); #takes $lambda_activated $time = shift(@_); #takes $total_time $factorial = factorial($num_events);
 * 1) This subroutine calculates the probability of events occurring

$probability = (($lambda*$time)**$num_events) * exp(-1*$lambda*$time); $probability = $probability/$factorial; #     (L^k)*(e^-L) # P = -- #          k!    # Where L = lambda*time #      k = number of events return $probability; }

sub factorial {   $number = shift(@_); $answer = 1; for($i = 1; $i <= $number; $i++) {       $answer *= $i; }   return $answer; }
 * 1) This subroutines calculates the factorial of a number

sub read_file {   my $file = shift(@_); my $arrayRef = shift(@_); # Error handling- make sure the files are functional and exist die "File $file not found!\n" unless -e $file; open(INPUT, $file) or die "Please input a functional $file\n";
 * 1) This subroutine builds an array with the values in a file.

# Builds an array @train composed of the times recorded in each file @$arrayRef = ; chomp(@$arrayRef); }

close($file1); close($file2); close($file3);

=hw7.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program utilizes the UPGMA algorithm to output a phylogenetic tree in Newick format.
 * 2) Syntax for use: perl -w hw7.pl ATP6.dismat

open(FILE, $ARGV[0]) or die("Please input a functional file.\n");

$row = 0; $i = 0; while ($line = ) {   chomp $line; $input = $line; $col = 0; while ($input =~ /\s(0\S*)/g) {       $matrix[$row][$col] = $1; $col++; }   if ($line =~ /([\w]+)/g) {       $names[$i] = $1; $i++; $row++; } }
 * 1) Reading the distance matrix

$nodes_to_merge = $row; while ($nodes_to_merge > 2) {   # This code will only analyze the bottom half of the matrix as divided by the diagonal of 0 $row = 2; # Start at row 2 because row 0 contains 0 and row 1 will be set as the starting value $low = $matrix[1][0]; @low_coordinates = (1, 0); while ($row < $nodes_to_merge) {       for ($col = 0; $col < $row; $col++) {           if ($matrix[$row][$col] < $low) {               $low = $matrix[$row][$col]; @low_coordinates = ($row, $col); }       }        $row++; }
 * 1) Implementing UPGMA

# Removing the node names and creating a name for the new node $new_dist = $low/2; $new_name = "($names[$low_coordinates[0]]: $new_dist, $names[$low_coordinates[1]]: $new_dist)"; for ($i = $low_coordinates[1]; $i < @names - 1; $i++) {       if (defined($names[$i + 1])) {           $names[$i] = $names[$i + 1]; }   }    pop @names; pop @names; $names[@names] = $new_name;

# Determining the distances from the new node to existing nodes for ($i = 0; $i < $nodes_to_merge; $i++) {       if ($i != $low_coordinates[0]) {           if ($i != $low_coordinates[1]) {               $new_matrix[$i] = $matrix[$low_coordinates[0]][$i] + $matrix[$low_coordinates[1]][$i]; $new_matrix[$i] = $new_matrix[$i]/2; }       }    }    push(@new_matrix, 0);

# Removing the rows of the removed nodes for ($i = $low_coordinates[1]; $i < $nodes_to_merge; $i++) {       if (defined($matrix[$i + 1][0])) {           for ($j = 0; $j < $nodes_to_merge; $j++) {               $matrix[$i][$j] = $matrix[$i+1][$j]; }       }    }    for ($i = 0; $i < $nodes_to_merge; $i++) {       delete($matrix[$nodes_to_merge - 1][$i]); }   for ($i = 0; $i < $nodes_to_merge; $i++) {       delete($matrix[$nodes_to_merge - 2][$i]); }

# Removing the columns of the removed nodes for ($j = $low_coordinates[1]; $j < $nodes_to_merge; $j++) {       for ($i = 0; $i < @names - 1; $i++) {           $matrix[$i][$j] = $matrix[$i][$j+1]; }   }    for ($i = 0; $i < $nodes_to_merge; $i++) {       delete($matrix[$i][$nodes_to_merge - 1]); }

# Adding the distances for the new node to the final row for ($i = 0; $i < $nodes_to_merge - 1; $i++) {       $matrix[$nodes_to_merge - 2][$i] = $new_matrix[$i]; }

# Adding the distances for the new node to the final column for ($i = 0; $i < $nodes_to_merge - 1; $i++) {       $matrix[$i][$nodes_to_merge - 2] = $new_matrix[$i]; }   $nodes_to_merge--;

}

if ($nodes_to_merge == 2) {   $low = $matrix[1][0]; $new_dist = $low/2; $name = "($names[0]: $new_dist, $names[1]: $new_dist);\n"; }

print "$name\n";

close FILE;

ATP6.dismat
homo/1-226        0 0.0738015 0.0625875 0.0625875  0.138846  0.185344   0.27362   0.25005   0.25005  0.283744  0.294668  0.303424  0.292037   0.25005   0.25005  0.344928  0.334664  0.378076  0.341188   0.55351  0.440883  0.664133  0.705446   0.67512   0.77852  0.684134  0.761917 gorilla/1-226 0.0738015        0 0.0968085 0.0783045   0.13734   0.21582  0.292975   0.26982  0.269218  0.331603  0.346329  0.336812  0.323687  0.277403  0.283883  0.342698  0.324822  0.422704  0.377899  0.582386  0.452769  0.677851  0.744492  0.669119  0.782786  0.681806  0.767679 chimp/1-226 0.0625875 0.0968085        0 0.0410783  0.139793  0.186775   0.25005   0.25005   0.25005  0.293951  0.304112  0.303731  0.301102   0.25005  0.239496  0.342796    0.3348  0.402059  0.358068  0.548358  0.417685  0.647813   0.70631  0.677308   0.78007  0.692176  0.747405 bonono/1-226 0.0625875 0.0783045 0.0410783        0  0.102451  0.157578  0.264377   0.25005  0.234074  0.291372  0.294032  0.296683  0.285373  0.232634  0.227393  0.341954  0.332151  0.344768  0.324751  0.512022  0.395755  0.642641   0.69459  0.668376  0.745967  0.647237  0.728454 gibbon/1-226 0.138846   0.13734  0.139793  0.102451         0  0.185819   0.25005  0.291285  0.282233  0.301522   0.32659  0.335026  0.337893  0.284042  0.277295  0.341237  0.337737   0.40302  0.362032  0.534442  0.412029  0.632707  0.680976  0.682861  0.740661  0.620834   0.76428 pongo/1-226 0.185344   0.21582  0.186775  0.157578  0.185819         0  0.282629  0.340551  0.322734  0.370284  0.390024   0.37295  0.360348  0.343724  0.323095  0.367001  0.346484  0.368288    0.3372  0.539299  0.421431  0.601509   0.62374  0.626618  0.724853  0.680402  0.745465 macaque/1-226  0.27362  0.292975   0.25005  0.264377   0.25005  0.282629         0  0.395325  0.375857   0.45126  0.477512  0.457833  0.415587  0.410536  0.369051  0.423669  0.410064  0.513552   0.47531  0.595626  0.448917   0.69683  0.739737  0.816946   0.75921  0.800688  0.871841 sheep/1-226  0.25005   0.26982   0.25005   0.25005  0.291285  0.340551  0.395325         0  0.143201   0.25005   0.26092  0.218019  0.175309  0.149277  0.145046  0.264917  0.266418  0.320933   0.32057  0.510882  0.407572  0.638235  0.643307  0.670779  0.808468  0.701233  0.680042 hippo/1-226  0.25005  0.269218   0.25005  0.234074  0.282233  0.322734  0.375857  0.143201         0  0.220343  0.230662  0.153618  0.168653  0.152375  0.136345   0.25005  0.238274  0.304775  0.276335  0.534692  0.409851  0.668358  0.650715  0.678558  0.784349  0.674483  0.692608 mouse/1-226 0.283744  0.331603  0.293951  0.291372  0.301522  0.370284   0.45126   0.25005  0.220343         0   0.05566  0.202598   0.19302  0.196536  0.150912  0.290454  0.301954  0.370677  0.343738  0.579136  0.466294  0.647063  0.705644  0.666976  0.809968  0.680251   0.77888 rat/1-226 0.294668  0.346329  0.304112  0.294032   0.32659  0.390024  0.477512   0.26092  0.230662   0.05566         0  0.190183  0.223247  0.224384  0.211132  0.344846  0.362487  0.385427  0.385136  0.597733     0.495  0.679659  0.732182  0.711417  0.837308  0.709975  0.793782 nzbat/1-226 0.303424  0.336812  0.303731  0.296683  0.335026   0.37295  0.457833  0.218019  0.153618  0.202598  0.190183         0   0.25005  0.192644  0.190625  0.328974  0.302095  0.345166  0.301889  0.568169     0.495  0.653267  0.686655  0.710242  0.829091  0.711419  0.744206 jbat/1-226 0.292037  0.323687  0.301102  0.285373  0.337893  0.360348  0.415587  0.175309  0.168653   0.19302  0.223247   0.25005         0  0.154239  0.151568  0.279936  0.276858  0.372326  0.325599  0.512018  0.475725  0.657509  0.686936  0.700738  0.798023  0.744414  0.765137 dog/1-226  0.25005  0.277403   0.25005  0.232634  0.284042  0.343724  0.410536  0.149277  0.152375  0.196536  0.224384  0.192644  0.154239         0  0.125075  0.263376  0.261798  0.343473   0.34042   0.49293  0.459547  0.620788  0.651385  0.685095   0.79036  0.708747  0.732911 horse/1-226  0.25005  0.283883  0.239496  0.227393  0.277295  0.323095  0.369051  0.145046  0.136345  0.150912  0.211132  0.190625  0.151568  0.125075         0  0.204345  0.229419  0.339417  0.322686  0.481083  0.413802   0.63924  0.672514  0.687837  0.767414  0.699996   0.68591 bwhale/1-226 0.344928  0.342698  0.342796  0.341954  0.341237  0.367001  0.423669  0.264917   0.25005  0.290454  0.344846  0.328974  0.279936  0.263376  0.204345         0 0.0553572  0.431796  0.382287  0.538051  0.473931  0.706893  0.724528  0.654348  0.754982  0.748748  0.762417 fwhale/1-226 0.334664  0.324822    0.3348  0.332151  0.337737  0.346484  0.410064  0.266418  0.238274  0.301954  0.362487  0.302095  0.276858  0.261798  0.229419 0.0553572         0  0.438977   0.37765  0.558118  0.474928  0.718646  0.709553   0.66775  0.775631  0.774051  0.790799 marsupial/1-226 0.378076  0.422704  0.402059  0.344768   0.40302  0.368288  0.513552  0.320933  0.304775  0.370677  0.385427  0.345166  0.372326  0.343473  0.339417  0.431796  0.438977         0  0.102451  0.595315  0.527027  0.613888  0.591789   0.60835  0.715894  0.650122   0.68022 opossum/1-226 0.341188  0.377899  0.358068  0.324751  0.362032    0.3372   0.47531   0.32057  0.276335  0.343738  0.385136  0.301889  0.325599   0.34042  0.322686  0.382287   0.37765  0.102451         0  0.582164  0.526104  0.619715   0.56556  0.616542   0.70414  0.645983  0.661645 elephant/1-230  0.55351  0.582386  0.548358  0.512022  0.534442  0.539299  0.595626  0.510882  0.534692  0.579136  0.597733  0.568169  0.512018   0.49293  0.481083  0.538051  0.558118  0.595315  0.582164         0  0.652143  0.722695  0.802194  0.782908  0.838024  0.871059  0.975715 monkey/1-226 0.440883  0.452769  0.417685  0.395755  0.412029  0.421431  0.448917  0.407572  0.409851  0.466294     0.495     0.495  0.475725  0.459547  0.413802  0.473931  0.474928  0.527027  0.526104  0.652143         0  0.766025  0.777113  0.737799  0.820082  0.709234  0.755398 falcon/1-227 0.664133  0.677851  0.647813  0.642641  0.632707  0.601509   0.69683  0.638235  0.668358  0.647063  0.679659  0.653267  0.657509  0.620788   0.63924  0.706893  0.718646  0.613888  0.619715  0.722695  0.766025         0  0.157714  0.371255       0.5  0.431898  0.464752 buteo/1-227 0.705446  0.744492   0.70631   0.69459  0.680976   0.62374  0.739737  0.643307  0.650715  0.705644  0.732182  0.686655  0.686936  0.651385  0.672514  0.724528  0.709553  0.591789   0.56556  0.802194  0.777113  0.157714         0  0.326191  0.481998  0.395332  0.436809 turtle/1-227  0.67512  0.669119  0.677308  0.668376  0.682861  0.626618  0.816946  0.670779  0.678558  0.666976  0.711417  0.710242  0.700738  0.685095  0.687837  0.654348   0.66775   0.60835  0.616542  0.782908  0.737799  0.371255  0.326191         0  0.470511  0.451024  0.535346 alligator/1-225  0.77852  0.782786   0.78007  0.745967  0.740661  0.724853   0.75921  0.808468  0.784349  0.809968  0.837308  0.829091  0.798023   0.79036  0.767414  0.754982  0.775631  0.715894   0.70414  0.838024  0.820082       0.5  0.481998  0.470511         0  0.616843  0.643341 salamander/1-227 0.684134  0.681806  0.692176  0.647237  0.620834  0.680402  0.800688  0.701233  0.674483  0.680251  0.709975  0.711419  0.744414  0.708747  0.699996  0.748748  0.774051  0.650122  0.645983  0.871059  0.709234  0.431898  0.395332  0.451024  0.616843         0  0.438723 frog/1-226 0.761917  0.767679  0.747405  0.728454   0.76428  0.745465  0.871841  0.680042  0.692608   0.77888  0.793782  0.744206  0.765137  0.732911   0.68591  0.762417  0.790799   0.68022  0.661645  0.975715  0.755398  0.464752  0.436809  0.535346  0.643341  0.438723         0

ATP6_functional.dismat
homo/1-226        0 0.0738015 0.0625875 0.0625875 gorilla/1-226 0.0738015        0 0.0968085 0.0783045 chimp/1-226 0.0625875 0.0968085        0 0.0410783 bonono/1-226 0.0625875 0.0783045 0.0410783        0

=hw8_simple.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program simulates random growth.  Initially, a point will be fixed at the center.  Then, points will be released from one of the four edges chosen at random and move until it bumps into a fixed pixel.  Pixels will be colored randomly.
 * 2) Syntax for use: perl -w hw8_simple.pl size output_file_name

use GD::Image; $size = $ARGV[0]; $image = new GD::Image($size, $size); $whiteIndex = $image->colorAllocate(255, 255, 255);

for ($i = 0; $i < 256; $i++) {   $color1 = int(rand(256)); $color2 = int(rand(256)); $color3 = int(rand(256)); if (!(($color1 == 255) & ($color2 == 255) & ($color3 == 255))) {       $colors[$i] = $image->colorAllocate($color1, $color2, $color3); } }
 * 1) Making a matrix of random colors

if ($size%2 == 0) {   $image->setPixel($size/2, $size/2, $colors[int(rand(256))]); } else {   $image->setPixel(($size - 1)/2, ($size - 1)/2, $colors[int(rand(256))]); }
 * 1) Setting the center pixel

$final_flag = 0; # $final_flag = 1 when a pixel becomes fixed on an edge while ($final_flag == 0) {   @coordinates = place_on_edge; $stick_flag = 0; # $stick_flag = 1 when a pixel becomes fixed while ($stick_flag == 0) {       # Moving randomly $move = int(rand(4)); # Determines in what direction the pixel will attempt to move if ($move == 0) # Move upwards {           if ($coordinates[1] != 0) {               $pixelColorIndex = $image->getPixel($coordinates[0], $coordinates[1] - 1); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0], $coordinates[1] - 1); }               else {                   $image->setPixel($coordinates[0], $coordinates[1], $colors[int(rand(256))]); $stick_flag = 1; }           }        }        elsif ($move == 1) # Move to left {           if ($coordinates[0] != 0) {               $pixelColorIndex = $image->getPixel($coordinates[0] - 1, $coordinates[1]); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0] - 1, $coordinates[1]); }               else {                   $image->setPixel($coordinates[0], $coordinates[1], $colors[int(rand(256))]); $stick_flag = 1; }           }        }        elsif ($move == 2) # Move to right {           if ($coordinates[0] != ($size - 1)) {               $pixelColorIndex = $image->getPixel($coordinates[0] + 1, $coordinates[1]); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0] + 1, $coordinates[1]); }               else {                   $image->setPixel($coordinates[0], $coordinates[1], $colors[int(rand(256))]); $stick_flag = 1; }           }        }        elsif ($move == 3) # Move downwards {           if ($coordinates[1] != ($size - 1)) {               $pixelColorIndex = $image->getPixel($coordinates[0], $coordinates[1] + 1); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0], $coordinates[1] + 1); }               else {                   $image->setPixel($coordinates[0], $coordinates[1], $colors[int(rand(256))]); $stick_flag = 1; }           }        }    }

# Determine if the stuck pixel is at the edge, thus terminating this while loop if ($coordinates[1] == 0) {       $final_flag = 1; }   elsif ($coordinates[0] == 0) {       $final_flag = 1; }   elsif ($coordinates[0] == ($size - 1)) {       $final_flag = 1; }   elsif ($coordinates[1] == ($size - 1)) {       $final_flag = 1; } }

$output = $ARGV[1]; open(OUTPUT, '>', $output); binmode(OUTPUT); print OUTPUT $image->png(0);
 * 1) Creating the output image

sub place_on_edge {   #Placing a pixel at a random point at the edge $rand = int(rand($size)); $case = int(rand(4)); if ($case == 0) #Along the top row {       return @coordinates = ($rand, 0); }   elsif($case == 1) #Along the left column {       return @coordinates = (0, $rand); }   elsif($case == 2) #Along the right column {       return @coordinates = ($size - 1, $rand); }   elsif($case == 3) #Along the bottom row {       return @coordinates = ($rand, $size - 1); } }
 * 1) This subroutine uses the rand function to randomly select an edge to begin a pixel's movement

=hw8_full.pl=
 * 1) ! /usr/bin/perl -w


 * 1) This program simulates random growth.  Initially, a point will be fixed at the center.  Then, points will be released from one of the four edges chosen at random and move until it bumps into a fixed pixel.  Pixels will be colored based on their distance from the center.
 * 2) Syntax for use: perl -w hw8_full.pl size output_file_name

use GD::Image; $size = $ARGV[0]; $image = new GD::Image($size, $size); $whiteIndex = $image->colorAllocate(255, 255, 255);

$i = 0; while ($i < 11) {   $color1 = int(rand(256)); $color2 = int(rand(256)); $color3 = int(rand(256)); if (!(($color1 == 255) & ($color2 == 255) & ($color3 == 255))) {       $colors[$i] = $image->colorAllocate($color1, $color2, $color3); $i++; } }
 * 1) Making a matrix of random colors other than white

if ($size%2 == 0) {   @center = ($size/2, $size/2); $image->setPixel($size/2, $size/2, $colors[int(rand(10))]); } else {   @center = (($size - 1)/2, ($size - 1)/2); $image->setPixel(($size - 1)/2, ($size - 1)/2, $colors[int(rand(10))]); }
 * 1) Setting the center pixel

@coordinates = (0, 0); $max_distance = calc_distance; if ($max_distance%10 != 0) {   # $region_dist = the distance between regions $region_dist = $max_distance - ($max_distance%10); $region_dist = $region_dist/10; } else {   # $region_dist = the distance between regions $region_dist = $max_distance/10; }
 * 1) Determining different regions of distance from the center

$final_flag = 0; while ($final_flag == 0) {   @coordinates = place_on_edge; $stick_flag = 0; # $stick_flag = 1 when a pixel becomes fixed while ($stick_flag == 0) {       #Moving randomly $move = int(rand(4)); # Determines in what direction the pixel will attempt to move if ($move == 0) # Move upwards {           if ($coordinates[1] != 0) {               $pixelColorIndex = $image->getPixel($coordinates[0], $coordinates[1] - 1); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0], $coordinates[1] - 1); }               else {                   $distance = calc_distance; $region = find_region; $image->setPixel($coordinates[0], $coordinates[1], $colors[$region]); $stick_flag = 1; }           }        }        elsif ($move == 1) # Move to left {           if ($coordinates[0] != 0) {               $pixelColorIndex = $image->getPixel($coordinates[0] - 1, $coordinates[1]); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0] - 1, $coordinates[1]); }               else {                   $distance = calc_distance; $region = find_region; $image->setPixel($coordinates[0], $coordinates[1], $colors[$region]); $stick_flag = 1; }           }        }        elsif ($move == 2) # Move to right {           if ($coordinates[0] != ($size - 1)) {               $pixelColorIndex = $image->getPixel($coordinates[0] + 1, $coordinates[1]); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0] + 1, $coordinates[1]); }               else {                   $distance = calc_distance; $region = find_region; $image->setPixel($coordinates[0], $coordinates[1], $colors[$region]); $stick_flag = 1; }           }        }        elsif ($move == 3) # Move downwards {           if ($coordinates[1] != ($size - 1)) {               $pixelColorIndex = $image->getPixel($coordinates[0], $coordinates[1] + 1); if ($pixelColorIndex == $whiteIndex) {                   @coordinates = ($coordinates[0], $coordinates[1] + 1); }               else {                   $distance = calc_distance; $region = find_region; $image->setPixel($coordinates[0], $coordinates[1], $colors[$region]); $stick_flag = 1; }           }        }    }
 * 1) Setting pixels and their movement

# Determine if the stuck pixel is at the edge, thus terminating this while loop if ($coordinates[1] == 0) {       $final_flag = 1; }   elsif ($coordinates[0] == 0) {       $final_flag = 1; }   elsif ($coordinates[0] == ($size - 1)) {       $final_flag = 1; }   elsif ($coordinates[1] == ($size - 1)) {       $final_flag = 1; } }

$output = $ARGV[1]; open(OUTPUT, '>', $output); binmode(OUTPUT); print OUTPUT $image->png(0);
 * 1) Creating the output image

sub place_on_edge {   #Placing a pixel at a random point at the edge $rand = int(rand($size)); $case = int(rand(4)); if ($case == 0) #Along the top row {       return @coordinates = ($rand, 0); }   elsif($case == 1) #Along the left column {       return @coordinates = (0, $rand); }   elsif($case == 2) #Along the right column {       return @coordinates = ($size - 1, $rand); }   elsif($case == 3) #Along the bottom row {       return @coordinates = ($rand, $size - 1); } }
 * 1) This subroutine uses the rand function to randomly select an edge to begin a pixel's movement

sub calc_distance {   $distance1 = $coordinates[0] - $center[0]; $distance1 **= 2; $distance2 = $coordinates[1] - $center[1]; $distance2 **= 2; $distance = $distance1 + $distance2; $distance = sqrt($distance); return $distance; }
 * 1) This subroutine calculates distance between two points using the formula:
 * 2) dist = sqrt[(x2 - x1)^2 + (y2 - y1)^2]

sub find_region {   if ((0 < $distance) & ($distance <= $region_dist)) {       $region = 0; }   if (($region_dist < $distance) & ($distance <= 2*$region_dist)) {       $region = 1; }   if ((2*$region_dist < $distance) & ($distance <= 3*$region_dist)) {       $region = 2; }   if ((3*$region_dist < $distance) & ($distance <= 4*$region_dist)) {       $region = 3; }   if ((4*$region_dist < $distance) & ($distance <= 5*$region_dist)) {       $region = 4; }   if ((5*$region_dist < $distance) & ($distance <= 6*$region_dist)) {       $region = 5; }   if ((6*$region_dist < $distance) & ($distance <= 7*$region_dist)) {       $region = 6; }   if ((7*$region_dist < $distance) & ($distance <= 8*$region_dist)) {       $region = 7; }   if ((8*$region_dist < $distance) & ($distance <= 9*$region_dist)) {       $region = 8; }   if ((9*$region_dist < $distance) & ($distance <= 10*$region_dist)) {       $region = 9; }   if ((10*$region_dist < $distance) & ($distance <= $max_distance)) {       $region = 10; }   return $region; }
 * 1) This subroutine uses a pixel's distance from the center to determine in which region a pixel stopped; the region number is used to access the matrix of colors