Skip to content

Commit 3b95d9d

Browse files
gperteagpertea
committed
merged 2.0.1 update
1 parent 8844944 commit 3b95d9d

5 files changed

Lines changed: 453 additions & 472 deletions

File tree

‎gclib/gff.h‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -651,11 +651,21 @@ class GffExon : public GSeg {
651651

652652
} //constructor
653653

654+
GffExon& operator=(const GffExon& o) { // copy operator : shallow!
655+
sharedAttrs=o.sharedAttrs;
656+
attrs=o.attrs;
657+
score=o.score;
658+
exontype=o.exontype;
659+
phase=o.phase;
660+
uptr=o.uptr;
661+
return *this;
662+
}
663+
654664
GffExon(const GffExon& ex):GSeg(ex.start, ex.end) { //copy constructor
655-
(*this)=ex; //use the default copy operator
656-
if (ex.attrs!=NULL) {
657-
attrs=new GffAttrs();
658-
attrs->copyAttrs(ex.attrs);
665+
(*this)=ex; //use the default (shallow!) copy operator
666+
if (ex.attrs!=NULL) { //make a deep copy here
667+
attrs=new GffAttrs();
668+
attrs->copyAttrs(ex.attrs);
659669
}
660670
}
661671
~GffExon() { //destructor

‎rlink.cpp‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -256,7 +256,7 @@ bool mismatch_anchor(CReadAln *rd,char *mdstr,int refstart, bam1_t *b) {
256256
rdlen=0;
257257
parsedlen=0;
258258
i=0;
259-
259+
260260
for (int j = 0; j < b->core.n_cigar; ++j) {
261261
int op = cigar[j]&0xf;
262262
if (op == BAM_CMATCH || op==BAM_CEQUAL ||
@@ -274,7 +274,7 @@ bool mismatch_anchor(CReadAln *rd,char *mdstr,int refstart, bam1_t *b) {
274274
}
275275
}
276276
}
277-
277+
278278
return false;
279279
}
280280

@@ -13286,7 +13286,9 @@ int print_predcluster(GList<CPrediction>& pred,int geneno,GStr& refname,
1328613286
}
1328713287

1328813288

13289-
GVec<int> bettercov[npred];
13289+
//GVec<int> bettercov[npred];
13290+
GVec<int> ivec;
13291+
GVec< GVec<int> > bettercov(npred, ivec);
1329013292

1329113293
predord.Sort(predordCmp); // sort predictions from the most highest coverage to lowest
1329213294
for(int i=0;i<predord.Count()-1;i++) {
@@ -13607,7 +13609,7 @@ int print_predcluster(GList<CPrediction>& pred,int geneno,GStr& refname,
1360713609
}
1360813610
}
1360913611

13610-
for(int n=0;n<npred;n++) if(pred[n]->flag){
13612+
for(int n=0;n<npred;n++) if(pred[n]->flag){
1361113613

1361213614
/*
1361313615
{ // DEBUG ONLY
@@ -13979,7 +13981,7 @@ uint find_midhash(int refstart,int start,int end,GVec<float>* bpcov) {
1397913981
*/
1398013982

1398113983

13982-
int printResults(BundleData* bundleData,int geneno, GStr& refname) {
13984+
int printResults(BundleData* bundleData, int geneno, GStr& refname) {
1398313985

1398413986
uint runoffdist=200;
1398513987
if(longreads) runoffdist=0;

‎stringtie.cpp‎

Lines changed: 39 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
#include "proc_mem.h"
1212
#endif
1313

14-
#define VERSION "2.01"
14+
#define VERSION "2.0.1"
1515

1616
//#define DEBUGPRINT 1
1717

@@ -22,11 +22,11 @@
2222
#define DBGPRINT4(a,b,c,d) GMessage(a,b,c,d)
2323
#define DBGPRINT5(a,b,c,d,e) GMessage(a,b,c,d,e)
2424
#else
25-
#define DBGPRINT(x)
26-
#define DBGPRINT2(a,b)
27-
#define DBGPRINT3(a,b,c)
28-
#define DBGPRINT4(a,b,c,d)
29-
#define DBGPRINT5(a,b,c,d,e)
25+
#define DBGPRINT(x)
26+
#define DBGPRINT2(a,b)
27+
#define DBGPRINT3(a,b,c)
28+
#define DBGPRINT4(a,b,c,d)
29+
#define DBGPRINT5(a,b,c,d,e)
3030
#endif
3131

3232
#define USAGE "StringTie v" VERSION " usage:\n\
@@ -251,9 +251,6 @@ bool noThreadsWaiting();
251251

252252
void workerThread(GThreadData& td); // Thread function
253253

254-
//check if a worker thread popped the data queue:
255-
bool queuePopped(GPVec<BundleData>& bundleQueue, int prevCount);
256-
257254
//prepare the next free bundle for loading
258255
int waitForData(BundleData* bundles);
259256
#endif
@@ -264,7 +261,7 @@ TInputFiles bamreader;
264261
int main(int argc, char* argv[]) {
265262

266263
// == Process arguments.
267-
GArgs args(argc, argv,
264+
GArgs args(argc, argv,
268265
"debug;help;version;conservative;keeptmp;bam;fr;rf;merge;exclude=zSEihvteuLx:n:j:s:D:G:C:l:m:o:a:j:c:f:p:g:P:M:Bb:A:F:T:");
269266
args.printError(USAGE, true);
270267

@@ -300,8 +297,9 @@ const char* ERR_BAM_SORT="\nError: the input alignment file is not sorted!\n";
300297
FILE* f=fopen(guidegff.chars(),"r");
301298
if (f==NULL) GError("Error: could not open reference annotation file (%s)!\n",
302299
guidegff.chars());
303-
// transcripts_only sort gffr->gfflst by loc?
304-
GffReader gffr(f, true, true); //loading only recognizable transcript features
300+
// transcripts_only sort by location?
301+
GffReader gffr(f, true, true); //loading only recognizable transcript features
302+
gffr.setRefAlphaSorted(); //alphabetical sorting of refseq IDs
305303
gffr.showWarnings(verbose);
306304
// keepAttrs mergeCloseExons noExonAttrs
307305
gffr.readAll(false, true, true);
@@ -337,8 +335,8 @@ const char* ERR_BAM_SORT="\nError: the input alignment file is not sorted!\n";
337335
if (skipGseq) continue;
338336
if (m->exons.Count()==0) {
339337
if (verbose)
340-
GMessage("Warning: exonless GFF object %s found, added implicit exon %d-%d.\n",
341-
m->getID(),m->start, m->end);
338+
GMessage("Warning: exonless GFF %s feature with ID %s found, added implicit exon %d-%d.\n",
339+
m->getFeatureName(), m->getID(), m->start, m->end);
342340
m->addExon(m->start, m->end); //should never happen!
343341
}
344342
//DONE: always keep a RC_TData pointer around, with additional info about guides
@@ -356,7 +354,6 @@ const char* ERR_BAM_SORT="\nError: the input alignment file is not sorted!\n";
356354
printTime(stderr);
357355
GMessage(" %d reference transcripts loaded.\n", gffr.gflst.Count());
358356
}
359-
//fclose(f); GffReader will close it anyway
360357
}
361358

362359
#ifdef GFF_DEBUG
@@ -398,6 +395,7 @@ const char* ERR_BAM_SORT="\nError: the input alignment file is not sorted!\n";
398395
if (ballgown)
399396
Ballgown_setupFiles(f_tdata, f_edata, f_idata, f_e2t, f_i2t);
400397
#ifndef NOTHREADS
398+
//model: one producer, multiple consumers
401399
#define DEF_TSTACK_SIZE 8388608
402400
int tstackSize=GThread::defaultStackSize();
403401
size_t defStackSize=0;
@@ -412,6 +410,8 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
412410
GThread* threads=new GThread[num_cpus]; //bundle processing threads
413411

414412
GPVec<BundleData> bundleQueue(false); //queue of loaded bundles
413+
//the consumers take (pop) bundles out of this queue for processing
414+
//the producer populates this queue with bundles, built from the BAM file
415415

416416
BundleData* bundles=new BundleData[num_cpus+1];
417417
//bundles[0..num_cpus-1] are processed by threads, loading bundles[num_cpus] first
@@ -436,7 +436,6 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
436436
bool chr_changed=false;
437437
int pos=0;
438438
const char* refseqName=NULL;
439-
//char strand=0;
440439
char xstrand=0;
441440
int nh=1;
442441
int hi=0;
@@ -554,12 +553,12 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
554553
}
555554
waitMutex.unlock();
556555
haveBundles.notify_one();
557-
this_thread::yield();
556+
current_thread::yield();
558557
queueMutex.lock();
559558
while (bundleQueue.Count()==qCount) {
560559
queueMutex.unlock();
561560
haveBundles.notify_one();
562-
this_thread::yield();
561+
current_thread::yield();
563562
queueMutex.lock();
564563
}
565564
queueMutex.unlock();
@@ -702,7 +701,7 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
702701
//countFragment(*bundle, *brec, hi,nh); // we count this in build_graphs to only include mapped fragments that we consider correctly mapped
703702
//fprintf(stderr,"fragno=%d fraglen=%lu\n",bundle->num_fragments,bundle->frag_len);if(bundle->num_fragments==100) exit(0);
704703
//if (!ballgown || ref_overlap)
705-
processRead(currentstart, currentend, *bundle, hashread, alndata);
704+
processRead(currentstart, currentend, *bundle, hashread, alndata);
706705
// *brec, strand, nh, hi);
707706
}
708707
} //for each read alignment
@@ -780,16 +779,12 @@ if(!mergeMode) {
780779
int nl;
781780
int istr;
782781
int tlen;
783-
float tcov;
784-
//float fpkm;
782+
float tcov; //do we need to increase precision here ? (double)
785783
float calc_fpkm;
786784
float calc_tpm;
787785
int t_id;
788786
while(fgetline(linebuf,linebuflen,ftmp_in)) {
789-
//sscanf(linebuf,"%d %d %d %g %g", &nl, &tlen, &t_id, &fpkm, &tcov);
790787
sscanf(linebuf,"%d %d %d %d %g", &istr, &nl, &tlen, &t_id, &tcov);
791-
//for the rare cases tcov < 0, invert it ??
792-
//if (tcov<0) tcov=-tcov; //should not happen
793788
if (tcov<0) tcov=0;
794789
calc_fpkm=tcov*1000000000/Frag_Len;
795790
calc_tpm=tcov*1000000/Cov_Sum;
@@ -912,6 +907,7 @@ void processOptions(GArgs& args) {
912907
args.printCmdLine(stderr);
913908
}
914909
//complete=!(args.getOpt('i')!=NULL);
910+
// trim=!(args.getOpt('t')!=NULL);
915911
includesource=!(args.getOpt('z')!=NULL);
916912
//EM=(args.getOpt('y')!=NULL);
917913
//weight=(args.getOpt('w')!=NULL);
@@ -1044,8 +1040,6 @@ void processOptions(GArgs& args) {
10441040
guidegff.chars());
10451041
}
10461042

1047-
1048-
10491043
enableNames=(args.getOpt('E')!=NULL);
10501044

10511045
retained_intron=(args.getOpt('i')!=NULL);
@@ -1222,11 +1216,11 @@ void noMoreBundles() {
12221216
if (areThreadsWaiting) {
12231217
DBGPRINT("##> NOTIFY ALL workers: no more data!\n");
12241218
haveBundles.notify_all();
1225-
this_thread::sleep_for(chrono::milliseconds(10));
1219+
current_thread::sleep_for(10);
12261220
waitMutex.lock();
12271221
areThreadsWaiting=(threadsWaiting>0);
12281222
waitMutex.unlock();
1229-
this_thread::sleep_for(chrono::milliseconds(10));
1223+
current_thread::sleep_for(10);
12301224
}
12311225
} while (areThreadsWaiting); //paranoid check that all threads stopped waiting
12321226
#else
@@ -1237,22 +1231,22 @@ void noMoreBundles() {
12371231
void processBundle(BundleData* bundle) {
12381232
if (verbose) {
12391233
#ifndef NOTHREADS
1240-
GLockGuard<GFastMutex> lock(logMutex);
1234+
GLockGuard<GFastMutex> lock(logMutex);
12411235
#endif
12421236
printTime(stderr);
1243-
GMessage(">bundle %s:%d-%d(%lu) (%d guides) loaded, begins processing...\n",
1244-
bundle->refseq.chars(), bundle->start, bundle->end, bundle->numreads,
1237+
GMessage(">bundle %s:%d-%d [%lu alignments (%d distinct), %d junctions, %d guides] begins processing...\n",
1238+
bundle->refseq.chars(), bundle->start, bundle->end, bundle->numreads, bundle->readlist.Count(), bundle->junction.Count(),
12451239
bundle->keepguides.Count());
1246-
#ifdef GMEMTRACE
1247-
double vm,rsm;
1248-
get_mem_usage(vm, rsm);
1249-
GMessage("\t\tstart memory usage: %6.1fMB\n",rsm/1024);
1250-
if (rsm>maxMemRS) {
1251-
maxMemRS=rsm;
1252-
maxMemVM=vm;
1253-
maxMemBundle.format("%s:%d-%d(%d)", bundle->refseq.chars(), bundle->start, bundle->end, bundle->readlist.Count());
1254-
}
1255-
#endif
1240+
#ifdef GMEMTRACE
1241+
double vm,rsm;
1242+
get_mem_usage(vm, rsm);
1243+
GMessage("\t\tstart memory usage: %6.1fMB\n",rsm/1024);
1244+
if (rsm>maxMemRS) {
1245+
maxMemRS=rsm;
1246+
maxMemVM=vm;
1247+
maxMemBundle.format("%s:%d-%d(%d)", bundle->refseq.chars(), bundle->start, bundle->end, bundle->readlist.Count());
1248+
}
1249+
#endif
12561250
}
12571251
#ifdef B_DEBUG
12581252
for (int i=0;i<bundle->keepguides.Count();++i) {
@@ -1318,8 +1312,8 @@ void processBundle(BundleData* bundle) {
13181312
fprintf(stderr,"Number of fragments in bundle: %g with sum %g\n",bundle->num_fragments,bundle->frag_len);
13191313
*/
13201314
printTime(stderr);
1321-
GMessage("^bundle %s:%d-%d(%d) done (%d processed potential transcripts).\n",bundle->refseq.chars(),
1322-
bundle->start, bundle->end, bundle->readlist.Count(), bundle->pred.Count());
1315+
GMessage("^bundle %s:%d-%d done (%d processed potential transcripts).\n",bundle->refseq.chars(),
1316+
bundle->start, bundle->end, bundle->pred.Count());
13231317
#ifdef GMEMTRACE
13241318
double vm,rsm;
13251319
get_mem_usage(vm, rsm);
@@ -1356,7 +1350,7 @@ void workerThread(GThreadData& td) {
13561350
queueMutex.unlock();
13571351
waitMutex.unlock();
13581352
haveThreads.notify_one(); //in case main thread is waiting
1359-
this_thread::yield();
1353+
current_thread::yield();
13601354
queueMutex.lock();
13611355
while (bundleWork && bundleQueue->Count()==0) {
13621356
haveBundles.wait(queueMutex);//unlocks queueMutex and wait until notified
@@ -1379,7 +1373,7 @@ void workerThread(GThreadData& td) {
13791373
dataClear.Push(readyBundle->idx);
13801374
dataMutex.unlock();
13811375
haveClear.notify_one(); //inform main thread
1382-
this_thread::yield();
1376+
current_thread::yield();
13831377
queueMutex.lock();
13841378
DBGPRINT2("---->> Thread%d locked back queueMutex\n", td.thread->get_id());
13851379

@@ -1390,15 +1384,6 @@ void workerThread(GThreadData& td) {
13901384
DBGPRINT2("---->> Thread%d DONE.\n", td.thread->get_id());
13911385
}
13921386

1393-
bool queuePopped(GPVec<BundleData>& bundleQueue, int prevCount) {
1394-
int c;
1395-
queueMutex.lock();
1396-
c=bundleQueue.Count();
1397-
queueMutex.unlock();
1398-
DBGPRINT3("##> post-notification check: qlen is now %d (was %d)\n", c, prevCount);
1399-
return (c==0 || c<prevCount);
1400-
}
1401-
14021387
//prepare the next available bundle slot for loading
14031388
int waitForData(BundleData* bundles) {
14041389
int bidx=-1;

0 commit comments

Comments
 (0)