/* mafOrtholog - find orthlogs in other species based on maf alignment and reference genePred */

/* Copyright (C) 2011 The Regents of the University of California 
 * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */

#include "common.h"
#include "linefile.h"
#include "hash.h"
#include "options.h"
#include "obscure.h"
#include "jksql.h"
#include "maf.h"
#include "hdb.h"
#include "hgMaf.h"
#include "genePred.h"
#include "genePredReader.h"



char *nibDir;
char *geneTable = "refGene";
char *altTable = "knownGene";
struct hash *geneListHash = NULL ; /* hash of list of genes in each species */

void usage()
/* Explain usage and exit. */
{
errAbort(
  "mafOrtholog - find orthlogs in other species based on maf alignment and reference genePred \n"
  "usage:\n"
  "   mafOrtholog database mafTable in.gp out.tab\n"
  "options:\n"
  "   -table=table - name of table to find orthlogs in other species [default: refGene]\n"
  "   -alt=table - name of alternate table to find orthlogs in other species if table above does not exist [default: knownGene]\n"
  "   -orgs=org.txt - File with list of databases/organisms in order\n"
  "   -thickOnly - Only extract subset between thickStart/thickEnd\n"
  "   -meFirst - Put native sequence first in maf\n"
  );
}

static struct optionSpec options[] = {
   {"orgs", OPTION_STRING},
   {"nibDir", OPTION_STRING},
   {"table", OPTION_STRING},
   {"alt", OPTION_STRING},
   {"thickOnly", OPTION_BOOLEAN},
   {"meFirst", OPTION_BOOLEAN},
   {NULL, 0},
};

boolean thickOnly = FALSE;
boolean meFirst = FALSE;

struct mafAli *mafFromBed12(char *database, char *track, struct bed *bed, 
	struct slName *orgList)
/* Construct a maf out of exons in bed. */
{
/* Loop through all block in bed, collecting a list of mafs, one
 * for each block.  While we're at make a hash of all species seen. */
struct hash *speciesHash = hashNew(0);
struct mafAli *mafList = NULL, *maf, *bigMaf;
struct mafComp *comp, *bigComp;
int totalTextSize = 0;
int i;
for (i=0; i<bed->blockCount; ++i)
    {
    int start = bed->chromStart + bed->chromStarts[i];
    int end = start + bed->blockSizes[i];
    if (thickOnly)
        {
	start = max(start, bed->thickStart);
	end = min(end, bed->thickEnd);
	}
    if (start < end)
        {
	maf = hgMafFrag(database, track, bed->chrom, start, end, '+',
	   database, NULL);
	slAddHead(&mafList, maf);
	for (comp = maf->components; comp != NULL; comp = comp->next)
	    hashStore(speciesHash, comp->src);
	totalTextSize += maf->textSize; 
	}
    }
slReverse(&mafList);

/* Add species in order list too */
struct slName *org;
for (org = orgList; org != NULL; org = org->next)
    hashStore(speciesHash, org->name);

/* Allocate memory for return maf that contains all blocks concatenated together. 
 * Also fill in components with any species seen at all. */
AllocVar(bigMaf);
bigMaf->textSize = totalTextSize;
struct hashCookie it = hashFirst(speciesHash);
struct hashEl *hel;
while ((hel = hashNext(&it)) != NULL)
    {
    AllocVar(bigComp);
    bigComp->src = cloneString(hel->name);
    bigComp->text = needLargeMem(totalTextSize + 1);
    memset(bigComp->text, '.', totalTextSize);
    bigComp->text[totalTextSize] = 0;
    bigComp->strand = '+';
    bigComp->srcSize = totalTextSize;	/* It's safe if a bit of a lie. */
    hel->val = bigComp;
    slAddHead(&bigMaf->components, bigComp);
    }

/* Loop through maf list copying in data. */
int textOffset = 0;
for (maf = mafList; maf != NULL; maf = maf->next)
    {
    for (comp = maf->components; comp != NULL; comp = comp->next)
        {
	bigComp = hashMustFindVal(speciesHash, comp->src);
	memcpy(bigComp->text + textOffset, comp->text, maf->textSize);
	bigComp->size += comp->size;
	}
    textOffset += maf->textSize;
    }

/* Cope with strand of darkness. */
if (bed->strand[0] == '-')
    {
    for (comp = bigMaf->components; comp != NULL; comp = comp->next)
	reverseComplement(comp->text, bigMaf->textSize);
    }

/* If got an order list then reorder components according to it. */
if (orgList != NULL)
    {
    struct mafComp *newList = NULL;
    for (org = orgList; org != NULL; org = org->next)
        {
	comp = hashMustFindVal(speciesHash, org->name);
	slAddHead(&newList, comp);
	}
    slReverse(&newList);
    bigMaf->components = newList;
    }

/* Rename our own component to bed name */
comp = hashMustFindVal(speciesHash, database);
freeMem(comp->src);
comp->src = cloneString(bed->name);


/* Clean up and go home. */
hashFree(&speciesHash);
mafAliFreeList(&mafList);
return bigMaf;
}

