//#include "tablemaker.h"

#include "rlink.h"

#include <numeric>

extern GStr ballgown_dir;

int rc_cov_inc(int i) {

return ++i;

}

/*

void rc_update_tdata(BundleData& bundle, GffObj& scaff,

double cov, double fpkm) {

if (bundle.rc_data==NULL) return;

RC_BundleData& rc = *(bundle.rc_data);

if (rc.exons.size()==0) return;

RC_ScaffData& q = (RC_ScaffData*)scaff.uptr;

set<RC_ScaffData>::iterator tdata = rc.tdata.find(q);

if (tdata==rc.tdata.end()) {

fprintf(stderr, "Error: cannot locate bundle ref. transcript %s (%s:%d-%d)!\n",

scaff.getID(), scaff.getGSeqName(), scaff.start, scaff.end);

return;

}

(*tdata).cov=cov;

(*tdata).fpkm=fpkm;

}

*/

void BundleData::rc_store_t(GffObj* t) {

//if (!rc_stage) return;

if (rc_data==NULL) {

rc_init(t);

}

rc_data->addTranscript(*t);

//check this read alignment against ref exons and introns

}

struct COvlSorter {

bool operator() (pair<int, const RC_Feature*> i,

pair<int, const RC_Feature*> j) {

return (i.first>j.first); //sort in decreasing order of overlap length

}

} OvlSorter;

void rc_updateExonCounts(const RC_Feature* exon, int nh) {

exon->rcount++;

exon->mrcount += (nh > 1) ? (1.0/nh) : 1;

if (nh<=1) exon->ucount++;

}

bool BundleData::rc_count_hit(GBamRecord& brec, char strand, int nh) { //, int hi) {

if (rc_data==NULL) return false; //no ref transcripts available for this reads' region

if (rc_data->tdata.Count()==0) return false; //nothing to do without transcripts

//check this read alignment against ref exons and introns

/*

int gstart=brec.start; //alignment start position on the genome

We NEVER update read-counting boundaries of the bundle based on reads - they should only be based on reference transcripts

if (rc_data->f_cov.size()==0 && gstart<rc_data->lmin) {

fprintf(stderr, "Warning: adjusting lmin coverage bundle from %d to %d !\n", int(rc_data->lmin), (int)gstart);

rc_data->lmin=gstart;

}

*/

if ((int)brec.end<rc_data->lmin || (int)brec.start>rc_data->rmax) {

return false; //hit outside coverage area

}

/*

int gpos=brec.start; //current genomic position

int rlen=0; //read length, obtained here from the cigar string

int segstart=gstart;

*/

vector<RC_Seg> rsegs;

vector<RC_Seg> rintrons;

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

rc_data->updateCov(strand, nh, brec.exons[i].start, brec.exons[i].len());

rsegs.push_back(RC_Seg(brec.exons[i].start, brec.exons[i].end) );

if (i>0) {

//add intron

rintrons.push_back(RC_Seg(brec.exons[i-1].end+1, brec.exons[i].start-1));

}

}

//now check rexons and rintrons with findExons() and findIntron()

for (size_t i=0;i<rintrons.size();++i) {

/*RC_FeatIt ri=rc_data->findIntron(rintrons[i].l, rintrons[i].r, strand);

if (ri!=rc_data->introns.end()) {

(*ri).rcount++;

(*ri).mrcount += (nh > 1) ? (1.0/nh) : 1;

if (nh==1) (*ri).ucount++;

}

*/

RC_Feature* ri=rc_data->findIntron(rintrons[i].l, rintrons[i].r, strand);

if (ri) {

ri->rcount++;

ri->mrcount += (nh > 1) ? (1.0/nh) : 1;

if (nh==1) ri->ucount++;

}

} //for each intron

