#include "rlink.h"

#include "GBitVec.h"

#include <float.h>

//import globals from main program:

//extern GffNames* gseqNames;

extern FILE *c_out; // file handle for the input transcripts that are fully covered by reads

extern bool specific;

extern bool trim;

extern bool partialcov;

extern bool complete;

//extern bool debugMode;

//extern bool verbose;

extern bool eonly;

extern int maxReadCov;

extern float isofrac;

extern float mcov;

extern int mintranscriptlen; // minimum number for a transcript to be printed

extern int sensitivitylevel;

extern int junctionsupport; // anchor length for junction to be considered well supported <- consider shorter??

extern int junctionthr; // number of reads needed to support a particular junction

extern float readthr; // read coverage per bundle bp to accept it; otherwise considered noise

extern uint bundledist; // reads at what distance should be considered part of separate bundles

// <- this is not addressed everywhere, e.g. in infer_transcripts -> look into this

extern bool includesource;

extern bool EM;

extern bool weight;

extern FILE* f_out;

extern GStr label;

void printTime(FILE* f) {

time_t ltime; /* calendar time */

ltime=time(NULL);

struct tm *t=localtime(&ltime);

fprintf(f, "[%02d/%02d %02d:%02d:%02d]",t->tm_mon+1, t->tm_mday,

t->tm_hour, t->tm_min, t->tm_sec);

}

void printBitVec(GBitVec& bv) {

for (uint i=0;i<bv.size();i++) {

fprintf(stderr, "%c", bv.test(i)?'1':'0');

}

}

CJunction* add_junction(int start, int end, int leftsupport, int rightsupport,

GList<CJunction>& junction, char strand, int nh) {

int oidx=-1;

CJunction *nj=NULL;

//CJunction* nj=junction.AddIfNew(new CJunction(start, end, strand), true, &oidx);

CJunction jn(start, end, strand);

if (junction.Found(&jn, oidx)) {

nj=junction.Get(oidx);

}

else {

nj=new CJunction(start, end, strand);

junction.Add(nj);

}

//if (nh==0) nh=1;

nj->nreads+=float(1)/nh;

if (leftsupport >= junctionsupport && rightsupport >=junctionsupport) {

nj->nreads_good+=float(1)/nh;

}

return nj;

}

void cov_add(GVec<float>& bpcov, int i, int j, float v) {

if (j>=bpcov.Count()) bpcov.Resize(j+1, 0);

for (int k=i;k<j;k++)

bpcov[k]+=v;

}

float getBCov(GVec<float>& bpcov, int p) {

//if (p<0) GMessage("Error: invalid bpcov index (%d)!\n", p);

if (p>=bpcov.Count()) return 0;

else return bpcov[p];

}

bool maxCovReached(int currentstart, GBamRecord& brec, BundleData& bdata) {

for (int i=0;i<brec.exons.Count();i++) {

if (getBCov(bdata.bpcov, brec.exons[i].start-currentstart)<=maxReadCov)

return false;

}

if (!bdata.covSaturated) {

GMessage("Warning: bundle %s:%d-%d(%d) (%djs) reached coverage saturation (%d) starting with read mapped at %d\n",

bdata.refseq.chars(), bdata.start, bdata.end, bdata.numreads, bdata.junction.Count(),

maxReadCov, brec.start);

bdata.covSaturated=true;

}

return true;

}

void countRead(BundleData& bdata, GBamRecord& brec, int hi) {

static uint32_t BAM_R2SINGLE = BAM_FREAD2 | BAM_FMUNMAP ;

if (hi==0) {

for (int i=0;i<brec.exons.Count();i++) {

bdata.frag_len+=brec.exons[i].len();

}

if (!brec.isPaired() || ((brec.flags()&BAM_FREAD1)!=0) ||

((brec.flags()&BAM_R2SINGLE)==BAM_R2SINGLE ) ) {

bdata.num_fragments++;

}

}

}