void moveMeToFirst(struct mafAli *maf, char *myName)
/* Find component matching myName, and move it to first. */
{
struct mafComp *comp;
for (comp = maf->components; comp != NULL; comp = comp->next)
    if (sameString(comp->src, myName))
        break;
assert(comp != NULL);
slRemoveEl(&maf->components, comp);
slAddHead(&maf->components, comp);
}

struct genePred *getOverlappingGeneDb(struct genePred **list, char *table, char *chrom, int cStart, int cEnd, char *name, int *retOverlap, char *db)
{
/* read all genes from a table find the gene with the biggest overlap. 
   Cache the list of genes to so we only read it once */

struct genePred *el = NULL, *bestMatch = NULL, *gp = NULL;
int overlap = 0 , bestOverlap = 0, i;
int *eFrames;

if (list == NULL)
    return NULL;

if (*list == NULL)
    {
    struct genePred *gpList = NULL;
    struct sqlConnection *conn = sqlConnect(db);
    struct genePredReader *gpr = NULL;
    if (!hTableExistsDb(db,table))
        table = altTable;
    if (!hTableExistsDb(db,table))
        {
        verbose(2,"no table %s in %s\n",table, db);
        return NULL;
        }
    gpr = genePredReaderQuery(conn, table, NULL);
    verbose(1,"Loading Predictions from %s in %s\n",table, db);
    gpList = genePredReaderAll(gpr);
    if (gpList != NULL)
        {
        hashAdd(geneListHash, db, gpList);
        *list = gpList;
        }
    sqlDisconnect(&conn);
    }
for (el = *list; el != NULL; el = el->next)
    {
    if (chrom != NULL && el->chrom != NULL)
        {
        overlap = 0;
        if ( sameString(chrom, el->chrom))
            {
            for (i = 0 ; i<(el->exonCount); i++)
                {
                overlap += positiveRangeIntersection(cStart,cEnd, el->exonStarts[i], el->exonEnds[i]) ;
                }
            if (overlap > 20 && sameString(name, el->name))
                {
                bestMatch = el;
                bestOverlap = overlap;
                *retOverlap = bestOverlap;
                }
            if (overlap > bestOverlap)
                {
                bestMatch = el;
                bestOverlap = overlap;
                *retOverlap = bestOverlap;
                }
            }
        }
    }
if (bestMatch != NULL)
    {
    /* Allocate genePred and fill in values. */
    AllocVar(gp);
    gp->name = cloneString(bestMatch->name);
    gp->chrom = cloneString(bestMatch->chrom);
    gp->strand[1] = bestMatch->strand[1];
    gp->strand[0] = bestMatch->strand[0];
    gp->txStart = bestMatch->txStart;
    gp->txEnd = bestMatch->txEnd;
    gp->cdsStart = bestMatch->cdsStart;
    gp->cdsEnd = bestMatch->cdsEnd;
    gp->exonCount = bestMatch->exonCount;
    AllocArray(gp->exonStarts, bestMatch->exonCount);
    AllocArray(gp->exonEnds, bestMatch->exonCount);
    for (i=0; i<bestMatch->exonCount; ++i)
        {
        gp->exonStarts[i] = bestMatch->exonStarts[i] ;
        gp->exonEnds[i] = bestMatch->exonEnds[i] ;
        }
    gp->optFields = bestMatch->optFields;
    gp->id = bestMatch->id;

    if (bestMatch->optFields & genePredName2Fld)
        gp->name2 = cloneString(bestMatch->name2);
    else
        gp->name2 = NULL;
    if (bestMatch->optFields & genePredCdsStatFld)
        {
        gp->cdsStartStat = bestMatch->cdsStartStat;
        gp->cdsEndStat = bestMatch->cdsEndStat;
        }
    if (bestMatch->optFields & genePredExonFramesFld)
        {
        gp->exonFrames = AllocArray(eFrames, bestMatch->exonCount);
        for (i = 0; i < bestMatch->exonCount; i++)
            gp->exonFrames[i] = bestMatch->exonFrames[i];
        }
    eFrames = gp->exonFrames;
    }

return gp;
}
void printOrthologs(FILE *f, struct mafAli *maf, struct genePred *gp)
{
struct mafComp *c;

fprintf(f, "%s", gp->name);
if (maf == NULL || maf->components == NULL)
    {
    fprintf(f,"\n");
    return;
    }

if (maf->components->next == NULL)
    {
    fprintf(f,"\n");
    return;
    }

for (c = maf->components->next; c != NULL ; c=c->next)
    {
    int end = c->start + c->size;
    int geneOverlap = 0;
    char *src[20];
    struct genePred *og = NULL;
    int wordCount;
    struct genePred *geneList = NULL;
    wordCount = chopString(c->src, ".", src, ArraySize(src));
    verbose(3, "searching in %s at chrom %s %d %d, %s\n",src[0], src[1], c->start, end, gp->name);
    assert(wordCount == 2);
    geneList = (struct genePred *)hashFindVal(geneListHash, src[0]);
    og = getOverlappingGeneDb(&geneList, geneTable, src[1], c->start, 
        end , gp->name, &geneOverlap, src[0]);
    if (og != NULL)
        fprintf(f,"\t%s\t%s\t%s\t%d\t%d",c->src, og->name, og->chrom, og->txStart, og->txEnd);
    }
    fprintf(f,"\n");
}