for (size_t i=0;i<rsegs.size();++i) {

RC_FeatPtrSet ovlex=rc_data->findExons(rsegs[i].l, rsegs[i].r, strand);

if (ovlex.size()==0) continue;

if (ovlex.size()>1) {

vector< pair<int, const RC_Feature*> > xovl; //overlapped exons sorted by decreasing overlap length

for (RC_FeatPtrSet::iterator ox=ovlex.begin();ox != ovlex.end(); ++ox) {

int ovlen=(*ox)->ovlen(rsegs[i].l, rsegs[i].r);

if (ovlen>=5)

xovl.push_back(pair<int, const RC_Feature*>(ovlen, *ox));

}

if (xovl.size()>1) {

sort(xovl.begin(), xovl.end(), OvlSorter); //larger overlaps first

//update the counts only for ref exons with max overlap to this segment

int max_ovl=xovl.begin()->first;

for (vector<pair<int, const RC_Feature*> >::iterator xo=xovl.begin();xo!=xovl.end();++xo) {

if (max_ovl - xo->first > 5 ) break; //more than +5 bases coverage for the other exons

rc_updateExonCounts(xo->second, nh);

}

} else if (xovl.size() == 1) {

rc_updateExonCounts(xovl.begin()->second, nh);

}

} else {

// 1 exon overlap only

int ovlen=(*ovlex.begin())->ovlen(rsegs[i].l, rsegs[i].r);

if (ovlen>=5) rc_updateExonCounts(*ovlex.begin(), nh);

}

} //for each read "exon"

return true;

}

FILE* rc_fwopen(const char* fname) {

if (strcmp(fname,"-")==0) return stdout;

GStr fpath(ballgown_dir); //always ends with '/'

//fpath += '.';

fpath += fname;

//fpath += "/";

//fpath += fname;

fpath += ".ctab";

FILE* fh=fopen(fpath.chars(), "w");

if (fh==NULL) {

fprintf(stderr, "Error: cannot create file %s\n",

fpath.chars());

exit(1);

}

return fh;

}

FILE* rc_frenopen(const char* fname) {

//if (strcmp(fname,"-")==0) return stdout;

GStr fpath(ballgown_dir);

//fpath += '.';

fpath += fname;

//fpath += "/";

//fpath += fname;

fpath += ".ctab";

GStr fren(fpath);

fren += ".tmp";

//rename fpath to fren and open fren for reading

if (rename(fpath.chars(), fren.chars())!=0) {

GError("Error: cannot rename %s to %s!\n", fpath.chars(), fren.chars());

}

FILE* fh=fopen(fren.chars(), "r");

if (fh==NULL) {

GError("Error: cannot open file %s\n", fren.chars());

}

return fh;

}

void rc_frendel(const char* fname) {

GStr fpath(ballgown_dir);

//fpath += '.';

fpath += fname;

fpath += ".ctab";

fpath += ".tmp";

if (remove(fpath.chars())!=0) {

GMessage("Warning: could not remove file %s!\n",fpath.chars());

}

}

void rc_write_f2t(FILE* fh, map<uint, set<uint> >& f2t) {

for (map<uint, set<uint> >::iterator m=f2t.begin(); m!=f2t.end(); ++m) {

uint f_id=(*m).first;

set<uint>& tset = (*m).second;

for (set<uint>::iterator it=tset.begin();it!=tset.end();++it) {

uint t_id = *it;

fprintf(fh, "%u\t%u\n", f_id, t_id);

}

}

fflush(fh);

}

void rc_update_exons(RC_BundleData& rc) {

//update stdev etc. for all exons in bundle

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

RC_Feature& exon = *(rc.exons[f]);

//assert( exon.l >= rc.lmin );

int L=exon.l-rc.lmin;

int xlen=exon.r-exon.l+1;

if (exon.l < rc.lmin) {

//shouldn't be here, bundle read-counting boundaries should be based on exons!

if (exon.r<rc.lmin) continue;

xlen-=(rc.lmin-exon.l+1);

L=0;

}

if (rc.rmax<exon.r) {

if (exon.l>rc.rmax) continue; //should never happen

xlen-=(exon.r-rc.rmax+1);

}

int R=L+xlen;

vector<int>::iterator xcov_begin;

vector<int>::iterator xcov_end;

vector<float>::iterator xmcov_begin;

vector<float>::iterator xmcov_end;

if (exon.strand=='+' || exon.strand=='.') {

xcov_begin = rc.f_cov.begin()+L;

xcov_end = rc.f_cov.begin()+R;

xmcov_begin = rc.f_mcov.begin()+L;

xmcov_end = rc.f_mcov.begin()+R;

} else {

xcov_begin = rc.r_cov.begin()+L;

xcov_end = rc.r_cov.begin()+R;

xmcov_begin = rc.r_mcov.begin()+L;

xmcov_end = rc.r_mcov.begin()+R;

}

double avg = (double)accumulate(xcov_begin, xcov_end, 0) / xlen;

vector<double> diff(xlen);

transform(xcov_begin, xcov_end, diff.begin(),

bind2nd( minus<double>(), avg));

double sq_sum = inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);