int processRead(int currentstart, int currentend, BundleData& bdata,

GHash<int>& hashread, GBamRecord& brec, char strand, int nh, int hi) {

GList<CReadAln>& readlist = bdata.readlist;

GList<CJunction>& junction = bdata.junction;

GVec<float>& bpcov = bdata.bpcov;

int readstart=brec.start;

CReadAln* readaln=NULL;

bool covSaturated=false;

if (bdata.end<currentend) {

bdata.start=currentstart;

bdata.end=currentend;

}

//bool first_mapping = (nh==1 || hi==0);

if (maxReadCov>0) {

covSaturated=maxCovReached(currentstart, brec, bdata);

//if (bdata.firstPass==2 && !covSaturated)

//single pass only -> unconditional read creation

if (!covSaturated)

readaln=new CReadAln(strand, nh, brec.start, brec.end);

}

//if (bdata.firstPass) { //first pass or singlePass use

bdata.numreads++;

//process cigar string

int leftsupport=0;

int rightsupport=brec.exons.Last().len();

int support=0;

int maxleftsupport=0;

int njunc=0;

GVec<int> leftsup;

GVec<int> rightsup;

GPVec<CJunction> newjunction(false);

for (int i=0;i<brec.exons.Count();i++) {

CJunction* nj=NULL;

if (i) {

//deal with introns

GSeg seg(brec.exons[i-1].end, brec.exons[i].start);

leftsupport=brec.exons[i-1].len();

if (leftsupport>maxleftsupport) {

maxleftsupport=leftsupport;

}

leftsup.Add(maxleftsupport);//on the left, always add current max support (?)

support=brec.exons[i].len();

rightsup.Add(support); //..but on the right, real support value is added

/* previously:

if (i==brec.exons.Count()-1)

nj=add_junction(seg.start, seg.end, leftsupport, support, junction, strand, nh);

else

nj=add_junction(seg.start, seg.end, maxleftsupport, support, junction, strand, nh);

*/

int maxrightsupport = support > rightsupport ? support : rightsupport;

nj=add_junction(seg.start, seg.end, maxleftsupport, maxrightsupport, junction, strand, nh);

newjunction.Add(nj);

}

cov_add(bpcov, brec.exons[i].start-currentstart,

brec.exons[i].end-currentstart, float(1)/nh);

if (readaln) {

if (nj) readaln->juncs.Add(nj);

readaln->segs.Add(brec.exons[i]);

}

}

njunc=newjunction.Count();

//--

if (njunc>1) { //2 or more introns

int maxrightsupport=rightsup[njunc-1];

for (int j=njunc-2;j>=0;j--) {

if (rightsup[j]>maxrightsupport) { maxrightsupport=rightsup[j];}

else { //adjust support on the right so that I don't exclude very short internal exons

if (rightsup[j]<junctionsupport && leftsup[j]>=junctionsupport && maxrightsupport>=junctionsupport){

newjunction[j]->nreads_good+=float(1)/nh;

}

}

}

}

if((int)brec.end>currentend) {

//return currentend;

currentend=brec.end;

bdata.end=currentend;

}

//if (bdata.firstPass == 1 || covSaturated) {

if (covSaturated) {

return(currentend); //for 1st-pass only, exit here

}

//} //<-- for single-pass or first pass

//--

//--> 2nd pass or single-pass run (firstPass==0 or firstPass==2)

GStr readname(brec.name());

if (readaln==NULL) {

//2nd pass

//TODO: we should try to see if we can collapse this read first

readaln=new CReadAln(strand, nh, brec.start, brec.end);

}

int n=readlist.Add(readaln);

if (brec.refId()==brec.mate_refId()) {

//only consider mate pairing data if mates are on the same chromosome/contig

//TODO: this should be reconsidered if we decide to care about FUSION transcripts!

GStr readname(brec.name());

int pairstart=brec.mate_start();

GStr id(readname); // init id with readname

id+=':';id+=readstart;id+=':';id+=hi;

if (readstart<=pairstart) {

//only add the first mate in a pair

if (!hashread[id.chars()])

hashread.Add(id.chars(), new int(n));

}

if (readstart>pairstart) {

//must have seen its pair earlier

GStr pairid(readname);

pairid+=':';pairid+=pairstart;pairid+=':';pairid+=hi;

const int* np=hashread[pairid.chars()];

if (np) {

readlist[n]->pair_idx=*np;

readlist[*np]->pair_idx=n;

//we can now discard this pair info

hashread.Remove(pairid.chars());

hashread.Remove(id.chars()); //just in case it exists

}

}

} //<-- if mate is mapped on the same chromosome

//int end=brec.end;

/*

if (bdata.firstPass==0) {//2nd pass only

for (int i=0;i<brec.exons.Count();i++) {

readaln->segs.Add(brec.exons[i]);

if (i) {

int jidx=-1;

CJunction jn(brec.exons[i-1].end, brec.exons[i].start, strand);

if (junction.Found(&jn, jidx)) {

readaln->juncs.Add(junction.Get(jidx));

} else GError("Error: junction %d-%d not found!\n", brec.exons[i-1].end, brec.exons[i].start);

}

}

}

*/

return currentend; //currentend already set earlier

//}

}

