-
Notifications
You must be signed in to change notification settings - Fork 80
/
Copy pathstringtie.cpp
1792 lines (1607 loc) · 58.4 KB
/
stringtie.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
//#define GFF_DEBUG 1 //debugging guides loading
#include "rlink.h"
#include "tmerge.h"
#ifndef NOTHREADS
#include "GThreads.h"
#endif
//#define GMEMTRACE 1
#ifdef GMEMTRACE
#include "proc_mem.h"
#endif
#define VERSION "3.0.1"
//#define DEBUGPRINT 1
#ifdef DEBUGPRINT
#define DBGPRINT(x) GMessage(x)
#define DBGPRINT2(a,b) GMessage(a,b)
#define DBGPRINT3(a,b,c) GMessage(a,b,c)
#define DBGPRINT4(a,b,c,d) GMessage(a,b,c,d)
#define DBGPRINT5(a,b,c,d,e) GMessage(a,b,c,d,e)
#else
#define DBGPRINT(x)
#define DBGPRINT2(a,b)
#define DBGPRINT3(a,b,c)
#define DBGPRINT4(a,b,c,d)
#define DBGPRINT5(a,b,c,d,e)
#endif
#define USAGE "StringTie v" VERSION " usage:\n\n\
stringtie <in.bam ..> [-G <guide_gff>] [-l <prefix>] [-o <out.gtf>] [-p <cpus>]\n\
[-v] [-a <min_anchor_len>] [-m <min_len>] [-j <min_anchor_cov>] [-f <min_iso>]\n\
[-c <min_bundle_cov>] [-g <bdist>] [-u] [-L] [-e] [--viral] [-E <err_margin>]\n\
[--ptf <f_tab>] [-x <seqid,..>] [-A <gene_abund.out>] [-h] {-B|-b <dir_path>}\n\
[--mix] [--conservative] [--rf] [--fr]\n\
Assemble RNA-Seq alignments into potential transcripts.\n\
Options:\n\
--version : print just the version at stdout and exit\n\
--conservative : conservative transcript assembly, same as -t -c 1.5 -f 0.05\n\
--mix : both short and long read data alignments are provided\n\
(long read alignments must be the 2nd BAM/CRAM input file)\n\
--rf : assume stranded library fr-firststrand\n\
--fr : assume stranded library fr-secondstrand\n\
-G reference annotation to use for guiding the assembly process (GTF/GFF)\n\
--ptf : load point-features from a given 4 column feature file <f_tab>\n\
-o output path/file name for the assembled transcripts GTF (default: stdout)\n\
-l name prefix for output transcripts (default: STRG)\n\
-f minimum isoform fraction (default: 0.01)\n\
-L long reads processing; also enforces -s 1.5 -g 0 (default:false)\n\
-R if long reads are provided, just clean and collapse the reads but\n\
do not assemble\n\
-m minimum assembled transcript length (default: 200)\n\
-a minimum anchor length for junctions (default: 10)\n\
-j minimum junction coverage (default: 1)\n\
-t disable trimming of predicted transcripts based on coverage\n\
(default: coverage trimming is enabled)\n\
-c minimum reads per bp coverage to consider for multi-exon transcript\n\
(default: 1)\n\
-s minimum reads per bp coverage to consider for single-exon transcript\n\
(default: 4.75)\n\
-v verbose (log bundle processing details)\n\
-g maximum gap allowed between read mappings (default: 50)\n\
-M fraction of bundle allowed to be covered by multi-hit reads (default:1)\n\
-p number of threads (CPUs) to use (default: 1)\n\
-A gene abundance estimation output file\n\
-N : use with non-polyA RNA-seq\n\
--nasc : output nascent transcripts\n\
-E define window around possibly erroneous splice sites from long reads to\n\
look out for correct splice sites (default: 25)\n\
-B enable output of Ballgown table files which will be created in the\n\
same directory as the output GTF (requires -G, -o recommended)\n\
-b enable output of Ballgown table files but these files will be \n\
created under the directory path given as <dir_path>\n\
-e only estimate the abundance of given reference transcripts (requires -G)\n\
--viral only relevant for long reads from viral data where splice sites\n\
do not follow consensus (default:false)\n\
-x do not assemble any transcripts on the given reference sequence(s)\n\
-u no multi-mapping correction (default: correction enabled)\n\
-h print this usage message and exit\n\
--ref/--cram-ref reference genome FASTA file for CRAM input\n\
\n\
Transcript merge usage mode: \n\
stringtie --merge [Options] { gtf_list | strg1.gtf ...}\n\
With this option StringTie will assemble transcripts from multiple\n\
input files generating a unified non-redundant set of isoforms. In this mode\n\
the following options are available:\n\
-G <guide_gff> reference annotation to include in the merging (GTF/GFF3)\n\
-o <out_gtf> output file name for the merged transcripts GTF\n\
(default: stdout)\n\
-m <min_len> minimum input transcript length to include in the merge\n\
(default: 50)\n\
-c <min_cov> minimum input transcript coverage to include in the merge\n\
(default: 0)\n\
-F <min_fpkm> minimum input transcript FPKM to include in the merge\n\
(default: 1.0)\n\
-T <min_tpm> minimum input transcript TPM to include in the merge\n\
(default: 1.0)\n\
-f <min_iso> minimum isoform fraction (default: 0.01)\n\
-g <gap_len> gap between transcripts to merge together (default: 250)\n\
-i keep merged transcripts with retained introns; by default\n\
these are not kept unless there is strong evidence for them\n\
-l <label> name prefix for output transcripts (default: MSTRG)\n\
"
/*
-C output a file with reference transcripts that are covered by reads\n\
-U unitigs are treated as reads and not as guides \n\ \\ not used now
-d disable adaptive read coverage mode (default: yes)\n\
-n sensitivity level: 0,1, or 2, 3, with 3 the most sensitive level (default 1)\n\ \\ deprecated for now
-O disable the coverage saturation limit and use a slower two-pass approach\n\
to process the input alignments, collapsing redundant reads\n\
-i the reference annotation contains partial transcripts\n\
-w weight the maximum flow algorithm towards the transcript with higher rate (abundance); default: no\n\
-y include EM algorithm in max flow estimation; default: no\n\
-z don't include source in the max flow algorithm\n\
-P output file with all transcripts in reference that are partially covered by reads
-M fraction of bundle allowed to be covered by multi-hit reads (paper uses default: 1)\n\
-c minimum bundle reads per bp coverage to consider for assembly (paper uses default: 3)\n\
-S more sensitive run (default: no) disabled for now \n\
-s coverage saturation threshold; further read alignments will be\n\ // this coverage saturation parameter is deprecated starting at version 1.0.5
ignored in a region where a local coverage depth of <maxcov> \n\
is reached (default: 1,000,000);\n\ \\ deprecated
-e (mergeMode) include estimated coverage information in the preidcted transcript\n\
-E (mergeMode) enable the name of the input transcripts to be included\n\
in the merge output (default: no)\n\
*/
//---- globals
FILE* f_out=NULL;
FILE* c_out=NULL;
//#define B_DEBUG 1
#ifdef B_DEBUG
FILE* dbg_out=NULL;
#endif
GStr outfname;
GStr out_dir;
GStr tmp_path;
GStr cram_ref; //reference genome FASTA for CRAM input
GStr tmpfname;
GStr genefname;
GStr traindir; // training directory for CDS option
bool guided=false;
bool trim=true;
bool viral=false;
bool eonly=false; // parameter -e ; for mergeMode includes estimated coverage sum in the merged transcripts
bool longreads=false;
bool rawreads=false;
bool nomulti=false;
bool enableNames=false;
bool includecov=false;
bool fr_strand=false;
bool rf_strand=false;
//bool complete=true; // set by parameter -i the reference annotation contains partial transcripts
bool retained_intron=false; // set by parameter -i for merge option
bool geneabundance=false;
//bool partialcov=false;
int num_cpus=1;
int mintranscriptlen=200; // minimum length for a transcript to be printed
//int sensitivitylevel=1;
uint junctionsupport=10; // anchor length for junction to be considered well supported <- consider shorter??
uint sserror=25; // window arround splice sites that we use to generate consensus in case of long read data
int junctionthr=1; // number of reads needed to support a particular junction
float readthr=1; // read coverage per bundle bp to accept it; // paper uses 3
float singlethr=4.75;
uint bundledist=50; // reads at what distance should be considered part of separate bundles
uint runoffdist=200;
float mcov=1; // fraction of bundle allowed to be covered by multi-hit reads paper uses 1
int allowed_nodes=1000;
//bool adaptive=true; // adaptive read coverage -> depends on the overall gene coverage
//GPVec<CDSparam> cds;
int no_xs=0; // number of records without the xs tag
float fpkm_thr=1;
float tpm_thr=1;
// different options of implementation reflected with the next three options
bool includesource=true;
//bool EM=false;
//bool weight=false;
float isofrac=0.01;
bool isunitig=true;
GStr label("STRG");
GStr ballgown_dir;
GFastaDb* gfasta=NULL;
GStr guidegff; // -G option
GStr ptff; // --ptf option (point features)
bool debugMode=false;
bool verbose=false;
bool ballgown=false;
bool havePtFeatures=false;
//int maxReadCov=1000000; //max local read coverage (changed with -s option)
//no more reads will be considered for a bundle if the local coverage exceeds this value
//(each exon is checked for this)
//bool forceBAM = false; //useful for stdin (piping alignments into StringTie)
bool mergeMode = false; //--merge option
bool keepTempFiles = false; //--keeptmp
bool genNascent = false; // -N/--nasc : internally generate synthetic nascent transcripts
bool printNascent = false; // output nascents
bool mixedMode = false; // both short and long read data alignments are provided
bool isnascent=false;
int GeneNo=0; //-- global "gene" counter
double Num_Fragments=0; //global fragment counter (aligned pairs)
double Frag_Len=0;
double Cov_Sum=0;
//bool firstPrint=true; //just for writing the GFF header before the first transcript is printed
GffNames* gseqNames=NULL; //used as a dictionary for reference sequence names and ids
int refseqCount=0; // number of reference sequences found in the guides file
#ifdef GMEMTRACE
double maxMemRS=0;
double maxMemVM=0;
GStr maxMemBundle;
#endif
#ifndef NOTHREADS
//single producer, multiple consumers
//main thread/program is always loading the producer
GMutex dataMutex; //manage availability of data records ready to be loaded by main thread
GVec<int> dataClear; //indexes of data bundles cleared for loading by main thread (clear data pool)
GConditionVar haveBundles; //will notify a thread that a bundle was loaded in the ready queue
//(or that no more bundles are coming)
int bundleWork=1; // bit 0 set if bundles are still being prepared (BAM file not exhausted yet)
// bit 1 set if there are Bundles ready in the queue
//GFastMutex waitMutex;
GMutex waitMutex; // controls threadsWaiting (idle threads counter)
int threadsWaiting; // idle worker threads
GConditionVar haveThreads; //will notify the bundle loader when a thread
//is available to process the currently loaded bundle
GConditionVar haveClear; //will notify when bundle buf space available
GMutex queueMutex; //controls bundleQueue and bundles access
GFastMutex printMutex; //for writing the output to file
GFastMutex logMutex; //only when verbose - to avoid mangling the log output
GFastMutex bamReadingMutex;
GFastMutex countMutex;
GFastMutex printCovMutex;
#endif
GStrSet<> excludeGseqs; //hash of chromosomes/contigs to exclude (e.g. chrM)
bool NoMoreBundles=false;
bool moreBundles(); //thread-safe retrieves NoMoreBundles
void noMoreBundles(); //sets NoMoreBundles to true
//--
void processOptions(GArgs& args);
int loadPtFeatures(FILE* f, GArray<GRefPtData>& refpts);
char* sprintTime();
void processBundle(BundleData* bundle);
//void processBundle1stPass(BundleData* bundle); //two-pass testing
void writeUnbundledGuides(GVec<GRefData>& refdata, FILE* fout, FILE* gout=NULL);
#ifndef NOTHREADS
bool noThreadsWaiting();
void workerThread(GThreadData& td); // Thread function
//prepare the next free bundle for loading
int waitForData(BundleData* bundles);
#endif
//#define DBG_ALN_DATA 1
#ifdef DBG_ALN_DATA
FILE* fdbgaln=NULL;
void dbg_waln(GSamRecord* b) {
if (fdbgaln==NULL) {
GStr fn=outfname;
if (fn.rindex('.')>0)
fn=fn.cut(fn.rindex('.'));
fn+="_alndbg.tab";
fdbgaln=fopen(fn.chars(), "w");
if (fdbgaln==NULL) GError("Error creating file %s\n", fn.chars());
}
// gseqname, flags, readname, start, end, cigar, nh, hi
char* pcigar=b->cigar();
fprintf(fdbgaln, "%s\t%d\t%s\t%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n", b->refName(), b->flags(),
b->name(), b->start, b->end, pcigar, (uint)b->tag_int("NM"), (uint)b->tag_int("NH"), (uint)b->tag_int("HI"),
b->mate_refId(), b->mate_start(), b->clipL, b->clipR);
GFREE(pcigar);
}
#endif
TInputFiles bamreader;
int main(int argc, char* argv[]) {
// == Process arguments.
GArgs args(argc, argv,
"debug;help;version;viral;conservative;mix;nasc;ref=;cram-ref=cds=;keeptmp;rseq=;ptf=;bam;fr;rf;merge;"
"exclude=zihvteuLRNx:n:j:s:D:G:C:S:l:m:o:a:j:c:f:p:g:P:M:Bb:A:E:F:T:");
args.printError(USAGE, true);
/**
* No code was selected, so no documentation can be generated.
*/
processOptions(args);
GVec<GRefData> refguides; // plain vector with guides loaded for each chromosome (if guided)
GArray<GRefPtData> refpts(true, true); // sorted,unique array of refseq point-features data
//table indexes for Ballgown Raw Counts data (-B/-b option)
GPVec<RC_TData> guides_RC_tdata(true); //raw count data or other info for all guide transcripts
GPVec<RC_Feature> guides_RC_exons(true); //raw count data for all guide exons
GPVec<RC_Feature> guides_RC_introns(true);//raw count data for all guide introns
GVec<int> alncounts(30); //keep track of the number of read alignments per chromosome [gseq_id]
int bamcount=bamreader.start(); //setup and open input files
//TODO: TEST ONLY!
//genNascent=true;
// ----
#ifndef GFF_DEBUG
if (bamcount<1) {
GError("%sError: no input files provided!\n",USAGE);
}
#endif
#ifdef DEBUGPRINT
verbose=true;
#endif
const char* ERR_BAM_SORT="\nError: the input alignment file is not sorted!\n";
if(guided) { // read guiding transcripts from input gff file
if (verbose) {
printTime(stderr);
GMessage(" Loading reference annotation (guides)..\n");
}
FILE* f=fopen(guidegff.chars(),"r");
if (f==NULL) GError("Error: could not open reference annotation file (%s)!\n",
guidegff.chars());
// transcripts_only sort by location?
GffReader gffr(f, true, true); //loading only recognizable transcript features
gffr.setRefAlphaSorted(); //alphabetical sorting of refseq IDs
gffr.showWarnings(verbose);
// keepAttrs mergeCloseExons noExonAttrs
gffr.readAll(false, true, true);
//the list of GffObj is in gffr.gflst, sorted by chromosome and start-end coordinates
//collect them in other data structures, if it's kept for later call gffobj->isUsed(true)
// (otherwise it'll be deallocated when gffr is destroyed due to going out of scope)
refseqCount=gffr.gseqtable.Count();
if (refseqCount==0 || gffr.gflst.Count()==0) {
GError("Error: could not any valid reference transcripts in %s (invalid GTF/GFF file?)\n",
guidegff.chars());
}
refguides.setCount(refseqCount); //maximum gseqid
uint c_tid=0; //counter for all transcripts
uint c_exon_id=0;
uint c_intron_id=0;
GList<RC_Feature> uexons(true, false, true); //sorted, free items, unique
GList<RC_Feature> uintrons(true, false, true);
//assign unique transcript IDs based on the sorted order
int last_refid=-1;
bool skipGseq=false;
for (int i=0;i<gffr.gflst.Count();i++) {
GffObj* m=gffr.gflst[i];
if (last_refid!=m->gseq_id) {
//chromosome switch
if (ballgown) { //prepare memory storage/tables for all guides on this chromosome
uexons.Clear();
uintrons.Clear();
}
last_refid=m->gseq_id;
skipGseq=excludeGseqs.hasKey(m->getGSeqName());
}
//sanity check: make sure there are no exonless "genes" or other
if (skipGseq) continue;
if (m->exons.Count()==0) {
if (verbose)
GMessage("Warning: exonless GFF %s feature with ID %s found, added implicit exon %d-%d.\n",
m->getFeatureName(), m->getID(), m->start, m->end);
m->addExon(m->start, m->end); //should never happen!
}
//DONE: always create a RC_TData for each guide in uptr, with additional tx data
RC_TData* tdata=new RC_TData(*m, ++c_tid);
m->uptr=tdata;
guides_RC_tdata.Add(tdata);
if (ballgown) { //already gather exon & intron info for all ref transcripts
tdata->rc_addFeatures(c_exon_id, uexons, guides_RC_exons,
c_intron_id, uintrons, guides_RC_introns);
}
GRefData& grefdata = refguides[m->gseq_id];
grefdata.add(&gffr, m); // add transcripts already sorted by location
}
if (verbose) {
printTime(stderr);
GMessage(" %d reference transcripts loaded.\n", gffr.gflst.Count());
}
}
gseqNames=GffObj::names; //might have been populated already by gff data
gffnames_ref(gseqNames); //initialize the names collection if not guided
// -- loading point-feature data
if (!ptff.is_empty()) {
FILE* f=fopen(ptff.chars(),"r");
if (f==NULL) GError("Error: could not open reference annotation file (%s)!\n",
ptff.chars());
// transcripts_only sort by location?
int numptf=loadPtFeatures(f, refpts); //adds to gseqNames->gseqs accordingly, populates refpts
havePtFeatures=(numptf>0);
fclose(f);
}
#ifdef GFF_DEBUG
for (int r=0;r<refguides.Count();++r) {
GRefData& grefdata = refguides[r];
for (int k=0;k<grefdata.rnas.Count();++k) {
GMessage("#transcript #%d : %s (%d exons)\n", k, grefdata.rnas[k]->getID(), grefdata.rnas[k]->exons.Count());
grefdata.rnas[k]->printGff(stderr);
}
}
GMessage("GFF Debug mode, exiting...\n");
exit(0);
#endif
// --- input processing
GHash<int> hashread; //read_name:pos:hit_index => readlist index
GList<GffObj>* guides=NULL; //list of transcripts on a specific reference
GList<GPtFeature>* refptfs=NULL; //list of point-features on a specific reference
int currentstart=0, currentend=0;
int ng_start=0;
int ng_end=-1;
int ptf_idx=0; //point-feature current index in the current (*refptfs)[]
int ng=0;
GStr lastref;
bool no_ref_used=true;
int lastref_id=-1; //last seen gseq_id
// int ncluster=0; used it for debug purposes only
Ref_RC_Data ref_rc; //will update ref_rc.refdata in loop
ref_rc.rc_tdata=&guides_RC_tdata;
ref_rc.rc_edata=&guides_RC_exons;
ref_rc.rc_idata=&guides_RC_introns;
//Ballgown files
FILE* f_tdata=NULL;
FILE* f_edata=NULL;
FILE* f_idata=NULL;
FILE* f_e2t=NULL;
FILE* f_i2t=NULL;
if (ballgown)
Ballgown_setupFiles(f_tdata, f_edata, f_idata, f_e2t, f_i2t);
#ifndef NOTHREADS
//model: one producer, multiple consumers
#define DEF_TSTACK_SIZE 8388608
size_t defStackSize=DEF_TSTACK_SIZE;
#ifdef _GTHREADS_POSIX_
size_t tstackSize=GThread::defaultStackSize();
if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
if (verbose) {
if (tstackSize<defStackSize){
size_t ssize=defStackSize;
GMessage("Default stack size for threads: %d (increased to %d)\n", tstackSize, ssize);
}
else GMessage("Default stack size for threads: %d\n", tstackSize);
}
#endif
GThread* threads=new GThread[num_cpus]; //bundle processing threads
GPVec<BundleData> bundleQueue(false); //queue of loaded bundles
//the consumers take (pop) bundles out of this queue for processing
//the producer populates this queue with bundles built from reading the BAM input
BundleData* bundles=new BundleData[num_cpus+1];
//bundles[0..num_cpus-1] are processed by threads, loading bundles[num_cpus] first
dataClear.setCapacity(num_cpus+1);
for (int b=0;b<num_cpus;b++) {
threads[b].kickStart(workerThread, (void*) &bundleQueue, defStackSize);
bundles[b+1].idx=b+1;
dataClear.Push(b);
}
BundleData* bundle = &(bundles[num_cpus]);
#else
BundleData bundles[1];
BundleData* bundle = &(bundles[0]);
#endif
GSamRecord* brec=NULL;
bool more_alns=true;
TAlnInfo* tinfo=NULL; // for --merge
int prev_pos=0;
bool skipGseq=false;
while (more_alns) {
bool chr_changed=false;
int pos=0;
const char* refseqName=NULL;
char xstrand=0;
int nh=1;
int hi=0;
int gseq_id=lastref_id; //current chr id
int bundle_last_kept_guide=-1;
bool new_bundle=false;
//delete brec;
if ((brec=bamreader.next())!=NULL) {
if (brec->isUnmapped()) continue;
if (brec->start<1 || brec->mapped_len<10) {
if (verbose) GMessage("Warning: invalid mapping found for read %s (position=%d, mapped length=%d)\n",
brec->name(), brec->start, brec->mapped_len);
continue;
}
#ifdef DBG_ALN_DATA
dbg_waln(brec);
#endif
refseqName=brec->refName();
xstrand=brec->spliceStrand(); // tagged strand gets priority
if(xstrand=='.' && (fr_strand || rf_strand)) { // set strand if stranded library
if(brec->isPaired()) { // read is paired
if(brec->pairOrder()==1) { // first read in pair
if((rf_strand && brec->revStrand())||(fr_strand && !brec->revStrand())) xstrand='+';
else xstrand='-';
}
else {
if((rf_strand && brec->revStrand())||(fr_strand && !brec->revStrand())) xstrand='-';
else xstrand='+';
}
}
else {
if((rf_strand && brec->revStrand())||(fr_strand && !brec->revStrand())) xstrand='+';
else xstrand='-';
}
}
/*
if (xstrand=='.' && brec->exons.Count()>1) {
no_xs++;
continue; //skip spliced alignments lacking XS tag (e.g. HISAT alignments)
}
// I might still infer strand later */
if (refseqName==NULL) GError("Error: cannot retrieve target seq name from BAM record!\n");
pos=brec->start; //BAM is 0 based, but GBamRecord makes it 1-based
chr_changed=(lastref.is_empty() || lastref!=refseqName);
if (chr_changed) {
skipGseq=excludeGseqs.hasKey(refseqName);
gseq_id=gseqNames->gseqs.addName(refseqName);
if (guided) {
if (gseq_id>=refseqCount) {
if (verbose)
GMessage("WARNING: no reference transcripts found for genomic sequence \"%s\"! (mismatched reference names?)\n",
refseqName);
}
else no_ref_used=false;
}
if (alncounts.Count()<=gseq_id) {
alncounts.Resize(gseq_id+1);
}
else if (alncounts[gseq_id]>0)
GError("%s\nAlignments (%d) already found for %s !\n",
ERR_BAM_SORT, alncounts[gseq_id], refseqName);
prev_pos=0;
}
if (pos<prev_pos) GError("%s\nread %s (start %d) found at position %d on %s when prev_pos=%d\n",
ERR_BAM_SORT, brec->name(), brec->start, pos, refseqName, prev_pos);
prev_pos=pos;
if (skipGseq) continue;
alncounts[gseq_id]++;
nh=brec->tag_int("NH");
if (nh==0) nh=1;
hi=brec->tag_int("HI");
if (mergeMode) {
//tinfo=new TAlnInfo(brec->name(), brec->tag_int("ZF"));
tinfo=new TAlnInfo(brec->name(), brec->uval);
GStr score(brec->tag_str("ZS"));
if (!score.is_empty()) {
GStr srest=score.split('|');
if (!score.is_empty())
tinfo->cov=score.asDouble();
score=srest.split('|');
if (!srest.is_empty())
tinfo->fpkm=srest.asDouble();
srest=score.split('|');
if (!score.is_empty())
tinfo->tpm=score.asDouble();
}
}
if (!chr_changed && currentend>0 && pos>currentend+(int)runoffdist) {
new_bundle=true;
}
}
else { //no more alignments
more_alns=false;
new_bundle=true; //fake a new start (end of last bundle)
}
if (new_bundle || chr_changed) { //read alignment starts a new bundle
hashread.Clear();
if (bundle->readlist.Count()>0) { // process reads in previous bundle
// (readthr, junctionthr, mintranscriptlen are globals)
if (refptfs) { //point-features defined for this reference
while (ptf_idx<refptfs->Count() && (int)(refptfs->Get(ptf_idx)->coord)<currentstart)
ptf_idx++;
//TODO: what if a PtFeature is nearby, just outside the bundle?
while (ptf_idx<refptfs->Count() && (int)(refptfs->Get(ptf_idx)->coord)<=currentend) {
bundle->ptfs.Add(refptfs->Get(ptf_idx)); //keep this PtFeature
ptf_idx++;
}
}
bundle->getReady(currentstart, currentend);
if (gfasta!=NULL) { //genomic sequence data requested
GFaSeqGet* faseq=gfasta->fetch(bundle->refseq.chars());
if (faseq==NULL) {
GError("Error: could not retrieve sequence data for %s!\n", bundle->refseq.chars());
}
bundle->gseq=faseq->copyRange(bundle->start, bundle->end, false, true);
}
#ifndef NOTHREADS
//push this in the bundle queue where it'll be picked up by the threads
DBGPRINT2("##> Locking queueMutex to push loaded bundle into the queue (bundle.start=%d)\n", bundle->start);
int qCount=0;
queueMutex.lock();
bundleQueue.Push(bundle);
bundleWork |= 0x02; //set bit 1
qCount=bundleQueue.Count();
queueMutex.unlock();
DBGPRINT2("##> bundleQueue.Count()=%d)\n", qCount);
//wait for a thread to pop this bundle from the queue
waitMutex.lock();
DBGPRINT("##> waiting for a thread to become available..\n");
while (threadsWaiting==0) {
haveThreads.wait(waitMutex);
}
waitMutex.unlock();
haveBundles.notify_one();
DBGPRINT("##> waitMutex unlocked, haveBundles notified, current thread yielding\n");
current_thread::yield();
queueMutex.lock();
DBGPRINT("##> queueMutex locked until bundleQueue.Count()==qCount\n");
while (bundleQueue.Count()==qCount) {
queueMutex.unlock();
DBGPRINT2("##> queueMutex unlocked as bundleQueue.Count()==%d\n", qCount);
haveBundles.notify_one();
current_thread::yield();
queueMutex.lock();
DBGPRINT("##> queueMutex locked again within while loop\n");
}
queueMutex.unlock();
#else //no threads
//Num_Fragments+=bundle->num_fragments;
//Frag_Len+=bundle->frag_len;
processBundle(bundle);
#endif
// ncluster++; used it for debug purposes only
} //have alignments to process
else { //no read alignments in this bundle?
#ifndef NOTHREADS
dataMutex.lock();
DBGPRINT2("##> dataMutex locked for bundle #%d clearing..\n", bundle->idx);
#endif
bundle->Clear();
#ifndef NOTHREADS
dataClear.Push(bundle->idx);
DBGPRINT2("##> dataMutex unlocking as dataClear got pushed idx #%d\n", bundle->idx);
dataMutex.unlock();
#endif
} //nothing to do with this bundle
if (chr_changed) {
if (guided) {
ng=0;
guides=NULL;
ng_start=0;
ng_end=-1;
if (refguides.Count()>gseq_id && refguides[gseq_id].rnas.Count()>0) {
guides=&(refguides[gseq_id].rnas);
ng=guides->Count();
ref_rc.refdata=&refguides[gseq_id];
}
}
if (havePtFeatures) {
ptf_idx=-1;
//setup refptf
refptfs=NULL;
GRefPtData rd(gseq_id);
int ridx=refpts.IndexOf(rd);
if (ridx>=0) {
refptfs=&(refpts[ridx].pfs);
ptf_idx=0;
}
}
lastref=refseqName;
lastref_id=gseq_id;
currentend=0;
}
if (!more_alns) {
if (verbose) {
#ifndef NOTHREADS
GLockGuard<GFastMutex> lock(logMutex);
#endif
if (Num_Fragments) {
printTime(stderr);
GMessage(" %g aligned fragments found.\n", Num_Fragments);
}
//GMessage(" Done reading alignments.\n");
}
noMoreBundles();
break;
}
#ifndef NOTHREADS
int new_bidx=waitForData(bundles);
if (new_bidx<0) {
//should never happen!
GError("Error: waitForData() returned invalid bundle index(%d)!\n",new_bidx);
break;
}
bundle=&(bundles[new_bidx]);
#endif
currentstart=pos;
currentend=brec->end;
if (guides) { //guided and guides!=NULL
bundle_last_kept_guide=bundle->keepguides.Count();
ng_start=ng_end+1;
while (ng_start<ng && (int)(*guides)[ng_start]->end < pos) {
// for now, skip guides which have no overlap with current read
ng_start++;
}
int ng_ovl=ng_start;
//add all guides w/ overlap with current read, + other guides that overlap them, transitively
while (ng_ovl<ng && (int)(*guides)[ng_ovl]->start<=currentend) { //while guide overlap
if (currentstart>(int)(*guides)[ng_ovl]->start)
currentstart=(*guides)[ng_ovl]->start;
if (currentend<(int)(*guides)[ng_ovl]->end)
currentend=(*guides)[ng_ovl]->end;
if (ng_ovl==ng_start && ng_ovl>0) {
// -- first time only, we have to check back all possible transitive guide overlaps
//char* geneid=(*guides)[ng_ovlstart]->getGeneID();
//if (geneid==NULL) geneid=(*guides)[ng_ovlstart]->getGeneName();
//if (geneid && !bgeneids.hasKey(geneid))
// bgeneids.shkAdd(geneid, &ng);
int g_back=ng_ovl; //start from the overlapping guide, going backwards
int g_ovl_start=ng_ovl;
while (g_back>ng_end+1) {
--g_back;
//if overlap, set g_back_start=g_back and update currentstart
if (currentstart<=(int)(*guides)[g_back]->end) {
g_ovl_start=g_back;
if (currentstart>(int)(*guides)[g_back]->start)
currentstart=(int)(*guides)[g_back]->start;
}
} //while checking previous guides that could be pulled in this bundle
//adding all guides that overlap the first added guide, in order
for (int gb=g_ovl_start;gb<=ng_ovl;++gb) {
bundle->keepGuide((*guides)[gb], ref_rc);
}
} //for 1st guide overlap: needed to check previous guides for overlaps
else // for the rest, just add the guide to the bundle
bundle->keepGuide((*guides)[ng_ovl], ref_rc);
ng_ovl++;
} //while guide overlap
ng_end=ng_ovl-1; //MUST update ng_end here, even if no overlaps were found
} //guides present on the current chromosome
bundle->refseq=lastref;
bundle->start=currentstart;
bundle->end=currentend;
} //<---- read alignment started a new bundle
if (currentend<(int)brec->end) { //read alignments extends the current bundle
//this might not happen if a longer guide had already been added to the bundle
currentend=brec->end;
if (guides) { //add any newly overlapping guides to bundle
if (bundle_last_kept_guide<0) bundle_last_kept_guide=bundle->keepguides.Count();
bool cend_changed;
do {
cend_changed=false;
while (ng_end+1<ng && (int)(*guides)[ng_end+1]->start<=currentend) {
++ng_end;
//more transcripts overlapping this bundle?
if ((int)(*guides)[ng_end]->end>=currentstart) {
//no transitive overlap check here, guide must overlap the bundle
bundle->keepGuide((*guides)[ng_end], ref_rc);
if(currentend<(int)(*guides)[ng_end]->end) {
currentend=(*guides)[ng_end]->end;
cend_changed=true;
}
}
}
} while (cend_changed);
}
} //adjusted currentend and checked for overlapping reference transcripts
if (guides && genNascent && bundle_last_kept_guide>=0)
bundle->generateAllNascents(bundle_last_kept_guide, ref_rc);
GReadAlnData alndata(brec, 0, nh, hi, tinfo);
bool ovlpguide=bundle->evalReadAln(alndata, xstrand);
if(!eonly || ovlpguide) { // in eonly case consider read only if it overlaps guide
//check for overlaps with ref transcripts which may set xstrand
if (xstrand=='+') alndata.strand=1;
else if (xstrand=='-') alndata.strand=-1;
//const char* bname=brec->name();
//GMessage("%s\t%c\t%d\thi=%d\n",bname, xstrand, alndata.strand,hi);
//countFragment(*bundle, *brec, hi,nh); // we count this in build_graphs to only include mapped fragments that we consider correctly mapped
//fprintf(stderr,"fragno=%d fraglen=%lu\n",bundle->num_fragments,bundle->frag_len);if(bundle->num_fragments==100) exit(0);
processRead(currentstart, currentend, *bundle, hashread, alndata,ovlpguide);
}
} //for each read alignment
//cleaning up
delete brec;
bamreader.stop(); //close all BAM files
if (guided && no_ref_used) {
GMessage("WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!\n"
"Please make sure the -G annotation file uses the same naming convention for the genome sequences.\n");
}
delete gfasta;
#ifndef NOTHREADS
for (int t=0;t<num_cpus;t++)
threads[t].join();
if (verbose) {
printTime(stderr);
GMessage(" All threads finished.\n");
}
delete[] threads;
delete[] bundles;
#else
if (verbose) {
printTime(stderr);
GMessage(" Done.\n");
}
#endif
#ifdef DBG_ALN_DATA
fclose(fdbgaln);
#endif
#ifdef B_DEBUG
fclose(dbg_out);
#endif
if (mergeMode && guided )
writeUnbundledGuides(refguides, f_out);
// clear refpts data, if loaded
if (refpts.Count()>0)
for (int i=0;i<refpts.Count();i++) {
refpts[i].pfs.setFreeItem(true);
}
fclose(f_out);
if (c_out && c_out!=stdout) fclose(c_out);
if(verbose && no_xs>0)
GMessage("Number spliced alignments missing the XS tag (skipped): %d\n",no_xs);
if(!mergeMode) {
if(verbose) {
GMessage("Total count of aligned fragments: %g\n", Num_Fragments);
if (Num_Fragments)
GMessage("Fragment coverage length: %g\n", Frag_Len/Num_Fragments);
}
f_out=stdout;
if(outfname!="stdout") {
f_out=fopen(outfname.chars(), "w");
if (f_out==NULL) GError("Error creating output file %s\n", outfname.chars());
}
fprintf(f_out,"# ");
args.printCmdLine(f_out);
fprintf(f_out,"# StringTie version %s\n",VERSION);
//fprintf(stderr,"cov_sum=%g frag_len=%g num_frag=%g\n",Cov_Sum,Frag_Len,Num_Fragments);
FILE *g_out=NULL;
if(geneabundance) {
g_out=fopen(genefname.chars(),"w");
if (g_out==NULL)
GError("Error creating gene abundance output file %s\n", genefname.chars());
fprintf(g_out,"Gene ID\tGene Name\tReference\tStrand\tStart\tEnd\tCoverage\tFPKM\tTPM\n");
}
FILE* ftmp_in=fopen(tmpfname.chars(),"rt");
if (ftmp_in!=NULL) {
char* linebuf=NULL;
int linebuflen=5000;
GMALLOC(linebuf, linebuflen);
int nl;
int istr;
int tlen;
float tcov; //do we need to increase precision here ? (double)
float calc_fpkm;
float calc_tpm;
int t_id;
while(fgetline(linebuf,linebuflen,ftmp_in)) {
sscanf(linebuf,"%d %d %d %d %g", &istr, &nl, &tlen, &t_id, &tcov);
if (tcov<0) tcov=0;
if (Frag_Len>0.001) calc_fpkm=tcov*1000000000/Frag_Len;
else calc_fpkm=0.0;
if (Cov_Sum>0.00001) calc_tpm=tcov*1000000/Cov_Sum;
else calc_tpm=0.0;
if(istr) { // this is a transcript
if (ballgown && t_id>0) {
guides_RC_tdata[t_id-1]->fpkm=calc_fpkm;
guides_RC_tdata[t_id-1]->cov=tcov;
}
for(int i=0;i<nl;i++) {
fgetline(linebuf,linebuflen,ftmp_in);
if(!i) {
//linebuf[strlen(line)-1]='\0';
fprintf(f_out,"%s",linebuf);
fprintf(f_out," FPKM \"%.6f\";",calc_fpkm);
fprintf(f_out," TPM \"%.6f\";",calc_tpm);
fprintf(f_out,"\n");
}
else fprintf(f_out,"%s\n",linebuf);
}
}
else { // this is a gene -> different file pointer
fgetline(linebuf, linebuflen, ftmp_in);
fprintf(g_out, "%s\t%.6f\t%.6f\n", linebuf, calc_fpkm, calc_tpm);
}
}
if (guided) {
writeUnbundledGuides(refguides, f_out, g_out);
}
fclose(f_out);
fclose(ftmp_in);
if(geneabundance) fclose(g_out);
GFREE(linebuf);
if (!keepTempFiles) {
int rv=remove(tmpfname.chars());
if (rv) {
perror(tmpfname.chars());
}
}
}
else {
fclose(f_out);
GError("No temporary file %s present!\n",tmpfname.chars());
}
//lastly, for ballgown, rewrite the tdata file with updated cov and fpkm
if (ballgown) {
rc_writeRC(guides_RC_tdata, guides_RC_exons, guides_RC_introns,
f_tdata, f_edata, f_idata, f_e2t, f_i2t);
}
}
if (!keepTempFiles) {
tmp_path.chomp('/');
// check the results of remove() and show error if it fails
if (Grmdir(tmp_path.chars())) {
GMessage("Error removing temporary directory: %s\n", tmp_path.chars());
}
}
gffnames_unref(gseqNames); //deallocate names collection
#ifdef GMEMTRACE
if(verbose) GMessage(" Max bundle memory: %6.1fMB for bundle %s\n", maxMemRS/1024, maxMemBundle.chars());