double stdev = sqrt(sq_sum / xlen);

double mavg = (double)accumulate(xmcov_begin, xmcov_end, 0) / xlen;

vector<double> mdiff(xlen);

transform(xmcov_begin, xmcov_end, mdiff.begin(),

bind2nd( minus<double>(), mavg));

sq_sum = inner_product(mdiff.begin(), mdiff.end(), mdiff.begin(), 0.0);

double mstdev = sqrt(sq_sum / xlen);

exon.avg=avg;

exon.stdev=stdev;

exon.mavg=mavg;

exon.mstdev=mstdev;

} //for each exon in bundle

}

void rc_write_RCfeature( GPVec<RC_ScaffData>& rcdata, GPVec<RC_Feature>& features, FILE*& fdata, FILE*& f2t,

bool is_exon=false) {

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

RC_Feature& f=*(features[i]);

const char* ref_name=rcdata[f.t_id-1]->scaff->getGSeqName();

if (is_exon) {

fprintf(fdata, "%u\t%s\t%c\t%d\t%d\t%d\t%d\t%.2f\t%.4f\t%.4f\t%.4f\t%.4f\n",

f.id, ref_name, f.strand, f.l, f.r, f.rcount,

f.ucount, f.mrcount, f.avg, f.stdev, f.mavg, f.mstdev);

}

else { //introns

fprintf(fdata,"%u\t%s\t%c\t%d\t%d\t%d\t%d\t%.2f\n",f.id, ref_name,

f.strand, f.l, f.r, f.rcount, f.ucount, f.mrcount);

}

// f2t -------

fprintf(f2t, "%u\t%u\n", f.id, f.t_id);

} //for each feature

fclose(fdata);

fclose(f2t);

}

/*void rc_write_counts(const char* refname, BundleData& bundle) {

RC_BundleData& rc = *bundle.rc_data;

//if (rc.exons.size()==0) return;

*

*/

void rc_writeRC(GPVec<RC_ScaffData>& RC_data,

GPVec<RC_Feature>& RC_exons,

GPVec<RC_Feature>& RC_introns,

FILE* &f_tdata, FILE* &f_edata, FILE* &f_idata,

FILE* &f_e2t, FILE* &f_i2t) {

for (int t=0;t<RC_data.Count();++t) {

//File: t_data.ctab

//t_id tname chr strand start end num_exons gene_id gene_name cufflinks_cov cufflinks_fpkm

const RC_ScaffData& sd=*RC_data[t];

const char* refname = sd.scaff->getGSeqName();

const char* genename= sd.scaff->getGeneName();

if (genename==NULL) genename=".";

fprintf(f_tdata, "%u\t%s\t%c\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t%f\t%f\n",

sd.t_id, refname, sd.strand, sd.l, sd.r, sd.scaff->getID(),

sd.num_exons, sd.eff_len, sd.scaff->getGeneID(),

genename, sd.cov, sd.fpkm);

}//for each transcript

//fflush(f_tdata);

fclose(f_tdata);

//File: e_data.ctab

//e_id chr gstart gend rcount ucount mrcount

rc_write_RCfeature(RC_data, RC_exons, f_edata, f_e2t, true);

//File: i_data.ctab

//i_id chr gstart gend rcount ucount mrcount

rc_write_RCfeature(RC_data, RC_introns, f_idata, f_i2t);

// feature-to-transcript link files

//rc_write_f2t(rc.fe2t, rc.e2t);

//rc_write_f2t(rc.fi2t, rc.i2t);

}