int get_min_start(CGroup **currgroup) {

int nextgr=0;

if(currgroup[0]!=NULL) {

if(currgroup[1]!=NULL) {

if(currgroup[2]!=NULL) {

int twogr = currgroup[0]->start < currgroup[1]->start ? 0 : 1;

nextgr = currgroup[twogr]->start < currgroup[2]->start ? twogr : 2;

return(nextgr);

}

else {

nextgr = currgroup[0]->start < currgroup[1]->start ? 0 : 1;

return(nextgr);

}

}

else {

if(currgroup[2]!=NULL) {

nextgr = currgroup[0]->start < currgroup[2]->start ? 0 : 2;

return(nextgr);

}

else {

return(0);

}

}

} // end if(currgroup[0]!=NULL)

else {

if(currgroup[1]!=NULL) {

if(currgroup[2]!=NULL) {

nextgr = currgroup[1]->start < currgroup[2]->start ? 1 : 2;

return(nextgr);

}

else {

return(1);

}

}

else {

return(2);

}

}

return(nextgr);

}

void set_strandcol(CGroup *prevgroup, CGroup *group, int grcol, GVec<int>& eqcol, GVec<int>& equalcolor){

int zerocol=eqcol[prevgroup->color];

if(zerocol>-1) {

while(eqcol[zerocol]!=-1 && eqcol[zerocol]!=zerocol) {

zerocol=eqcol[zerocol];

}

int tmpcol=zerocol;

while(equalcolor[zerocol]!=zerocol) {

zerocol=equalcolor[zerocol];

}

eqcol[prevgroup->color]=zerocol;

eqcol[tmpcol]=zerocol;

if(zerocol<grcol) {

equalcolor[grcol]=zerocol;

group->color=zerocol;

}

else if(grcol<zerocol) {

equalcolor[zerocol]=grcol;

eqcol[prevgroup->color]=grcol;

}

} // if(zerocol>-1)

else {

eqcol[prevgroup->color]=grcol;

}

}

void add_group_to_bundle(CGroup *group, CBundle *bundle, GPVec<CBundlenode>& bnode){

CBundlenode *currlastnode=bnode[bundle->lastnodeid];

int bid=bnode.Count();

if(group->start > currlastnode->end) { // group after last bnode

CBundlenode *currbnode=new CBundlenode(group->start,group->end,group->cov_sum,bid);

currlastnode->nextnode=currbnode;

bnode.Add(currbnode);

bundle->lastnodeid=bid;

bundle->len+=group->end-group->start+1;

bundle->cov+=group->cov_sum;

bundle->nread+=group->nread;

bundle->multi+=group->multi;

}

else { // group overlaps bnode

if(currlastnode->end < group->end) {

bundle->len+= group->end - currlastnode->end;

currlastnode->end= group->end;

}

bundle->cov+=group->cov_sum;

bundle->nread+=group->nread;

bundle->multi+=group->multi;

currlastnode->cov+=group->cov_sum;

}

}

int create_bundle(GPVec<CBundle>& bundle,CGroup *group,GPVec<CBundlenode>& bnode) {

int bid=bnode.Count();

int bno=bundle.Count();

CBundlenode *startbnode=new CBundlenode(group->start,group->end,group->cov_sum,bid);

CBundle *newbundle=new CBundle(group->end-group->start+1,group->cov_sum,group->nread,group->multi,bid,bid);

bundle.Add(newbundle);

bnode.Add(startbnode);

return(bno);

}

int setCmp(const pointer p1, const pointer p2) {

CTrInfo *a=(CTrInfo*)p1;

CTrInfo *b=(CTrInfo*)p2;

if(a->trno<b->trno) return -1;

if(a->trno>b->trno) return 1;

return 0;

}

int capCmp(const pointer p1, const pointer p2) {

CTrInfo *a=(CTrInfo*)p1;

CTrInfo *b=(CTrInfo*)p2;

if(a->abundance<b->abundance) return -1;

if(a->abundance>b->abundance) return 1;

return 0;

}

int trCmp(const pointer p1, const pointer p2) {

CTransfrag *a=(CTransfrag*)p1;

CTransfrag *b=(CTransfrag*)p2;

if(a->pattern.count()<b->pattern.count()) return 1;

if(a->pattern.count()>b->pattern.count()) return -1;

if(!a->real && b->real) return -1;

if(a->real && !b->real) return 1;

if(a->abundance<b->abundance) return 1;

if(a->abundance>b->abundance) return -1;

return 0;

}

int edgeCmp(const pointer p1, const pointer p2) {

CNetEdge *a=(CNetEdge*)p1;

CNetEdge *b=(CNetEdge*)p2;

if(a->rate<b->rate) return 1;

if(a->rate>b->rate) return -1;

return 0;

}