void mafOrtholog(char *database, char *track, char *genePredFile, char *outFile)
/* mafOrtholog - find orthlogs in other species based on maf alignment and reference genePred */
{
struct slName *orgList = NULL;
FILE *f = mustOpen(outFile, "w");
struct genePredReader *gpr = genePredReaderFile(genePredFile, NULL);
struct genePred *gpList = genePredReaderAll(gpr), *gp = NULL;
struct sqlConnection *conn = hAllocConn();

if (optionExists("nibDir"))
    nibDir = optionVal("nibDir", NULL);
if (optionExists("orgs"))
    {
    char *orgFile = optionVal("orgs", NULL);
    char *buf;
    readInGulp(orgFile, &buf, NULL);
    orgList = stringToSlNames(buf);
    }

for (gp = gpList ; gp != NULL ; gp=gp->next)
    {
    struct mafAli *maf = NULL;
    if (thickOnly)
        maf = mafLoadInRegion(conn, track,
        gp->chrom, gp->cdsStart, gp->cdsEnd);
    else
        maf = mafLoadInRegion(conn, track,
        gp->chrom, gp->txStart, gp->txEnd);
    if (meFirst)
        moveMeToFirst(maf, database);
    printOrthologs(f, maf, gp);
    mafAliFree(&maf);
    }
carefulClose(&f);
}

int main(int argc, char *argv[])
/* Process command line. */
{
optionInit(&argc, argv, options);
if (argc != 5)
    usage();
hSetDb(argv[1]);
geneListHash = hashNew(0);
geneTable = optionVal("table", geneTable);
altTable = optionVal("alt", altTable);
thickOnly = optionExists("thickOnly");
meFirst = optionExists("meFirst");
mafOrtholog(argv[1], argv[2], argv[3], argv[4]);
return 0;
}