#include "rlink.h"

#include "GBitVec.h"

#include <float.h>

//#define GMEMTRACE 1 //debugging mem allocation

#ifdef GMEMTRACE

#include "proc_mem.h"

#endif

//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; // this is false if I use guides that are fragments, but I want them to be true

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 uint 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 bool geneabundance; // need to compute the gene abundance

extern float fpkm_thr;

extern float tpm_thr;

extern bool enableNames;

extern bool includecov;

extern bool retained_intron;

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;

}

*/

CJunction* add_junction(int start, int end, GList<CJunction>& junction, char strand) {

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);

}

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;

}

*/

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

bool neutral=false;

if(sno!=1) neutral=true; // why is neutral true here: because if the sno is -/+ than I want to add their counts to bpcov[1] too

if (j>=bpcov[sno].Count())

for(int s=0;s<3;s++) bpcov[s].Resize(j+1, 0);

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

bpcov[sno][k]+=v;

if(neutral) bpcov[1][k]+=v; // neutral (stranded) gets added twice here

}

}

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) { // coverage saturation is reached if read spans no portion of saturation less than maxReadCov

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 countFragment(BundleData& bdata, GBamRecord& brec, int nh) {

static uint32_t BAM_R2SINGLE = BAM_FREAD2 | BAM_FMUNMAP ;

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

bdata.frag_len+=float(1)*brec.exons[i].len()/nh;

}

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

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

bdata.num_fragments+=float(1)/nh;

}

/*

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++;

}

}

else if (hi==1) {

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

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

}

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

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

bdata.num_fragments1++;

}

}

*/

}

bool exonmatch(GVec<GSeg> &prevexons, GVec<GSeg> &exons) {

if(prevexons.Count() != exons.Count()) return false;

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

if(prevexons[i].end!=exons[i].end || prevexons[i].start!=exons[i].start) return false;

}

return true;

}

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

GHash<int>& hashread, GReadAlnData& alndata) { // some false positives should be eliminated here in order to break the bundle

GList<CReadAln>& readlist = bdata.readlist; // list of reads gathered so far

GList<CJunction>& junction = bdata.junction; // junctions added so far

GBamRecord& brec=*(alndata.brec); // bam record

char strand=alndata.strand;

int nh=alndata.nh;

int hi=alndata.hi;

int readstart=brec.start;

CReadAln* readaln=NULL; // readaln is initialized with NULL

//bool covSaturated=false; // coverage is set to not saturated

bool match=false; // true if current read matches a previous read

int n=readlist.Count()-1;

if(!mergeMode) while(n>-1 && readlist[n]->start==brec.start) {

if(strand==readlist[n]->strand) match=exonmatch(readlist[n]->segs,brec.exons);

if(match) break; // this way I make sure that I keep the n of the matching readlist

n--;

}

else if((alndata.tinfo->cov>=0 && alndata.tinfo->cov<readthr) || (alndata.tinfo->fpkm>=0 && alndata.tinfo->fpkm<fpkm_thr) ||

(alndata.tinfo->tpm>=0 && alndata.tinfo->tpm < tpm_thr)) return; // do not store 'read' if it doesn't meet minimum criteria

else { //mergeMode but the read is above thresholds

if(alndata.tinfo->cov>=0) bdata.covflags |= IS_COV_FLAG;

if(alndata.tinfo->fpkm>=0) bdata.covflags |= IS_FPKM_FLAG;

if(alndata.tinfo->tpm>=0) bdata.covflags |= IS_TPM_FLAG;

}

if (bdata.end<currentend) {// I am not sure why this is done here?

bdata.start=currentstart;

bdata.end=currentend;

}

bdata.numreads++; // number of reads gets increased no matter what

//bdata.wnumreads+=float(1)/nh;

if (!match) { // if this is a new read I am seeing I need to set it up

if(mergeMode && mintranscriptlen) {

int len=0;

for (int i=0;i<brec.exons.Count();i++) len+=brec.exons[i].len();

if(len<mintranscriptlen) return;

}

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

alndata.tinfo=NULL; //alndata.tinfo was passed to CReadAln

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

readaln->len+=brec.exons[i].len();

if(i) {

CJunction* nj=add_junction(brec.exons[i-1].end, brec.exons[i].start, junction, strand);

if (alndata.juncs.Count())

nj->guide_match=alndata.juncs[i-1]->guide_match;

if (nj) {

readaln->juncs.Add(nj);

}

}

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

}

n=readlist.Add(readaln); // reset n for the case there is no match

}

else { //redundant read alignment matching a previous alignment

// keep shortest nh so that I can see for each particular read the multi-hit proportion

if(nh<readlist[n]->nh) readlist[n]->nh=nh;

/*

//for mergeMode, we have to free the transcript info:

if (alndata.tinfo!=NULL) {

delete alndata.tinfo;

alndata.tinfo=NULL;

}

*/

}

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