int pointCmp(const pointer p1, const pointer p2) {

CTrimPoint *a=(CTrimPoint*)p1;

CTrimPoint *b=(CTrimPoint*)p2;

if(a->start && !b->start) return -1;

if(!a->start && b->start) return 1;

if(a->abundance<b->abundance) return 1;

if(a->abundance>b->abundance) return -1;

return 0;

}

int edgeCmpEM(const pointer p1, const pointer p2) {

CNetEdge *a=(CNetEdge*)p1;

CNetEdge *b=(CNetEdge*)p2;

if(a->fake && !b->fake) return 1; // check if this is right -> I want the fake one to come last

if(!a->fake && b->fake) return -1;

if(a->rate<b->rate) return 1;

if(a->rate>b->rate) return -1;

return 0;

}

int guideabundCmp(const pointer p1, const pointer p2) {

CTransfrag *a=(CTransfrag*)p1;

CTransfrag *b=(CTransfrag*)p2;

if(a->abundance<b->abundance) return 1;

if(a->abundance>b->abundance) return -1;

if(a->pattern.count()<b->pattern.count()) return 1;

if(a->pattern.count()>b->pattern.count()) return -1;

return 0;

}

int guidedabundCmp(const pointer p1, const pointer p2) {

CGuide *a=(CGuide*)p1;

CGuide *b=(CGuide*)p2;

if(a->trf->abundance<b->trf->abundance) return 1;

if(a->trf->abundance>b->trf->abundance) return -1;

if(a->trf->pattern.count()<b->trf->pattern.count()) return 1;

if(a->trf->pattern.count()>b->trf->pattern.count()) return -1;

return 0;

}

int juncCmpEnd(const pointer p1, const pointer p2) {

CJunction* a=(CJunction*)p1;

CJunction* b=(CJunction*)p2;

if (a->end<b->end) return -1;

if (a->end>b->end) return 1;

if (a->start<b->start) return -1;

if (a->start>b->start) return 1;

return 0;

}

void merge_fwd_groups(GPVec<CGroup>& group, CGroup *group1, CGroup *group2, GVec<int>& merge, GVec<int>& eqcol) {

// get end of group (group1 is assumed to come before group2)

group1->end=group2->end;

// get smallest color of group

while(eqcol[group1->color]!=group1->color) {

group1->color=eqcol[group1->color];

}

while(eqcol[group2->color]!=group2->color) {

group2->color=eqcol[group2->color];

}

if(group1->color<group2->color) {

eqcol[group2->color]=group1->color;

// print STDERR "in merge set eqcol[",$$group2[2],"]=",$$eqcol{$$group1[2]},"\n";

}

else if(group1->color>group2->color) {

eqcol[group1->color]=group2->color;

// print STDERR "in merge set eqcol[",$$group1[2],"]=",$$group2[2],"\n";

group1->color=group2->color;

}

group1->cov_sum+=group2->cov_sum;

group1->next_gr=group2->next_gr; // this is possible because group1->next_gr=group2

merge[group2->grid]=group1->grid;

group1->nread+=group2->nread;

group1->multi+=group2->multi;

// delete group2

group.freeItem(group2->grid);

}

