##perl routines for process_vwf 

#############################
sub complement {
    $_[0] =~ y/CGAT/GCTA/;
    return $_[0];
}

###############################
sub translate_seq_to_protein {
    my($string) = @_;
    my %DNA_code = (
	'GCT' => 'A', 'GCC' => 'A', 'GCA' => 'A', 'GCG' => 'A', 'TTA' => 'L',
	'TTG' => 'L', 'CTT' => 'L', 'CTC' => 'L', 'CTA' => 'L', 'CTG' => 'L',
	'CGT' => 'R', 'CGC' => 'R', 'CGA' => 'R', 'CGG' => 'R', 'AGA' => 'R',
        'AGG' => 'R', 'AAA' => 'K', 'AAG' => 'K', 'AAT' => 'N', 'AAC' => 'N',
        'ATG' => 'M', 'GAT' => 'D', 'GAC' => 'D', 'TTT' => 'F', 'TTC' => 'F',
        'TGT' => 'C', 'TGC' => 'C', 'CCT' => 'P', 'CCC' => 'P', 'CCA' => 'P',
        'CCG' => 'P', 'CAA' => 'Q', 'CAG' => 'Q', 'TCT' => 'S', 'TCC' => 'S',
        'TCA' => 'S', 'TCG' => 'S', 'AGT' => 'S', 'AGC' => 'S', 'GAA' => 'E',
        'GAG' => 'E', 'ACT' => 'T', 'ACC' => 'T', 'ACA' => 'T', 'ACG' => 'T',
        'GGT' => 'G', 'GGC' => 'G', 'GGA' => 'G', 'GGG' => 'G', 'TGG' => 'W',
        'CAT' => 'H', 'CAC' => 'H', 'TAT' => 'Y', 'TAC' => 'Y', 'ATT' => 'I',
        'ATC' => 'I', 'ATA' => 'I', 'GTT' => 'V', 'GTC' => 'V', 'GTA' => 'V',
        'GTG' => 'V', 'TAA' => 'X', 'TAG' => 'B', 'TGA' => 'X');

#    my %DNA_code = (
#	'GCT' => 'A', 'GCC' => 'A', 'GCA' => 'A', 'GCG' => 'A', 'TTA' => 'L',
#	'TTG' => 'L', 'CTT' => 'L', 'CTC' => 'L', 'CTA' => 'L', 'CTG' => 'L',
#	'CGT' => 'R', 'CGC' => 'R', 'CGA' => 'R', 'CGG' => 'R', 'AGA' => 'R',
#        'AGG' => 'R', 'AAA' => 'K', 'AAG' => 'K', 'AAT' => 'N', 'AAC' => 'N',
#        'ATG' => 'M', 'GAT' => 'D', 'GAC' => 'D', 'TTT' => 'F', 'TTC' => 'F',
#        'TGT' => 'C', 'TGC' => 'C', 'CCT' => 'P', 'CCC' => 'P', 'CCA' => 'P',
#        'CCG' => 'P', 'CAA' => 'Q', 'CAG' => 'Q', 'TCT' => 'S', 'TCC' => 'S',
#        'TCA' => 'S', 'TCG' => 'S', 'AGT' => 'S', 'AGC' => 'S', 'GAA' => 'E',
#        'GAG' => 'E', 'ACT' => 'T', 'ACC' => 'T', 'ACA' => 'T', 'ACG' => 'T',
#        'GGT' => 'G', 'GGC' => 'G', 'GGA' => 'G', 'GGG' => 'G', 'TGG' => 'W',
#        'CAT' => 'H', 'CAC' => 'H', 'TAT' => 'Y', 'TAC' => 'Y', 'ATT' => 'I',
#        'ATC' => 'I', 'ATA' => 'I', 'GTT' => 'V', 'GTC' => 'V', 'GTA' => 'V',
#        'GTG' => 'V', 'TAA' => 'X', 'TAG' => 'X', 'TGA' => 'X');
    my($protein) = '';
    my($end_pos) = int(length($string)/3) * 3;
    for ($tt=0;$tt<$end_pos;$tt+=3) {
	$p4 = uc(substr($string,$tt,3));
	if (length($p4) < 2) {
	    print "P4 p1 " . $p1 . " p2 " . $p2 . " p3 " . $p3 . "\n";
	}
	if ($tmp = $DNA_code{$p4}) {
	    $protein .= $tmp;
	} else {
	    $protein .= "*";
	}
    }

    return $protein;
}