currentend=brec.end;

bdata.end=currentend;

}

float rdcount=float(1)/nh;

readlist[n]->read_count+=rdcount; // increase single count just in case I don't see the pair

// store the mismatch count per junction so that I can eliminate it later

int nm=brec.tag_int("NM"); // read mismatch

if(!nm) {

nm=brec.tag_int("nM"); // paired mismatch : big problem with STAR alignments

if(brec.isPaired()) nm/=2;

}

// this counts reads over allowed fraction of mismatches and multi-mapped reads

if(nh>1 || nm/readlist[n]->len>mismatchfrac) for(int i=0;i<readlist[n]->juncs.Count();i++) {

readlist[n]->juncs[i]->nm+=rdcount;

}

// now set up the pairing

if (brec.refId()==brec.mate_refId()) { //only consider mate pairing data if mates are on the same chromosome/contig

int pairstart=brec.mate_start();

if (currentstart<=pairstart) { // if pairstart is in a previous bundle I don't care about it

GStr readname(brec.name());

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

//id+=':';id+=readstart; // this shouldn't be needed since I have the HI tag

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

if(pairstart<=readstart) { // if I've seen the pair already

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

if(np) { // the pair was stored --> why wouldn't it be? : only in the case that the pair starts at the same position

bool notfound=true;

for(int i=0;i<readlist[*np]->pair_idx.Count();i++)

if(readlist[*np]->pair_idx[i]==n) {

readlist[*np]->pair_count[i]+=rdcount;

notfound=false;

break;

}

if(notfound) { // I didn't see the pairing before

readlist[*np]->pair_idx.Add(n);

readlist[*np]->pair_count.Add(rdcount);

}

notfound=true;

for(int i=0;i<readlist[n]->pair_idx.Count();i++)

if(readlist[n]->pair_idx[i]==*np) {

readlist[n]->pair_count[i]+=rdcount;

notfound=false;

break;

}

if(notfound) { // I didn't see the pairing before

int i=*np;

readlist[n]->pair_idx.Add(i);

readlist[n]->pair_count.Add(rdcount);

}

hashread.Remove(id.chars());

}

/* do not allow pairings starting at the same place --> way too confusing and I don't gain much info from these

else if(readstart==pairstart) { // I have not seen the pair because it starts at the same place -> need to store it

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

}

*/

}

else { // I might still see the pair in the future

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

}

}

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

}

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, uint bundledist){

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

int bid=bnode.Count();

if(group->start > currlastnode->end + bundledist) { // 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 within bundledist

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 mgtrnodeCmp(const pointer p1, const pointer p2) {

CMTransfrag *a=(CMTransfrag*)p1;

CMTransfrag *b=(CMTransfrag*)p2;

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

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

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

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

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

if(a->transfrag->abundance>b->transfrag->abundance) return -1;*/

return 0;

}

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

CMTransfrag *a=(CMTransfrag*)p1;

CMTransfrag *b=(CMTransfrag*)p2;

if(!a->transfrag->real && b->transfrag->real) return 1; // guides come first

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

if(a->transfrag->abundance<b->transfrag->abundance) return 1; // most abundant transcript comes first

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

if(a->transfrag->pattern.count()<b->transfrag->pattern.count()) return 1; // the one with largest pattern comes first

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

// something to decide equal cases

if(a->transfrag->nodes.Count()<b->transfrag->nodes.Count()) return 1; // the one with most nodes comes first

if(a->transfrag->nodes.Count()>b->transfrag->nodes.Count()) return -1;

int i=0;

while(i<a->transfrag->nodes.Count()) {

if(a->transfrag->nodes[i]<b->transfrag->nodes[i]) return 1; // the transcript starting more to the left comes first

if(a->transfrag->nodes[i]>b->transfrag->nodes[i]) return -1;

i++;

}

return 0;

}

int partguideCmp(const pointer p1, const pointer p2) { // sorting partial guides in order of preference

CPartGuide *a=(CPartGuide*)p1;

CPartGuide *b=(CPartGuide*)p2;

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

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

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

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

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

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

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

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

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

if(a->glen>b->glen) 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 merge_read_to_group(int n,int np, int p, float readcov, int sno,int readcol,GList<CReadAln>& readlist,int color,GPVec<CGroup>& group,CGroup **allcurrgroup,

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

CGroup **startgroup,GVec<int> *readgroup,GVec<int>& eqcol,GVec<int>& merge,int *usedcol) {

//fprintf(stderr,"merge readcol=%d\n",readcol);

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+=readcov*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[p]=-1;

for(int j=0;j<readlist[np]->pair_idx.Count();j++)

if(readlist[np]->pair_idx[j]==n) {

readlist[np]->pair_idx[j]=-1;

break;

}

readcol=usedcol[sno]; // readcol gets back the new color

if(readcol<0) {

usedcol[sno]=color;

readcol=color;

eqcol.Add(color);

color++;

}

}

}

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

thisgroup->nread+=readcov;

if(readlist[n]->nh>1) thisgroup->multi+=readcov; // 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+=(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)*readcov; // 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[p]=-1;

for(int j=0;j<readlist[np]->pair_idx.Count();j++)

if(readlist[np]->pair_idx[j]==n) {

readlist[np]->pair_idx[j]=-1;

break;

}

readcol=usedcol[sno]; // readcol gets back the new color

if(readcol<0) {

usedcol[sno]=color;

readcol=color;

eqcol.Add(color);

color++;

}

}

