#!/usr/bin/perl -w
use strict;
#  A perl script to map bisulfite sequencing reads from Illumina GA to padlock captured targets
# USAGE: mapBisReads.pl raw_read_file target_CG_file target_TG_File > sample_mapped.txt
#   raw_read_file: raw sequencing reads, one line per read, no quality scores
#   target_CG_file: the bisulfite converted sequence of all captured targets assuming complete methylation
#   target_TG_file: the bisulfite converted sequence of all captured targets assuming complete unmethylation
#   
# Dependencies:
#   (1) Bioperl (available from CPAN)
#   (2) SOAP (available from soap.genomics.org.cn)
#
#
# Written by Kun Zhang (kzhang@bioeng.ucsd.edu), Aug. 2008
#

# THE FOLLOWING TWO LINES NEED TO MODIFIED BASED THE LOCATIONS OF BIOPERL AND SOAP.
use lib qw(/home/kunzhang/PerlLib);
my $SOAP_EXE = "/home/kunzhang/softwares/soap_1.11/soap.short";

use Bio::SeqIO;
my $setSize = 400000;

my %mappedReadInfo;
my %CpgInfo;
my %cInfo;
my %conversionStat;
my %targetLen;

my $readFile = $ARGV[0];
my $targetCgFile = $ARGV[1];
my $targetTgFile = $ARGV[2];

my %rcTable;
$rcTable{'A'}='T';
$rcTable{'T'}='A';
$rcTable{'G'}='C';
$rcTable{'C'}='G';
$rcTable{'N'}='N';

my $tmpCgFaFile = $readFile . ".CG.subset.fa";
my $tmpTgFaFile = $readFile . ".TG.subset.fa";

&loadTargetCpGs($targetCgFile);

open(READ_FILE, "$readFile") || die("Error in opening file $readFile.\n");
open(CG_FA_FILE, ">$tmpCgFaFile") || die("Error in opening file $tmpCgFaFile.\n");
open(TG_FA_FILE, ">$tmpTgFaFile") || die("Error in opening file $tmpTgFaFile.\n");

my ($totalReads, $mappedReads);
my $read_number = 1;
while(my $long_line = <READ_FILE>){
        next if($long_line =~ /[N.]/);
        chop($long_line);
        my $seqID = 'c'.$read_number;
        my @words = split(/\t/, $long_line);
        my $line = ($words[4] && $words[4] =~ /[ATGC]/) ? $words[4] : $words[0];
	$line = &revComp($line) if(&guess_strand($line) eq'rev');
	$mappedReadInfo{$seqID}->{'oS'}=$line;
	print CG_FA_FILE ">$seqID\n$line\n";
	$line =~ s/C/T/gi;
	print TG_FA_FILE ">$seqID\n$line\n";
	$read_number ++;
	if($read_number >= $setSize){
		$totalReads+=$read_number;
		close(CG_FA_FILE);
		close(TG_FA_FILE);
		$read_number = 1;
		&mapOneSet($tmpCgFaFile, $tmpTgFaFile,  $targetCgFile, $targetTgFile);		
		&reportMapping();
		open(CG_FA_FILE, ">$tmpCgFaFile") || die("Error in opening file $tmpCgFaFile.\n");
		open(TG_FA_FILE, ">$tmpTgFaFile") || die("Error in opening file $tmpTgFaFile.\n");
		undef(%mappedReadInfo);
	}
}
close(READ_FILE);

if($read_number > 0){
		close(CG_FA_FILE);
		close(TG_FA_FILE);
		$totalReads+=$read_number;
		$read_number = 1;
		&mapOneSet($tmpCgFaFile, $tmpTgFaFile,  $targetCgFile, $targetTgFile);		
		&reportMapping();
}
unlink("soap.log");
unlink($tmpCgFaFile);
unlink($tmpTgFaFile);
print "$mappedReads of $totalReads were mapped.\n";

sub mapOneSet(){
	my $readCGFile = shift;
	my $readTGFile = shift;
	my $targetCGFile = shift;
	my $targetTGFile = shift;
	
	#map reads to unmethylated targets
	my $TGsoapOutFile = $readTGFile . ".soap.out";	
	my $TG_cmd = "$SOAP_EXE -a $readTGFile -d $targetTGFile -o $TGsoapOutFile -v 1 -e 10000 -c 0 -f 8 -r 0 -s 12 -n 1 > $readFile." . "soap.log";
	system($TG_cmd);
	open(INFILE, "$TGsoapOutFile") || die("Error in opening file $TGsoapOutFile!\n");
	while(my $line = <INFILE>){
		#parse SOAP output
		chop($line);
		my @fields = split(/\t/, $line);
		my $seqID = $fields[0];
		my $querySeq = $fields[1];
		my $targetID = $fields[7];
		my $strand = $fields[6];
		my $start = $fields[8];
		my $nMismatch = $fields[9];
		#no need to update the mapping info if the mapping results are consistent between two type of templates.
		next if($mappedReadInfo{$seqID}->{"tID"} && $targetID eq $mappedReadInfo{$seqID}->{"tID"} &&
				$mappedReadInfo{$seqID}->{"start"} && $mappedReadInfo{$seqID}->{"start"} == $start);
		#choose the one that has fewer mismatches.
		next if($mappedReadInfo{$seqID}->{'nMM'} && $mappedReadInfo{$seqID}->{'nMM'} <$nMismatch);
		$mappedReadInfo{$seqID}->{"qS"} = $querySeq;
		$mappedReadInfo{$seqID}->{"tID"} = $targetID;
		$mappedReadInfo{$seqID}->{"start"} = $start;
		$mappedReadInfo{$seqID}->{"strand"} = $strand;
	}
	close(INFILE);
	unlink($TGsoapOutFile);
	return %mappedReadInfo;
}