int add_read_to_group(int n,GList<CReadAln>& readlist,int color,GPVec<CGroup>& group,CGroup **allcurrgroup,

CGroup **startgroup,GVec<int> *readgroup,GVec<int>& eqcol,GVec<int>& merge,float& fraglen,uint& fragno) {

int sno=readlist[n]->strand+1; // 0: negative strand; 1: zero strand; 2: positive strand (starting from -1,0,1)

int readcol=color; // start with the new color -> might change

// check if I've seen read's pair and if yes get its readcol; at the least get read's pair strand if available

int np=readlist[n]->pair_idx; // pair read number

if(np>-1 && readlist[np]->nh) { // read pair exists and it wasn't deleted

if(np<n) { // there is a pair and it came before the current read in sorted order of position

// first group of pair read is: $$readgroup[$np][0]

int grouppair=readgroup[np][0];

while( merge[grouppair]!=grouppair) {

grouppair=merge[grouppair];

}

readgroup[np][0]=grouppair;

readcol=group[readgroup[np][0]]->color; // readcol gets assigned the color of the pair's group

while(eqcol[readcol]!=readcol) { // get smallest color

readcol=eqcol[readcol];

}

//print STDERR "Adjust color of group ",$$readgroup[$np][0]," to $readcol\n";

group[readgroup[np][0]]->color=readcol;

}

else { // it's the first time I see the read in the fragment

fragno++;

// see if I have the correct read strand

char snop=readlist[np]->strand+1; // snop is the strand of pair read

if(sno!=snop) { // different strands for read and pair

if(sno==1 && snop!=1) { // read n is on zero (neutral) strand, but pair has strand

readlist[n]->strand=readlist[np]->strand;

sno=snop; // assign strand of pair to read

}

else if(snop==1 && sno!=1) { // pair read np is on zero strand

readlist[np]->strand=readlist[n]->strand;

}

else { // conflicting strands -> un-pair reads in the hope that one is right

readlist[n]->pair_idx=-1;

readlist[np]->pair_idx=-1;

}

}

eqcol.Add(color); // read will keep the new color so we need to add it to the eqcol

color++;

}

} // if(np>-1 && readlist[np]->nh) : read pair exists and it wasn't deleted

else {

fragno++;

eqcol.Add(color); // read will keep the new color so we need to add it to the eqcol

color++;

}

CGroup *currgroup=allcurrgroup[sno];

if(currgroup != NULL) { // this type of group - negative, unknown, or positive - was created before

//set currgroup first

CGroup *lastgroup=NULL;

while(currgroup!=NULL && readlist[n]->start > currgroup->end) { // while read start after the current group's end advance group -> I might have more groups leaving from current group due to splicing

lastgroup=currgroup;

currgroup=currgroup->next_gr;

}

if(currgroup==NULL || readlist[n]->segs[0].end < currgroup->start) // currgroup is null only if we reached end of currgroup list because currgroup is not NULL initially

currgroup=lastgroup;

// now process each group of coordinates individually

CGroup *thisgroup=currgroup;

int ncoord=readlist[n]->segs.Count(); // number of "exons" in read

int lastpushedgroup=-1;

bool added=false;

for(int i=0;i<ncoord;i++) {

fraglen+=readlist[n]->segs[i].len(); // this might be useful to have if I decide not to use the HI tag anymore

// skip groups that are left behind

while(thisgroup!=NULL && readlist[n]->segs[i].start > thisgroup->end) { // find the group where "exon" fits

lastgroup=thisgroup;

thisgroup=thisgroup->next_gr;

}

if(thisgroup && readlist[n]->segs[i].end >= thisgroup->start) { // read overlaps group

// I need to split pairs here if color didn't reach this group: it means there is a gap between these groups and no reads joining them

if(!i && np>-1 && readlist[np]->nh && np<n) { // I only consider first exon here because the rest of the groups need to get the same color

int grouppair=thisgroup->grid;

while( merge[grouppair]!=grouppair) {

grouppair=merge[grouppair];

}

thisgroup->grid=grouppair;

int thiscol=thisgroup->color;

while(eqcol[thiscol]!=thiscol) { // get smallest color

thiscol=eqcol[thiscol];

}

thisgroup->color=thiscol;

if(thiscol!=readcol) { // pair color didn't reach this group

//fprintf(stderr,"Split pairs: %d-%d and %d-%d on strand %d\n",readlist[np]->start,readlist[np]->end,readlist[n]->start,readlist[n]->end,readlist[n]->strand);

readlist[n]->pair_idx=-1;

readlist[np]->pair_idx=-1;

readcol=color; // readcol gets back the new color

eqcol.Add(color);

color++;

}

}

if(!added) { // read gets added only to first group - why ??

thisgroup->nread+=(float)1/readlist[n]->nh;

if(readlist[n]->nh>1) thisgroup->multi+=(float)1/readlist[n]->nh; // this will probably need to be coded differently if I do super-reads, or collapse more than one read into the same nh

added=true;

}

if(readlist[n]->segs[i].start<thisgroup->start) {

thisgroup->start=readlist[n]->segs[i].start;

}

// find end of new group

CGroup *nextgroup=thisgroup->next_gr;

while(nextgroup!=NULL && readlist[n]->segs[i].end >= nextgroup->start) {

merge_fwd_groups(group,thisgroup,nextgroup,merge,eqcol);

nextgroup=thisgroup->next_gr;

}

if(readlist[n]->segs[i].end > thisgroup->end) {

thisgroup->end=readlist[n]->segs[i].end;

}

// get smallest color of group

while(eqcol[thisgroup->color]!=thisgroup->color) {

thisgroup->color=eqcol[thisgroup->color];

}

if(readcol!=thisgroup->color) { // read color is different from group color

if(readcol<thisgroup->color) { // set group color to current read color

eqcol[thisgroup->color]=readcol;

thisgroup->color=readcol;

}

else { // read color is bigger than group

eqcol[readcol]=thisgroup->color;

readcol=thisgroup->color;

}

}

if(thisgroup->grid != lastpushedgroup) {

readgroup[n].Add(thisgroup->grid); // readgroup for read n gets the id of group

lastpushedgroup=thisgroup->grid;

}

thisgroup->cov_sum+=(float)(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)/readlist[n]->nh; // coverage is different than number of reads

} // end if(thisgroup && readlist[n]->segs[i].end >= thisgroup->start)

else { // read is at the end of groups, or read is not overlapping other groups -> lastgroup is not null here because currgroup was not null

// I need to split pairs here because I have no overlap => color didn't reach here for sure if I am at the first exon

if(!i && np>-1 && readlist[np]->nh && np<n) {

//fprintf(stderr,"Split pairs: %d-%d and %d-%d on strand %d\n",readlist[np]->start,readlist[np]->end,readlist[n]->start,readlist[n]->end,readlist[n]->strand);

readlist[n]->pair_idx=-1;

readlist[np]->pair_idx=-1;

readcol=color;

eqcol.Add(color);

color++;

}

int ngroup=group.Count();

float nread=0;

float multi=0;

if(!added) {

nread=(float)1/readlist[n]->nh;

if(readlist[n]->nh>1) multi=(float)1/readlist[n]->nh;

added=true;

}

CGroup *newgroup=new CGroup(readlist[n]->segs[i].start,readlist[n]->segs[i].end,readcol,ngroup,(float)(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)/readlist[n]->nh,nread,multi);

group.Add(newgroup);

merge.Add(ngroup);

lastgroup->next_gr=newgroup; // can lastgroup be null here -> no from the way I got here

newgroup->next_gr=thisgroup;

readgroup[n].Add(ngroup);

lastpushedgroup=ngroup;

thisgroup=lastgroup;

}

} // for(int i=0;i<ncoord;i++)

} // if(currgroup != NULL)

else { // create new group of this type

// I need to split pairs here because I have no overlap => color didn't reach here for sure

if(np>-1 && readlist[np]->nh && np<n) {

//fprintf(stderr,"Split pairs: %d-%d and %d-%d on strand %d\n",readlist[np]->start,readlist[np]->end,readlist[n]->start,readlist[n]->end,readlist[n]->strand);

readlist[n]->pair_idx=-1;

readlist[np]->pair_idx=-1;

readcol=color;

eqcol.Add(color);

color++;

}

int ncoord=readlist[n]->segs.Count();

CGroup *lastgroup=NULL;

float nread=(float)1/readlist[n]->nh;

float multi=0;

if(readlist[n]->nh>1) multi=(float)1/readlist[n]->nh;

for(int i=0;i<ncoord;i++) {

fraglen+=readlist[n]->segs[i].len();

int ngroup=group.Count();

CGroup *newgroup=new CGroup(readlist[n]->segs[i].start,readlist[n]->segs[i].end,readcol,ngroup,(float)(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)/readlist[n]->nh,nread,multi);

nread=0;

multi=0;

group.Add(newgroup);

merge.Add(ngroup);

if(lastgroup!=NULL) {

lastgroup->next_gr=newgroup;

}

else {

currgroup=newgroup;

}

lastgroup=newgroup;

readgroup[n].Add(ngroup);

}

}

allcurrgroup[sno]=currgroup;

if(startgroup[sno]==NULL) startgroup[sno]=currgroup;

return color;

}

