#!/usr/bin/perl -w
# A perl script for assemblying capturing arms for padlock probes into Agilent oligonucleotides
# USAGE: ./probeAssembler.pl ProbeSetName_H1H2_list.txt > ProbeSetName_AgilentOligos.txt
#   ProbeSetName_H1H2_list.txt: the output of getCpgProbes.pl
# Written by Kun Zhang (kzhang@bioeng.ucsd.edu), Aug. 2008
#

use strict;
my %cmpTable = ("A", "T", "T", "A", "G", "C", "C", "G","R","Y","Y", "R", "M", "K", "K", "M", "S", "S", "W", "W",  "N", "N");
my $AP1 = "GTAGACTGGAAGAGCACTGTT"; #V4
my $AP2 = "GATCGGATACGCATGAGGCTA"; #V8
my $linker = "GTTGGAGGCTCATCGTTCCTATTCAGCTGCAGATGTTATCGAGGTCCGAC"; #V8
my %seqDB;
my $maxLen = 58;

my @bases=("A", "T", "C", "G");
open(INFILE,$ARGV[0]) || die("Error in opening input file!\n");
my $seqFile = $ARGV[0] . ".seq";
open(OUTFILE, ">$seqFile") || die("Error in opening ouput file!\n");
while(my $line = <INFILE>){
	next if($line !~ /^chr/);
	chop($line);
	my @fields = split(/[\t ]+/, $line);
	my $LP = $fields[3];
	$LP =~ s/C/Y/i;
	my $RP = $fields[5];
	$RP =~ s/C/Y/i;
	my $Cpg = $fields[0];
	my $target_start = $fields[1]-50;
	my $target_end = $fields[2]-50;
	my $len = length($LP.$RP);

	next if($len > $maxLen);
	my $fill_in="";
	while($len < $maxLen){
		$fill_in = $fill_in . $bases[int(rand(2))];
		$len ++; 
	}
	my $padlock = $AP1. &RevComp($LP) . $fill_in . $linker . &RevComp($RP) . $AP2;
	my @fragments = split(/R/, $padlock);
	my @variants;
	push(@variants, $fragments[0]);
	for(my $i=0; $i<scalar(@fragments)-1; $i++){
		my $size = scalar(@variants);
		for(my $j=0; $j < $size; $j++){
			push(@variants, $variants[$j] . 'G' . $fragments[$i+1]);
		}
		for(my $j=0; $j < $size; $j++){
			$variants[$j] = $variants[$j] . 'A' . $fragments[$i+1];
		}
	}
    for(my $i=1; $i<=scalar(@variants); $i++){
	    print OUTFILE "$Cpg:$target_start-$target_end:$i\t", $variants[$i-1], "\n";
    }
}
close(INFILE);
close(OUTFILE);

#-------------------
sub RevComp(){
#-------------------
    my $seq = shift;
	my $seqLen = length($seq);
	my $revcom = $cmpTable{substr($seq,0,1)};
	for(my $i = 1; $i < $seqLen; $i++){
		$revcom = $cmpTable{substr($seq, $i, 1)} . $revcom;
	}
    return $revcom;
}