sub reportMapping(){
	foreach my $seqID (keys(%mappedReadInfo)){
		next if(!$mappedReadInfo{$seqID}->{"tID"});
		$mappedReads ++;
		my $querySeq =  ($mappedReadInfo{$seqID}->{"strand"} eq '+') ? 
						$mappedReadInfo{$seqID}->{"oS"} :
						&revComp($mappedReadInfo{$seqID}->{"oS"});
		my $qLen = length($querySeq);
		my $tID = $mappedReadInfo{$seqID}->{"tID"};
		my ($chr, $range, $offset, $strand) = split(/:/, $tID);
		my ($start, $end) = split(/\-/, $range);
		my $targetOffset;;
		my $id = $chr.":".$start."-".$end;
		#check whether the mapped sequence contains any CpG site or not.
		my $containCpG = 0;
		
		for(my $i=0; $i<$qLen; $i++){
			if($strand =~ /W/i){
				$targetOffset = $offset +  $mappedReadInfo{$seqID}->{"start"} + $i -1;
			}else{
				$targetOffset = $offset + $targetLen{$tID} - $mappedReadInfo{$seqID}->{"start"} - $i -1 ; 
			}
			
			if($CpgInfo{$id}->{$targetOffset}){
				$containCpG=1; 
			}			
		}		
		
		if ($containCpG){				
			print "$querySeq\t$id\t$strand\t";
			if($strand =~ /W/i){
				print $mappedReadInfo{$seqID}->{"start"}+$offset;
			}else{
				print $offset + $targetLen{$tID} - $mappedReadInfo{$seqID}->{"start"};
			}
				
			for(my $i=0; $i<$qLen; $i++){
					if($strand =~ /W/i){
						$targetOffset = $offset +  $mappedReadInfo{$seqID}->{"start"} + $i -1;
					}else{
						$targetOffset = $offset + $targetLen{$tID} - $mappedReadInfo{$seqID}->{"start"} - $i -1 ; 
					}				
					if($CpgInfo{$id}->{$targetOffset}){
						print "\t$targetOffset\t", substr($querySeq,$i,1); 
					}		
			}
			print "\n";
		}
	}
}


sub loadTargetCpGs(){
	my $targetSeqFile = shift;
	my $seq_in = Bio::SeqIO->new(-file => $targetSeqFile , -format =>'Fasta');	
	while (my $seqObj = $seq_in->next_seq()){
		my $tSeq = $seqObj->seq();
		my $tID = $seqObj->id();
		my ($chr, $range, $offset, $strand) = split(/:/, $tID);
		$chr =~ s/>//;
		my ($start, $end) = split(/\-/, $range);
		my $tRange = length($tSeq)-1;
		for(my $i=0; $i<$tRange; $i++){
			if(substr($tSeq,$i,2) eq "CG"){
				my $cpgOffset = $strand eq "W" ? 
				        $offset + $i : 
				        $offset + $tRange - $i - 1;
				$CpgInfo{$chr.":".$start."-".$end}->{$cpgOffset}=$strand;
			}			
		}
		$targetLen{$tID}=$seqObj->length();
	}
}

sub revComp(){
	my $seq = shift;
	my $rcSeq='';
	for(my $i=0; $i<length($seq); $i++){
		$rcSeq = $rcTable{substr($seq,$i,1)} . $rcSeq;
	}
	return $rcSeq;
}

sub guess_strand(){
	my $seq = shift;
	my %baseCounts;
	$baseCounts{'A'}=0.001;
	$baseCounts{'T'}=0.001;
	$baseCounts{'G'}=0.001;
	$baseCounts{'C'}=0.001;
	while(my $base = chop($seq)){
		$baseCounts{$base}++;
	}
	if($baseCounts{'T'}/$baseCounts{'C'} >$baseCounts{'A'}/$baseCounts{'G'}) {
		return "fwd";
	}else{
		return "rev";
	}
}