CGraphnode *create_graphnode(int s, int g, uint start,uint end,int nodeno,CBundlenode *bundlenode,

GVec<CGraphinfo> **bundle2graph,GPVec<CGraphnode> **no2gnode) {

/*

{ // DEBUG ONLY

fprintf(stderr,"create_graphnode[%d][%d]:%d-%d\n",s,g,start,end);

}

*/

CGraphnode* gnode=new CGraphnode(start,end,nodeno);

CGraphinfo ginfo(g,nodeno);

bundle2graph[s][bundlenode->bid].Add(ginfo);

no2gnode[s][g].Add(gnode);

return(gnode);

}

float compute_chi(GArray<float>& winleft, GArray<float>& winright, float sumleft, float sumright) {

float chi=0;

for(int j=0;j<CHI_WIN;j++) {

float mul=(winleft[j]+ winright[j])/(sumleft+sumright);

float mur=mul;

mul*=sumleft;

mur*=sumright;

chi+= (winleft[j]-mul)/mul+(winright[j]-mur)/mur;

}

return(chi);

}

void find_trims(int refstart,uint start,uint end,GVec<float>& bpcov,uint& sourcestart,float& maxsourceabundance,uint& sinkend,

float& maxsinkabundance){

if(end-start<2*CHI_WIN-1) return;

float sumleft=0;

float sumright=0;

//float maxsinkchi=0;

//float maxsourcechi=0;

float sinkabundance=0;

float sourceabundance=0;

GArray<float> winleft(CHI_WIN,false); // not auto-sort

GArray<float> winright(CHI_WIN,false); // not auto-sort

for(uint i=start;i<=end;i++) {

if(i-start<2*CHI_WIN-1) { // I have to compute the sumleft and sumright first

if(i-start<CHI_WIN) {

sumleft+=bpcov[i-refstart];

winleft.Add(bpcov[i-refstart]);

}

else {

sumright+=bpcov[i-refstart];

winright.Add(bpcov[i-refstart]);

if(i-start==2*CHI_WIN-2) {

winleft.setSorted(true);

winright.setSorted(true);

}

}

}

else { // I can do the actual sumleft, sumright comparision

sumright+=bpcov[i-refstart];

winright.Add(bpcov[i-refstart]);

float chi=0;

if(sumleft!=sumright) chi=compute_chi(winleft,winright,sumleft,sumright);

if(chi>CHI_THR) { // there is a significant difference

if(sumleft>sumright) { // possible drop (sink cut)

//if(chi>maxsinkchi) {

sinkabundance=(sumleft-sumright)/CHI_WIN;

if(maxsinkabundance<sinkabundance) {

sinkend=i-CHI_WIN;

//maxsinkchi=chi;

maxsinkabundance=sinkabundance;

}

}

else if(sumright>sumleft) { // increase (source cut)

//if(chi>maxsourcechi) {

sourceabundance=(sumright-sumleft)/CHI_WIN;

if(maxsourceabundance<sourceabundance) {

sourcestart=i-CHI_WIN+1;

maxsourceabundance=sourceabundance;

//maxsourcechi=chi;

}

}

}

sumleft-=bpcov[i-refstart-2*CHI_WIN+1];

int idx=winleft.IndexOf(bpcov[i-refstart-2*CHI_WIN+1]);

winleft.Delete(idx);

sumleft+=bpcov[i-refstart-CHI_WIN+1];

winleft.Add(bpcov[i-refstart-CHI_WIN+1]);

sumright-=bpcov[i-refstart-CHI_WIN+1];

idx=winright.IndexOf(bpcov[i-refstart-CHI_WIN+1]);

winright.Delete(idx);

}

}

}

