/* @source newcpgseek application
**
** Reports ALL cpg rich regions in a sequence
** @author Copyright (C) Alan Bleasby (ableasby@hgmp.mrc.ac.uk)
** @author Modified by Rodrigo Lopez (rls@ebi.ac.uk)
** @@
**
** Original program "CPGREPORT" by Rodrigo Lopez (EGCG 1995)
**  CpG island finder. Larsen et al Genomics 13 1992 p1095-1107
**  "usually defined as >200bp with %GC > 50% and obs/exp CpG >
**  0.6". Here use running sum rather than window to score: if not CpG
**  at position i, then decrement runSum counter, but if CpG then runSum
**  += CPGSCORE.     Spans > threshold are searched
**  for recursively.
**
** This program is free software; you can redistribute it and/or
** modify it under the terms of the GNU General Public License
** as published by the Free Software Foundation; either version 2
** of the License, or (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
******************************************************************************/

#include "emboss.h"
#include <math.h>
#include <stdlib.h>




static void newcpgseek_cpgsearch(AjPFile *outf, ajint s, ajint len,
				 const char *seq,
				 const char *name, ajint begin, ajint score);
static void newcpgseek_calcgc(ajint from, ajint to, const char *p, ajint *dcg,
			      ajint *dgc, ajint *gc);




/* @prog newcpgseek ***********************************************************
**
** Reports CpG rich regions
**
******************************************************************************/

int main(int argc, char **argv)
{
    AjPSeqall seqall;
    AjPSeq seq    = NULL;
    AjPFile outf  = NULL;
    AjPStr strand = NULL;
    AjPStr substr = NULL;

    ajint begin;
    ajint end;
    ajint len;
    ajint score;

    embInit("newcpgseek",argc,argv);

    seqall = ajAcdGetSeqall("sequence");
    score  = ajAcdGetInt("score");
    outf   = ajAcdGetOutfile("outfile");


    substr = ajStrNew();


    while(ajSeqallNext(seqall, &seq))
    {
	begin = ajSeqallGetseqBegin(seqall);
	end   = ajSeqallGetseqEnd(seqall);

	strand = ajSeqGetSeqCopyS(seq);
	ajStrFmtUpper(&strand);

	ajStrAssignSubC(&substr,ajStrGetPtr(strand),begin-1,end-1);

	len=ajStrGetLen(substr);

	ajFmtPrintF(outf,"\n\nNEWCPGSEEK of %s from %d to %d\n",
		    ajSeqGetNameC(seq),begin,begin+len-1);
	ajFmtPrintF(outf,"with score > %d \n\n",score);

	ajFmtPrintF(outf," Begin    End  Score");
	ajFmtPrintF(outf,"        CpG  %%CG  CG/GC\n");

	newcpgseek_cpgsearch(&outf,0,len,ajStrGetPtr(substr),ajSeqGetNameC(seq),
			     begin,score);
	ajFmtPrintF(outf,"-------------------------------------------\n");

	ajStrDel(&strand);
    }


    ajSeqDel(&seq);
    ajStrDel(&substr);
    ajFileClose(&outf);

    ajSeqallDel(&seqall);
    ajFileClose(&outf);

    embExit();

    return 0;
}




/* @funcstatic newcpgseek_cpgsearch *******************************************
**
** Perform cpg search
**
** @param [u] outf [AjPFile*] outfile
** @param [r] from [ajint] start pos
** @param [r] to [ajint] end pos
** @param [r] p [const char*] sequence
** @param [r] name [const char*] sequence name
** @param [r] begin [ajint] start in sequence
** @param [r] score [ajint] score
** @@
******************************************************************************/

static void newcpgseek_cpgsearch(AjPFile *outf, ajint from, ajint to,
				 const char *p, const char *name, ajint begin,
				 ajint score)
{
    ajint i;
    ajint c;
    ajint z;

    ajint sum;
    ajint ssum;
    ajint lsum;
    ajint t;
    ajint top;

    ajint dcg;
    ajint dgc;
    ajint gc;


    for(i=from,c=to-1,sum=ssum=t=top=0,lsum=-1,z=begin-1;i<to;++i,ssum=sum)
    {
	if(p[i]=='C' && p[i+1]=='G' && c-i)
	    sum+=score+1;
	--sum;

	if(sum<0)
	    sum=0;

	if(!sum && ssum)
	{
	    newcpgseek_calcgc(lsum+1,t+2,p,&dcg,&dgc,&gc);
	    if(dgc)
	    {
		ajFmtPrintF(*outf,"%6d %6d %5d ",lsum+2+z,
			    t+2+z,top);
		ajFmtPrintF(*outf,"     %5d %5.1f %6.2f\n",
			    dcg,(float)gc*100.0/(float)(t+1-lsum),
			    (float)(dcg/dgc));
	    }

	    newcpgseek_cpgsearch(outf,t+2,i,p,name,begin,score);
	    sum = ssum = lsum = t = top =0;
	}

	if(sum>top)
	{
	    t = i;
	    top=sum;
	}

	if(!sum)
	    lsum = i;
    }


    if(sum)
    {
	newcpgseek_calcgc(lsum+1,t+2,p,&dcg,&dgc,&gc);
	if(dgc)
	{
	    ajFmtPrintF(*outf,"*%6d %6d %5d ",lsum+2+z,t+2+z,
			top);
	    ajFmtPrintF(*outf,"     %5d %5.1f %6.2f\n",
			dcg,(float)gc*100.0/(float)(t+1-lsum),
			((float)dcg/(float)dgc));
	}

	newcpgseek_cpgsearch(outf,t+2,to,p,name,begin,score);
    }

    return;
}




/* @funcstatic newcpgseek_calcgc **********************************************
**
** Calculate gc content
**
** @param [r] from [ajint] from
** @param [r] to [ajint] to
** @param [r] p [const char*] sequence
** @param [w] dcg [ajint*] number cg
** @param [w] dgc [ajint*] number gc
** @param [w] gc [ajint*] number c or g
** @@
******************************************************************************/

static void newcpgseek_calcgc(ajint from, ajint to, const char *p, ajint *dcg,
			      ajint *dgc, ajint *gc)
{

    ajint i;
    ajint c;

    c = to-1;

    for(i=from,*gc=*dgc=*dcg=0;i<=c;++i)
    {
	if(p[i]=='G' || p[i]=='C')
	    ++*gc;

	if(p[i]=='C' && p[i+1]=='G' && c-i)
	    ++*dcg ;

	if(p[i]=='G' && p[i+1]=='C' && c-i )
	    ++*dgc ;
    }

    return;
}