##########################################
sub quality_cutoff {
    my($string) = @_;
    my($last_pos) = 0;
    my($cutoff) = 20;
#    print "q string " . $string . "\n";
    for ($tt=0;$tt<=length($string);$tt++) {
        $tmp = substr($string,$tt,1);
        $q_value = ord($tmp) - 33 ;
#           print "q value " . $q_value . "," . " chr " . $tmp . " last " . $last_pos. "\n" ;                      
        if ($q_value > $cutoff) {
            $last_pos = $tt;
        } else {
            last;
        }
    }
    $last_pos ++;
    return $last_pos;
}

#######################################
sub reverse_compliment {
    my($string) = @_;
    my($string_comp) = '';
    $string = reverse($string);
    for ($tt=0;$tt<=length($string);$tt++) {
        $string_comp .= complement(substr($string,$tt,1));
    }
    return $string_comp;
}

#######################################
sub string_compliment {
    my($string) = @_;
    my($string_comp) = '';
    for ($tt=0;$tt<=length($string);$tt++) {
        $string_comp .= complement(substr($string,$tt,1));
    }
    return $string_comp;
}

#####################################
sub compare_reads {
    my($pos1_start,$pos1_end,$pos2_start,$pos2_end,$seq1,$seq2) = @_;
    my($mis) = 0;
    my($r1_start,$r1_end,$r2_start);
    my($tt);
    my($r1_string,$r2_string) = "";
    if ($pos1_start == $pos2_start) {
	$r1_start = 0;
	$r2_start = 0;
    } elsif ($pos1_start > $pos2_start) {
	$r1_start = 0;
	$r2_start = $pos1_start - $pos2_start;
    } else {
	$r1_start = $pos2_start - $pos1_start;
	$r2_start = 0;
    }
    $r1_length = $pos1_end - $pos1_start;
    $r2_length = $pos2_end - $pos2_start;
    $r1_end = min($r1_length,$r2_length);
    for ($tt=$r1_start;$tt<=$r1_end;$tt++) {
	if (!substr($seq1,$tt,1) || !substr($seq2,$r2_start,1) ) {
	    last;
	} else {
	    $r1_string .= substr($seq1,$tt,1);
	    $r2_string .= substr($seq2,$r2_start,1);
	}
	if (substr($seq1,$tt,1) ne substr($seq2,$r2_start,1)) {
	    $mis ++;
	}
	$r2_start ++;
    }
    if ($display) {
	print "\nIn compare_read\n";
	print "pos1_start " . $pos1_start . "\n";
	print "pos1_end " . $pos1_end . "\n";
	print "pos2_start " . $pos2_start . "\n";
	print "pos2_end " . $pos2_end . "\n";
	print "r1_start " . $r1_start . "\n";
	print "r1_end " . $r1_end . "\n";
	print "r2_start " . $r2_start . "\n";
	print "seq1 " . $seq1 . "\n";
	print "seq2 " . $seq2 . "\n";
	print "mis " . $mis . "\n";
	print "r1_string " . $r1_string . "\n";
	print "r2_string " . $r2_string . "\n";
    }
    return ($mis,$r1_string,$r2_string);
}

##################################################### 
## average      
#####################################################                                                                                                                                                      
sub average {
    my(@fields) = @_;
    my $sum;
    my $sumsq;
    my $n = 0;
    for $tt (0 .. $#fields){
	$sum += $fields[$tt];
        $sumsq += $fields[$tt]**2;
        $n++;
    }
    if ($n) {
        return $sum/$n;
    } else {
        return 0.0;
    }
}

#########################################                                                                                                                                             
sub deviation {
    my(@array) = @_;
    if ($#array == 0) {
        return 0.0;
    } elsif ($#array > 0) {
        undef(@variance);
        $avg = average(@array);
        for $tt (0 .. $#array) {
            $tmp = ($array[$tt] - $avg)**2;
            push(@variance,$tmp);
        }
        $sum = 0;
        for $tt (0 .. $#variance) {
            $sum += $variance[$tt];
        }
        ##$deviation = sqrt($sum/$#assay);
	$deviation = sqrt(average(@variance));
        return $deviation;
    } else {
        return 0.0;
    }
}



1;