CGraphnode *add_trim_to_graph(int s, int g,uint lastpos,CTrimPoint& mytrim,CGraphnode *graphnode,CGraphnode *source,CGraphnode *sink,GVec<float>& futuretr,

int& graphno,CBundlenode *bundlenode,GVec<CGraphinfo> **bundle2graph,GPVec<CGraphnode> **no2gnode) {

if(mytrim.start) { // this is a source link

float tmp=graphno-1;

CGraphnode *prevnode=NULL;

if(mytrim.pos>graphnode->start) { // there is place for another node

uint prevend=graphnode->end;

graphnode->end=mytrim.pos-1;

prevnode=graphnode;

graphnode=create_graphnode(s,g,mytrim.pos,prevend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

}

source->child.Add(graphnode->nodeid); // this node is the child of source

graphnode->parent.Add(source->nodeid); // this node has source as parent

if(prevnode) {

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

}

futuretr.Add(tmp);

futuretr.cAdd(0.0);

futuretr.Add(mytrim.abundance);

}

else { // this is a link to sink

float tmp=graphno-1;

if(mytrim.pos<lastpos) { // there is still place for a new node

uint prevend=graphnode->end;

graphnode->end=mytrim.pos;

CGraphnode *prevnode=graphnode;

graphnode=create_graphnode(s,g,mytrim.pos+1,prevend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

sink->parent.Add(prevnode->nodeid); // prevnode is the parent of sink

}

else {

sink->parent.Add(graphnode->nodeid); // prevnode is the parent of sink

}

futuretr.Add(tmp);

futuretr.cAdd(1.0);

futuretr.Add(mytrim.abundance);

}

return(graphnode);

}

CGraphnode *trimnode(int s, int g, int refstart,uint newend, CGraphnode *graphnode,CGraphnode *source, CGraphnode *sink, GVec<float>& bpcov,

GVec<float>& futuretr, int& graphno,CBundlenode *bundlenode,GVec<CGraphinfo> **bundle2graph,GPVec<CGraphnode> **no2gnode) {

uint sourcestart=0;

uint sinkend=0;

float sinkabundance=0;

float sourceabundance=0;

find_trims(refstart,graphnode->start,newend,bpcov,sourcestart,sourceabundance,sinkend,sinkabundance);

if(sourcestart < sinkend) { // source trimming comes first

if(sourcestart) { // there is evidence of graphnode trimming from source

graphnode->end=sourcestart-1;

CGraphnode *prevnode=graphnode;

graphnode=create_graphnode(s,g,sourcestart,newend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

source->child.Add(graphnode->nodeid); // this node is the child of source

graphnode->parent.Add(source->nodeid); // this node has source as parent

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

float tmp=graphno-1;

futuretr.cAdd(0.0);

futuretr.Add(tmp);

sourceabundance+=trthr;futuretr.Add(sourceabundance);

tmp=prevnode->nodeid;futuretr.Add(tmp);

tmp=graphnode->nodeid;futuretr.Add(tmp);

tmp=trthr;futuretr.Add(tmp);

}

// sinkend is always positive since it's bigger than sourcestart

float tmp=graphno-1;

graphnode->end=sinkend;

CGraphnode *prevnode=graphnode;

graphnode=create_graphnode(s,g,sinkend+1,newend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

sink->parent.Add(prevnode->nodeid); // prevnode is the parent of sink

// remember to create transfrag as well -> I don't know the gno yet, so I can not create it here

futuretr.Add(tmp);

futuretr.cAdd(-1.0);

sinkabundance+=trthr;futuretr.Add(sinkabundance);

tmp=prevnode->nodeid;futuretr.Add(tmp);

tmp=graphnode->nodeid;futuretr.Add(tmp);

tmp=trthr;futuretr.Add(tmp);

}

else if(sourcestart > sinkend) { // sink trimming comes first

if(sinkend) {

graphnode->end=sinkend;

CGraphnode *prevnode=graphnode;

graphnode=create_graphnode(s,g,sinkend+1,newend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

sink->parent.Add(prevnode->nodeid); // prevnode is the parent of sink

// remember to create transfrag as well -> I don't know the gno yet, so I can not create it here

float tmp=graphno-2;

futuretr.Add(tmp);

futuretr.cAdd(-1.0);

sinkabundance+=trthr;futuretr.Add(sinkabundance);

tmp=prevnode->nodeid;futuretr.Add(tmp);

tmp=graphnode->nodeid;futuretr.Add(tmp);

tmp=trthr;futuretr.Add(tmp);

}

// sourcestart is positive since it's bigger than sinkend

graphnode->end=sourcestart-1;

CGraphnode *prevnode=graphnode;

graphnode=create_graphnode(s,g,sourcestart,newend,graphno,bundlenode,bundle2graph,no2gnode);

graphno++;

source->child.Add(graphnode->nodeid); // this node is the child of source

graphnode->parent.Add(source->nodeid); // this node has source as parent

prevnode->child.Add(graphnode->nodeid); // this node is the child of previous node

graphnode->parent.Add(prevnode->nodeid); // this node has as parent the previous node

float tmp=graphno-1;

futuretr.cAdd(0.0);

futuretr.Add(tmp);

sourceabundance+=trthr;futuretr.Add(sourceabundance);

tmp=prevnode->nodeid;futuretr.Add(tmp);

tmp=graphnode->nodeid;futuretr.Add(tmp);

tmp=trthr;futuretr.Add(tmp);

}

// else both source and sink trimming are not present

return(graphnode);

}

inline int edge(int min, int max, int gno) {

//return((gno-1)*min-min*(min-1)/2+max-min); // this should be changed if source to node edges are also stored

return((gno-1)*(min+1)-min*(min-1)/2+max-min); // this includes source to node edges

}

GBitVec traverse_dfs(int s,int g,CGraphnode *node,CGraphnode *sink,GBitVec parents,int gno, GVec<bool>& visit,

GPVec<CGraphnode> **no2gnode,GPVec<CTransfrag> **transfrag){

if(visit[node->nodeid]) {

node->parentpat = node->parentpat | parents;

for(int n=0;n<gno;n++) {

if(parents[n]) // add node's children to all parents of node

no2gnode[s][g][n]->childpat = no2gnode[s][g][n]->childpat | node->childpat;

else if(node->childpat[n])

no2gnode[s][g][n]->parentpat = no2gnode[s][g][n]->parentpat | node->parentpat;

}

}

else {

node->childpat.resize(1+gno*(gno+1)/2);

node->parentpat.resize(1+gno*(gno+1)/2);

node->parentpat = node->parentpat | parents;

visit[node->nodeid]=true;

parents[node->nodeid]=1; // add the node to the parents

if(node->parent.Count()==1 && !node->parent[0]) { // node has source only as parent -> add transfrag from source to node

GBitVec trpat(1+gno*(gno+1)/2);

trpat[0]=1;

trpat[node->nodeid]=1;

trpat[edge(0,node->nodeid,gno)]=1;

GVec<int> nodes;

nodes.cAdd(0);

nodes.Add(node->nodeid);

CTransfrag *tr=new CTransfrag(nodes,trpat,trthr);

/*

{ // DEBUG ONLY

fprintf(stderr,"Add source transfrag[%d][%d]= %d and pattern",s,g,transfrag[s][g].Count());