← Back to branch summary

~ubuntu-branches/ubuntu/trusty/mira/trusty-proposed

« back to all changes in this revision

Viewing changes to src/mira/readpool.C

  • Committer: Package Import Robot
  • Author(s): Thorsten Alteholz, Thorsten Alteholz, Andreas Tille
  • Date: 2014-02-02 22:51:35 UTC
  • mfrom: (7.1.1 sid)
  • Revision ID: package-import@ubuntu.com-20140202225135-nesemzj59jjgogh0
Tags: 4.0-1
http://bugs.debian.org/735798
[ Thorsten Alteholz ]
* New upstream version 
* debian/rules: add boost dir in auto_configure (Closes: #735798)

[ Andreas Tille ]
* cme fix dpkg-control
* debian/patches/{make.patch,spelling.patch}: applied upstream (thus removed)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mira/readpool.C
26
26
 */
27
27
 
28
28
 
29
 
#include "readpool.H"
30
 
 
31
 
#include <iostream>
32
 
#include <fstream>
33
 
#include <ctype.h>
 
29
#include "mira/readpool.H"
34
30
 
35
31
#include "errorhandling/errorhandling.H"
36
 
#include "util/fileanddisk.H"
37
 
#include "util/progressindic.H"
38
32
 
39
33
#include "io/generalio.H"
40
34
#include "io/fasta.H"
42
36
#include "io/phd.H"
43
37
#include "io/ncbiinfoxml.H"
44
38
 
 
39
#include "util/fileanddisk.H"
 
40
#include "util/progressindic.H"
 
41
 
 
42
#include "caf/caf.H"
 
43
 
45
44
#include "mira/gbf_parse.H"
46
45
#include "mira/gff_parse.H"
47
46
#include "mira/maf_parse.H"
48
 
#include "caf/caf.H"
49
47
 
50
48
#include <boost/unordered_map.hpp>
51
49
#include <boost/unordered_set.hpp>
252
250
  int32 acttid=0;
253
251
  int32 validreads=0;
254
252
 
255
 
  uint32 outlines=0;
 
253
  uint32 outlines_fatal=0;
 
254
  uint32 outlines_warn=0;
 
255
  uint32 ol_emptytn=0;
 
256
  uint32 ol_gt2reads=0;
 
257
  uint32 ol_unknownsegment=0;
 
258
  uint32 ol_unknownsegmentfixed=0;
 
259
  uint32 ol_samesegment=0;
256
260
  for(uint32 rid=0; rid<size();++rid){
257
261
      //cout << "acttid: " << acttid << "\t";
258
262
    if(getRead(rid).hasValidData()==false) continue;
263
267
      ++acttid;
264
268
      if(getRead(rid).getReadGroupID().hasTemplateInfo()){
265
269
        // Ooooops ... empty template but readgroupinfo?
266
 
        if(outlines<200){
 
270
        if(ol_emptytn<200){
 
271
          ++outlines_fatal;
267
272
          cout << "Read " << getRead(rid).getName() << " has template info given, but the internal template name is empty? This is fishy!" << endl;
268
 
          if(++outlines==200){
 
273
          if(++ol_emptytn==200){
269
274
            cout << "More than 200 cases like the above, will not report more.\n";
270
275
          }
271
276
        }
278
283
        if(++tidcounter[tnI->second]==2){
279
284
          getRead(rid).setTemplatePartnerID(firstpartner);
280
285
          getRead(firstpartner).setTemplatePartnerID(rid);
 
286
 
 
287
          // look out for unknown segments ("0"), try to fix if possible, else dump error message
 
288
          uint8 numunknownseg=(getRead(rid).getTemplateSegment()==0)+(getRead(firstpartner).getTemplateSegment()==0);
 
289
          if(numunknownseg){
 
290
            bool couldfix=false;
 
291
            if(numunknownseg==1){
 
292
              Read * tofix=&getRead(rid);
 
293
              Read * other=&getRead(firstpartner);
 
294
              uint8 segmenttoset=0;
 
295
              if(tofix->getTemplateSegment()>0) swap(tofix,other);
 
296
              if(other->getTemplateSegment()==1){
 
297
                segmenttoset=255;
 
298
              }else if(other->getTemplateSegment()==255){
 
299
                segmenttoset=1;
 
300
              }
 
301
              if(segmenttoset){
 
302
                couldfix=true;
 
303
                tofix->setTemplateSegment(segmenttoset);
 
304
                ol_unknownsegmentfixed++;
 
305
                if(ol_unknownsegmentfixed<200){
 
306
                  cout << "Read " << tofix->getName() << ": fixed unrecognised template segment.\n";
 
307
                }
 
308
                if(++ol_unknownsegmentfixed==200){
 
309
                  cout << "More than 200 cases like the above, will not report more.\n";
 
310
                }
 
311
              }
 
312
            }
 
313
            if(!couldfix){
 
314
              if(ol_unknownsegment<200){
 
315
                ++outlines_fatal;
 
316
                if(getRead(rid).getTemplateSegment()==0){
 
317
                  cout << "DNA template " << tnI->first << ", read " << getRead(rid).getName() << ": template segment not recognised.\n";
 
318
                }
 
319
                if(getRead(firstpartner).getTemplateSegment()==0){
 
320
                  cout << "DNA template " << tnI->first << ", read " << getRead(firstpartner).getName() << ": template segment not recognised.\n";
 
321
                }
 
322
                if(++ol_unknownsegment==200){
 
323
                  cout << "More than 200 cases like the above, will not report more.\n";
 
324
                }
 
325
              }
 
326
            }
 
327
          }
 
328
 
 
329
          // the folowing is almost impossible when reading FASTA/FASTQ and other simple files,
 
330
          //  but it may happen in CAF/MAF (and maybe SAM?)
 
331
          if(getRead(rid).getTemplateSegment()==getRead(firstpartner).getTemplateSegment()){
 
332
            if(ol_samesegment<200){
 
333
              ++outlines_fatal;
 
334
              cout << "Reads " << getRead(rid).getName() << " and " << getRead(firstpartner).getName() << " have the same template segments.\n";
 
335
              if(++ol_samesegment==200){
 
336
                cout << "More than 200 cases like the above, will not report more.\n";
 
337
              }
 
338
            }
 
339
          }
 
340
 
281
341
        }else{
282
 
          // Ooooops ... not good: more than one read for this template
 
342
          // Ooooops ... not good: more than two reads for this template
283
343
          getRead(rid).setTemplatePartnerID(-1);
284
344
          getRead(firstpartner).setTemplatePartnerID(-1);
285
 
          if(outlines<200){
 
345
          if(ol_gt2reads<200){
 
346
            ++outlines_fatal;
286
347
            cout << tidcounter[tnI->second] << " ";
287
348
            cout << "DNA template " << tnI->first << " has more than two reads, template info not used. Read at fault: " << getRead(rid).getName() << endl;
288
 
            if(++outlines==200){
 
349
            if(++ol_gt2reads==200){
289
350
              cout << "More than 200 cases like the above, will not report more.\n";
290
351
            }
291
352
          }
302
363
    }
303
364
  }
304
365
 
305
 
  if(outlines>0 && (*REP_miraparams)[0].getNagAndWarnParams().nw_stop_templateproblems){
306
 
    MIRANOTIFY(Notify::FATAL,"Problems found with the template data, see output log for more info. This points some serious problem either with read naming (like unrecognised read naming scheme) or broken template information, please fix your input files!\nOr switch off this warning with -NW:sote=no (but you'll do this at own risk!)")
 
366
  if(outlines_fatal>0 && (*REP_miraparams)[0].getNagAndWarnParams().nw_check_templateproblems){
 
367
    static string emsg="Problems found with the template data, see output log for more info. This points some serious problem either with read naming (like unrecognised read naming scheme) or broken template information, please fix your input files!\nOr switch off this warning with -NW:ctp=no (but you'll do this at own risk!)";
 
368
    if((*REP_miraparams)[0].getNagAndWarnParams().nw_check_templateproblems==NWSTOP) {
 
369
      MIRANOTIFY(Notify::FATAL,emsg);
 
370
    }else{
 
371
      cout << "WARNING!\n" << emsg << endl;
 
372
    }
307
373
  }
308
374
 
309
375
  if(verbose){
315
381
}
316
382
 
317
383
 
 
384
/*************************************************************************
 
385
 *
 
386
 *
 
387
 *
 
388
 *************************************************************************/
 
389
 
 
390
void ReadPool::checkTemplateIDs(const string & errmsg)
 
391
{
 
392
  FUNCSTART("void ReadPool::checkTemplateIDs(string & errmsg)");
 
393
  for(uint32 ri=0;ri< REP_thepool3.size(); ++ri){
 
394
    if(REP_thepool3.getRead(ri).getTemplatePartnerID()>=0
 
395
       && REP_thepool3.getRead(ri).getTemplateID() != REP_thepool3.getRead(REP_thepool3.getRead(ri).getTemplatePartnerID()).getTemplateID()){
 
396
      cout << "Ouch, template problem for read " << ri << " " << REP_thepool3.getRead(ri).getName() << ", dumping readpool for debug\n";
 
397
      dumpAs(cout,Read::AS_TEXT,true);
 
398
      BUGIFTHROW(true, errmsg);
 
399
    }
 
400
  }
 
401
}
318
402
 
319
403
/*************************************************************************
320
404
 *
349
433
    }
350
434
  }
351
435
 
352
 
  if(outlines>0 && (*REP_miraparams)[0].getNagAndWarnParams().nw_stop_duplicatereadnames){
353
 
    MIRANOTIFY(Notify::FATAL,"Some read names were found more than once (see log above). This usually hints to a serious problem with your input and should really, really be fixed. You can choose to ignore this error with '-NW:sodrn=no', but this will almost certainly lead to problems with result files (ACE and CAF for sure, maybe also SAM) and probably to other unexpected effects.");
 
436
  if(outlines>0 && (*REP_miraparams)[0].getNagAndWarnParams().nw_check_duplicatereadnames){
 
437
    static string emsg="Some read names were found more than once (see log above). This usually hints to a serious problem with your input and should really, really be fixed. You can choose to ignore this error with '-NW:cdrn=no', but this will almost certainly lead to problems with result files (ACE and CAF for sure, maybe also SAM) and probably to other unexpected effects.";
 
438
    if((*REP_miraparams)[0].getNagAndWarnParams().nw_check_duplicatereadnames==NWSTOP){
 
439
      MIRANOTIFY(Notify::FATAL,emsg);
 
440
    }else{
 
441
      cout << "WARNING!\n" << emsg << endl;
 
442
    }
354
443
  }
355
444
 
356
 
 
357
445
  FUNCEND();
358
446
  return allok;
359
447
}
383
471
    loadDataFromFASTQ_rgid(filename1, fqqualoffset, rgid, countonly, callback);
384
472
  }else if(filetype=="fasta"){
385
473
    loadDataFromFASTA_rgid(filename1, rgid, true, optfilename2, countonly, callback);
386
 
  }else if(filetype=="fna" || filetype == "fastanoqual"){ // FASTA, but no qual. fna to be used by users, fastanoqual internally
 
474
  }else if(filetype=="fna" || filetype == "fastanoqual" || filetype=="fa"){ // FASTA, but no qual. fna to be used by users, fastanoqual internally
387
475
    loadDataFromFASTA_rgid(filename1, rgid, false, "", countonly, callback);
388
476
  }else if(filetype=="gbf"){
389
477
    loadDataFromGBF_rgid(filename1, rgid, countonly, callback);
390
478
  }else if(filetype=="gff3"){
391
479
    loadDataFromGFF3_rgid(filename1, rgid, countonly, callback);
392
480
  }else if(filetype=="fofnexp"){
393
 
    loadDataFromEXPs_rgid(filename1, rgid, countonly, callback);
 
481
    loadDataFromFOFNEXP_rgid(filename1, rgid, countonly, callback);
 
482
  }else if(filetype=="exp"){
 
483
    loadDataFromEXP_rgid(filename1, rgid, countonly, callback);
394
484
  }else if(filetype=="caf"){
395
485
    loadDataFromCAF_rgid(filename1, rgid, countonly, callback);
396
486
  }else if(filetype=="maf"){
420
510
 *
421
511
 *************************************************************************/
422
512
 
423
 
size_t ReadPool::loadDataFromEXPs_rgid(const string & fofn, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))
 
513
size_t ReadPool::loadDataFromFOFNEXP_rgid(const string & fofn, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))
424
514
{
425
 
  FUNCSTART("size_t ReadPool::loadEXPs(const string & fofn, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))");
 
515
  FUNCSTART("size_t ReadPool::loadDataFromFOFNEXP(const string & fofn, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))");
426
516
 
427
517
  string justfilename, justpath;
428
518
  splitFullPathAndFileName(fofn,justpath,justfilename);
536
626
 *
537
627
 *
538
628
 *
 
629
 *
 
630
 *************************************************************************/
 
631
 
 
632
size_t ReadPool::loadDataFromEXP_rgid(const string & filename, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))
 
633
{
 
634
  FUNCSTART("size_t ReadPool::loadEXP(const string & filename, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))");
 
635
 
 
636
  uint32 no_files_ok=0;
 
637
  try{
 
638
    Read & actread = getRead(provideEmptyRead());
 
639
    actread.setReadGroupID(rgid);
 
640
 
 
641
    actread.loadDataFromEXP(filename,"");
 
642
    actread.transferSVTagsToClip(20,60);
 
643
 
 
644
    no_files_ok++;
 
645
 
 
646
    if(callback!=nullptr) {
 
647
      (*callback)(*this);
 
648
    }
 
649
  }
 
650
  catch(Flow){
 
651
  }
 
652
  catch(Notify n){
 
653
    n.handleError(THISFUNC);
 
654
  }
 
655
 
 
656
  FUNCEND();
 
657
 
 
658
  return no_files_ok;
 
659
}
 
660
 
 
661
 
 
662
/*************************************************************************
 
663
 *
 
664
 *
 
665
 *
539
666
 *************************************************************************/
540
667
size_t ReadPool::loadDataFromFASTQ_rgid(const string & filename, const ReadGroupLib::ReadGroupID rgid, bool countonly, void (*callback)(ReadPool &))
541
668
{
1869
1996
/*************************************************************************
1870
1997
 *
1871
1998
 * like above, but using save()
1872
 
 * This is for convert_project
 
1999
 * This is for miraconvert
1873
2000
 *
1874
2001
 *************************************************************************/
1875
2002
void ReadPool::saveAsMAF(ostream & ostr, bool alsoinvalids) const