void RC_ScaffData::addFeature(int fl, int fr, GPVec<RC_Feature>& fvec,

uint& f_id, set<RC_ScaffSeg>& fset, set<RC_ScaffSeg>::iterator& fit,

GPVec<RC_Feature>& fdata) {

//f_id is the largest f_id inserted so far in fset

bool add_new = true;

RC_ScaffSeg newseg(fl,fr,this->strand);

//RC_Feature* newfeature=NULL;

if (fset.size()>0) {

if (fit == fset.end()) --fit;

if (newseg < (*fit)) {

bool eq=false;

while (newseg < (*fit) || (eq = (newseg==(*fit)))) {

if (eq) {

add_new = false;

newseg.id = fit->id;

break;

}

if (fit==fset.begin()) {

break;

}

--fit;

}

}

else { //newseg >= *fit

bool eq=false;

while ((*fit) < newseg || (eq = (newseg==(*fit)))) {

if (eq) {

add_new = false;

newseg.id = fit->id;

break;

}

++fit;

if (fit==fset.end()) {

--fit;

break;

}

}

}

}

if (add_new) {

newseg.id = ++f_id;

pair< set<RC_ScaffSeg>::iterator, bool> ret = fset.insert(newseg);

//ret.second = was_inserted (indeed new)

if (!ret.second) {

GError("Error: feature %d-%d (%c) already in segment set!\n",

newseg.l, newseg.r, newseg.strand);

//newseg.id = ret.first->id;

}

fdata.Add(new RC_Feature(newseg, this->t_id));

#ifdef DEBUG

if (fdata.Count()!=(int)f_id) {

GMessage("Error: fdata.Count=%d, f_id=%d\n", fdata.Count(), f_id);

}

#endif

GASSERT((uint)fdata.Count()==f_id);

}

//fvec.push_back(newseg);

GASSERT(fdata[newseg.id-1]->id==newseg.id);

fvec.Add(fdata[newseg.id-1]);

}

void Ballgown_setupFiles(FILE* &f_tdata, FILE* &f_edata, FILE* &f_idata,

FILE* &f_e2t, FILE* &f_i2t) {

if (f_tdata == NULL) {

//first call, create the files

f_tdata = rc_fwopen("t_data");

fprintf(f_tdata, "t_id\tchr\tstrand\tstart\tend\tt_name\tnum_exons\tlength\tgene_id\tgene_name\tcov\tFPKM\n");

f_edata = rc_fwopen("e_data");

fprintf(f_edata, "e_id\tchr\tstrand\tstart\tend\trcount\tucount\tmrcount\tcov\tcov_sd\tmcov\tmcov_sd\n");

f_idata = rc_fwopen("i_data");

fprintf(f_idata, "i_id\tchr\tstrand\tstart\tend\trcount\tucount\tmrcount\n");

f_e2t = rc_fwopen("e2t");

fprintf(f_e2t, "e_id\tt_id\n");

f_i2t = rc_fwopen("i2t");

fprintf(f_i2t, "i_id\tt_id\n");

}

}

void RC_ScaffData::rc_addFeatures(uint& c_e_id, set<RC_ScaffSeg>& fexons, GPVec<RC_Feature>& edata,

uint& c_i_id, set<RC_ScaffSeg>& fintrons, GPVec<RC_Feature>& idata) {

GASSERT(scaff);

GffObj& m = *(scaff);

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

set<RC_ScaffSeg>::iterator eit=fexons.end();

set<RC_ScaffSeg>::iterator iit=fintrons.end();

addFeature(m.exons[i]->start, m.exons[i]->end, t_exons, c_e_id, fexons, eit, edata);

if (i>0) { //store intron

addFeature(m.exons[i-1]->end+1, m.exons[i]->start-1, t_introns, c_i_id, fintrons, iit, idata);

}

} //for each exon

}