int ngroup=group.Count();

float nread=0;

float multi=0;

if(!added) {

nread=readcov;

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

added=true;

}

CGroup *newgroup=new CGroup(readlist[n]->segs[i].start,readlist[n]->segs[i].end,readcol,ngroup,(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)*readcov,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[p]=-1;

for(int j=0;j<readlist[np]->pair_idx.Count();j++)

if(readlist[np]->pair_idx[j]==n) {

readlist[np]->pair_idx[j]=-1;

break;

}

readcol=usedcol[sno]; // readcol gets back the new color

if(readcol<0) {

usedcol[sno]=color;

readcol=color;

eqcol.Add(color);

color++;

}

}

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

CGroup *lastgroup=NULL;

float nread=readcov;

float multi=0;

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

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

//fraglen+=readcov*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,(readlist[n]->segs[i].end-readlist[n]->segs[i].start+1)*readcov,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;

}

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) {

CGroup **startgroup,GVec<int> *readgroup,GVec<int>& eqcol,GVec<int>& merge) {

int usedcol[3]={-1,-1,-1};

float single_count=readlist[n]->read_count; // need to compute read single count

if(!mergeMode) for(int p=0;p<readlist[n]->pair_idx.Count();p++) {

// for each pair I have to pretend is independent of the previous one since I grouped together several reads which might be independent

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

// 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[p]; // pair read number

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

// 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) { // read n is on zero (neutral) strand, but pair has strand

// readlist[n]->strand=readlist[np]->strand; I can not update the read strand anymore -> REMEMBER this later

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

}

else if(snop!=1) { // conflicting strands -> un-pair reads in the hope that one is right

readlist[n]->pair_idx[p]=-1;

for(int j=0;j<readlist[np]->pair_idx.Count();j++)

if(readlist[np]->pair_idx[j]==n) {

readlist[np]->pair_idx[j]=-1;

break;

}

np=-1;

}

}

if(np>-1) { // reads didn't get split

single_count-=readlist[n]->pair_count[p]; // update single count of read here

int readcol=usedcol[sno];

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+=readlist[n]->pair_count[p];

if(usedcol[sno]<0) { // I didn't use the color yet

usedcol[sno]=color;

readcol=color;

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

color++;

}

}

color=merge_read_to_group(n,np,p,readlist[n]->pair_count[p],sno,readcol,readlist,color,group,

allcurrgroup,startgroup,readgroup,eqcol,merge,usedcol);

} // this ends if(np>-1)

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

} // for(int p=0,...)

// now I need to deal with the single count

if(mergeMode || single_count>epsilon) { // my way of controlling for rounding errors

//fragno+=single_count;

int readcol=usedcol[readlist[n]->strand+1];

if(readcol<0) {

readcol=color;

usedcol[readlist[n]->strand+1]=color;

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

color++;

}

color=merge_read_to_group(n,-1,-1,single_count,readlist[n]->strand+1,readcol,readlist,color,group,

allcurrgroup,startgroup,readgroup,eqcol,merge,usedcol);

}

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;

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

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

}

return(chi);

}

void find_trims(int refstart,int sno,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

float cov;

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

cov=bpcov[sno][i-refstart];

if(bpcov[1][i-refstart]>cov) { // it means bpcov[1]>0 since bpcov can not be negative and sno!=1 because it can't be neutral group

//cov+=(bpcov[1][i-refstart]-bpcov[0][i-refstart]-bpcov[2][i-refstart])/bpcov[1][i-refstart];

if(bpcov[2-sno][i-refstart]) // I have stranded coverage on different strand

cov+=(bpcov[1][i-refstart]-bpcov[0][i-refstart]-bpcov[2][i-refstart])*cov/(cov+bpcov[2-sno][i-refstart]);

else cov=bpcov[1][i-refstart];

}

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

if(i-start<CHI_WIN) {

sumleft+=cov;

winleft.Add(cov);

}

else {

sumright+=cov;

winright.Add(cov);

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

winleft.setSorted(true);

winright.setSorted(true);

}

}

}

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

sumright+=cov;