← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/mira/contig.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/contig.C
265
265
  CON_us_steps_iric.resize(USCLOIRIC_END,0);
266
266
  CON_us_steps_drfc.clear();
267
267
  CON_us_steps_drfc.resize(USCLODRFC_END,0);
 
268
  CON_us_steps_cons.clear();
 
269
  CON_us_steps_cons.resize(USCLOCONS_END,0);
268
270
  CON_track_numins=0;
269
271
  CON_track_numdels=0;
270
272
 
618
620
    CON_stats.max_covperst[i]=0;
619
621
    CON_stats.avg_covperst[i]=0.0;
620
622
    CON_stats.readsperst[i]=0;
621
 
    CON_stats.readlenperst[i]=0;
 
623
    CON_stats.totalbasesperst[i]=0;
622
624
  }
623
625
}
624
626
 
646
648
 
647
649
  CON_stats.conlength=CON_counts.size();
648
650
 
649
 
  // avg. coverage of solexa and solid must be made differently than
650
 
  //  for sanger & 454 (because of merged reads)
651
 
  uint64 actcovsxa=0;
652
 
  uint64 sumcovsxa=0;
653
 
  uint64 actcovsid=0;
654
 
  uint64 sumcovsid=0;
655
 
 
656
651
  // using a pointer to CON_counts is (a lot) faster than accessing it
657
652
  //  via [] each time
658
653
  auto ccI=CON_counts.begin();
712
707
 
713
708
    CON_stats.max_coverage=max(CON_stats.max_coverage,
714
709
                               static_cast<uint32>(ccI->total_cov));
715
 
    for(uint32 i=0; i<ReadGroupLib::SEQTYPE_END; i++){
 
710
    for(uint32 i=0; i<ReadGroupLib::SEQTYPE_END; ++i){
 
711
      CON_stats.totalbasesperst[i]+=static_cast<uint64>(ccI->seqtype_cov[i]);
716
712
      CON_stats.max_covperst[i]=max(CON_stats.max_covperst[i],
717
713
                                    static_cast<uint32>(ccI->seqtype_cov[i]));
718
714
    }
719
715
 
720
 
    // 100% mapped Solexa/SOLiD reads must be accounted for a bit differently
721
 
 
722
 
    actcovsxa=static_cast<uint64>(ccI->seqtype_cov[ReadGroupLib::SEQTYPE_SOLEXA]);
723
 
    sumcovsxa+=actcovsxa;
724
 
 
725
 
    CON_stats.max_covperst[ReadGroupLib::SEQTYPE_SOLEXA]=
726
 
      max(CON_stats.max_covperst[ReadGroupLib::SEQTYPE_SOLEXA],static_cast<uint32>(actcovsxa));
727
 
 
728
 
    actcovsid=static_cast<uint64>(ccI->seqtype_cov[ReadGroupLib::SEQTYPE_ABISOLID]);
729
 
    sumcovsid+=actcovsid;
730
 
 
731
 
    CON_stats.max_covperst[ReadGroupLib::SEQTYPE_ABISOLID]=
732
 
      max(CON_stats.max_covperst[ReadGroupLib::SEQTYPE_ABISOLID],static_cast<uint32>(actcovsid));
733
 
 
734
716
    CON_stats.starInR+=ccI->star;
735
717
    CON_stats.NinR+=ccI->N;
736
718
  }
737
719
 
738
720
  // ok, in the section above, rail reads would have been also counted in
739
721
  //  (in .coverage)
740
 
  // compute reads numbers
 
722
  // compute reads statistics and remove rail bases from base statistics
741
723
 
742
724
  {
743
725
    auto cI=CON_reads.begin();
744
726
    for(;cI!=CON_reads.end();++cI){
745
 
      if(!cI->isRail()){
 
727
      if(cI->isRail()){
 
728
        CON_stats.totalbasesperst[cI->getSequencingType()]-=cI->getLenClippedSeq();
 
729
      }else{
746
730
 
747
731
        // Nono, backbones must be counted in!
748
732
        // Scenario: assembly at 6x used as backbone, assembled further 1x
751
735
        //  would want
752
736
        // && !cI->read.isBackbone()) {
753
737
 
754
 
        CON_stats.total_reads++;
 
738
        ++CON_stats.total_reads;
755
739
        CON_stats.readsperst[cI->getSequencingType()]++;
756
 
        if(cI.getORPID()>=0){
757
 
          CON_stats.readlenperst[cI->getSequencingType()]+=CON_readpool->getRead(cI.getORPID()).getLenClippedSeq();
758
 
          if(cI->hasQuality()){
759
 
            CON_stats.numreads_withqual++;
760
 
          }else{
761
 
            CON_stats.numreads_noqual++;
762
 
          }
 
740
        if(cI->hasQuality()){
 
741
          CON_stats.numreads_withqual++;
763
742
        }else{
764
 
          CON_stats.readlenperst[cI->getSequencingType()]+=cI->getLenClippedSeq();
 
743
          CON_stats.numreads_noqual++;
765
744
        }
766
745
      }
767
746
    }
768
747
  }
769
748
 
770
 
  CON_stats.readlenperst[ReadGroupLib::SEQTYPE_SOLEXA]=static_cast<uint32>(sumcovsxa);
771
 
  CON_stats.readlenperst[ReadGroupLib::SEQTYPE_ABISOLID]=static_cast<uint32>(sumcovsid);
772
 
 
 
749
  CON_stats.avg_coverage=0.0;
773
750
  if(!CON_counts.empty()){
774
 
    CON_stats.avg_coverage=0.0;
775
 
    for(uint32 i=0; i<ReadGroupLib::SEQTYPE_END; i++){
776
 
      if(CON_stats.readlenperst[i]>0.0){
777
 
        // subtracting the number of uncovered bases in the following firmula is a bad, bad idea and
778
 
        //  gives totally wrong numbers.
779
 
        //CON_stats.avg_covperst[i]=static_cast<double>(CON_stats.readlenperst[i])/(CON_counts.size()- CON_stats.numnocoverage);
780
 
        CON_stats.avg_covperst[i]=static_cast<double>(CON_stats.readlenperst[i])/(CON_counts.size());
781
 
        CON_stats.avg_coverage+=CON_stats.avg_covperst[i];
782
 
      }
 
751
    uint64 totalbases=0;
 
752
    for(uint32 sti=0; sti<ReadGroupLib::SEQTYPE_END; ++sti){
 
753
      CON_stats.avg_covperst[sti]=static_cast<double>(CON_stats.totalbasesperst[sti])/CON_counts.size();
 
754
      totalbases+=CON_stats.totalbasesperst[sti];
783
755
    }
 
756
    CON_stats.avg_coverage=static_cast<double>(totalbases)/CON_counts.size();
784
757
  }
785
758
 
786
759
  {
787
 
    vector<base_quality_t>::const_iterator I=consqual.begin();
788
760
    double qual=0.0;
789
 
    while(I!=consqual.end()){
790
 
      qual+=static_cast<double>(*I);
791
 
      I++;
792
 
    }
 
761
    for(auto & cqe : consqual) qual+=static_cast<double>(cqe);
793
762
    CON_stats.avg_conqual=static_cast<uint16>((qual/consseq.size())+.5);
794
763
  }
795
764
 
800
769
  // IUPACs are perhaps not set yet
801
770
  updateStatsFromConsensusTags(true, true, false, true, true);
802
771
 
 
772
  CON_stats.contains_long_repeats_only=CON_contains_long_repeats_only;
 
773
  CON_stats.islargecontig=0; // undefined, to be set by assembly_info
 
774
 
803
775
  CON_stats.statsvalid=true;
804
776
 
805
777
  FUNCEND();
1132
1104
/*************************************************************************
1133
1105
 *
1134
1106
 * addRead_wrapped might change some alignment parameters
 
1107
 *
 
1108
 * returns:
 
1109
 *   errstat for communicating errors
 
1110
 *   templateguess for communicating measure template info
 
1111
 *
1135
1112
 * as it returns from a number of different points within the loop,
1136
1113
 *  setting back original this function takes care to set everything back
1137
1114
 *  to normal if needed.
1140
1117
 
1141
1118
//#define BUGHUNT
1142
1119
 
1143
 
void Contig::addRead(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, errorstatus_t & errstat)
 
1120
void Contig::addRead(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, templateguessinfo_t & templateguess, errorstatus_t & errstat)
1144
1121
{
1145
 
  FUNCSTART("void Contig::addRead(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, errorstatus_t & errstat)");
 
1122
  FUNCSTART("void Contig::addRead(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, templateguess_t & templateguess, errorstatus_t & errstat)");
1146
1123
 
1147
1124
  //setCEBUGFlag(newid,refid);
1148
1125
 
1149
1126
  paranoiaBUGIF(aligncache.size()!=ReadGroupLib::SEQTYPE_END,
1150
1127
                MIRANOTIFY(Notify::INTERNAL, "aligncache is not size of available sequencing types???"));
1151
1128
 
 
1129
  templateguess.rgid.resetLibId();
 
1130
 
1152
1131
  // this is ugly: save the align params because
1153
1132
  //  addRead_wrapped() might change some *sigh*
1154
1133
  align_parameters oldalignparams=
1179
1158
                    direction_frnid,
1180
1159
                    newid_ismulticopy,
1181
1160
                    forcegrow,
 
1161
                    templateguess,
1182
1162
                    errstat);
1183
1163
 
1184
1164
    if(errstat.code!=ENOERROR){
 
1165
      // remove an eventual guess for template placement
 
1166
      templateguess.rgid.resetLibId();
 
1167
      templateguess.splace_seen=ReadGroupLib::SPLACE_UNKNOWN;
 
1168
      templateguess.tsize_seen=0;
 
1169
 
1185
1170
      // check if refid has been put into reads_affected by the contig
1186
1171
      // if not, warn and dump
1187
1172
      bool isin=false;
1244
1229
 
1245
1230
#define CEBUG(bla)
1246
1231
 
1247
 
void Contig::addRead_wrapped(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, errorstatus_t & errstat)
 
1232
void Contig::addRead_wrapped(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, templateguessinfo_t & templateguess, errorstatus_t & errstat)
1248
1233
{
1249
1234
  FUNCSTART("void Contig::addRead_wrapped(vector<Align> & aligncache, const AlignedDualSeqFacts * initialadsf, int32 refid, int32 newid, int32 direction_frnid, bool newid_ismulticopy, int32 forcegrow, errorstatus_t & errstat)");
1250
1235
 
1839
1824
 
1840
1825
  bool havematchingtemplatepartner=false;
1841
1826
 
1842
 
  if(CON_readpool->getRead(newid).hasTemplateInfo() > 0){
 
1827
  if(CON_readpool->getRead(newid).getTemplatePartnerID() != -1){
1843
1828
    // 1st: direction
1844
1829
#ifndef PUBLICQUIET
1845
1830
    cout << " tmplhand1 ";
1846
1831
    cout << CON_readpool->getRead(newid).getTemplatePartnerID() << " ";
1847
1832
    cout.flush();
1848
1833
#endif
1849
 
    auto pcrI=CON_reads.getIteratorOfReadpoolID(CON_readpool->getRead(newid).getTemplatePartnerID());
 
1834
    auto tppcrI=CON_reads.getIteratorOfReadpoolID(CON_readpool->getRead(newid).getTemplatePartnerID());
1850
1835
   // is partner already in contig?
1851
 
    if(pcrI==CON_reads.end()){
 
1836
    if(tppcrI==CON_reads.end()){
1852
1837
#ifndef PUBLICQUIET
1853
 
      cout << "nic ";
 
1838
      cout << "nic("
 
1839
           << CON_readpool->getRead(newid).getTemplatePartnerID();
 
1840
      if(CON_readpool->getRead(newid).getTemplatePartnerID()>=0){
 
1841
        cout << " - " << CON_readpool->getRead(CON_readpool->getRead(newid).getTemplatePartnerID()).getName();
 
1842
      }
 
1843
      cout << ") ";
1854
1844
      cout.flush();
1855
1845
#endif
1856
1846
    }else{
1871
1861
 
1872
1862
#ifndef PUBLICQUIET
1873
1863
        cout << "\ndni: " << static_cast<int32>(direction_newid_incontig);
1874
 
        cout << "\npdi: " << static_cast<int32>(pcrI.getReadDirection());
 
1864
        cout << "\npdi: " << static_cast<int32>(tppcrI.getReadDirection());
1875
1865
        cout << "\ntbd: " << static_cast<int32>(CON_readpool->getRead(newid).getTemplateBuildDirection());
1876
1866
#endif
1877
 
        if(direction_newid_incontig*pcrI.getReadDirection() != CON_readpool->getRead(newid).getTemplateBuildDirection()){
 
1867
        if(direction_newid_incontig*tppcrI.getReadDirection() != CON_readpool->getRead(newid).getTemplateBuildDirection()){
1878
1868
          // not direction wanted, not good
1879
1869
#ifndef PUBLICQUIET
1880
1870
          cout << "templ in wrong dir";
1899
1889
        }
1900
1890
      }
1901
1891
 
 
1892
      // for segment placement and template size check, we need the positions in the contig
 
1893
 
 
1894
      bool newreadisleft=true;
 
1895
 
 
1896
      // Note: play it safe while adding right extend ... take only
 
1897
      //  part of it
 
1898
      int32 leftrx=xcut;
 
1899
      int32 leftry=ycut;
 
1900
      if(direction_newid_incontig>0) {
 
1901
        leftry+=CON_readpool->getRead(newid).getRightExtend()*2/3;
 
1902
      }else{
 
1903
        leftrx-=CON_readpool->getRead(newid).getRightExtend()*2/3;
 
1904
      }
 
1905
 
 
1906
      int32 rightrx=tppcrI.getReadStartOffset();
 
1907
      int32 rightry=tppcrI.getReadStartOffset()+tppcrI->getLenClippedSeq();
 
1908
      if(direction_refid>0) {
 
1909
        rightry+=tppcrI->getRightExtend()*2/3;
 
1910
      }else{
 
1911
        rightrx-=tppcrI->getRightExtend()*2/3;
 
1912
      }
 
1913
 
 
1914
      if(rightrx<leftrx) {
 
1915
        swap(rightrx,leftrx);
 
1916
        swap(rightry,leftry);
 
1917
        newreadisleft=false;  // well, yeah, it is on the right
 
1918
      }
 
1919
 
 
1920
      int32 actinsertsize=rightry - leftrx;
 
1921
      if(rightry<leftry) actinsertsize=leftry - leftrx;
 
1922
 
 
1923
      // we now have all info needed for storing measured template info
 
1924
      bool guesstemplatesegplace=true;
 
1925
      templateguess.tsize_seen=0; // just as default
 
1926
      templateguess.splace_seen=ReadGroupLib::SPLACE_UNKNOWN; // just as default
 
1927
 
 
1928
      // Due to clipping at ends of reads,
 
1929
      // E.g.:
 
1930
      //         ----------->
 
1931
      //             <----------
 
1932
      // which can become
 
1933
      //                ---->
 
1934
      //             <----
 
1935
      //
 
1936
      // one should only guess on disjunct (non-overlapping) reads
 
1937
      // or if overlapping, the right read must start >= 10 bp
 
1938
      //  later than the left read
 
1939
      if(leftry<rightrx){
 
1940
        // if overlapping
 
1941
        if(leftrx+10>=rightrx){
 
1942
          // and start of reads <= 10bp
 
1943
          guesstemplatesegplace=false;
 
1944
        }
 
1945
      }
 
1946
 
 
1947
      if(guesstemplatesegplace){
 
1948
        templateguess.rgid=CON_readpool->getRead(newid).getReadGroupID();  // this validates the entry (else its  rgid: 0)
 
1949
        templateguess.tsize_seen=actinsertsize;
 
1950
        // placement code is hardest
 
1951
        if(direction_newid_incontig*tppcrI.getReadDirection() < 0){
 
1952
          // FR or RF
 
1953
          if(newreadisleft){
 
1954
            if(direction_newid_incontig>0){
 
1955
              templateguess.splace_seen=ReadGroupLib::SPLACE_FR;
 
1956
            }else{
 
1957
              templateguess.splace_seen=ReadGroupLib::SPLACE_RF;
 
1958
            }
 
1959
          }else{
 
1960
            if(direction_newid_incontig>0){
 
1961
              templateguess.splace_seen=ReadGroupLib::SPLACE_RF;
 
1962
            }else{
 
1963
              templateguess.splace_seen=ReadGroupLib::SPLACE_FR;
 
1964
            }
 
1965
          }
 
1966
        }else{
 
1967
          // same direction ...
 
1968
          //
 
1969
          // I'm sure that some boolean logic could get through this, but I'm not inclined to think about it now
 
1970
          templateguess.splace_seen=ReadGroupLib::SPLACE_SU; // just as default
 
1971
          if((direction_newid_incontig>0 && newreadisleft==true)
 
1972
             || (direction_newid_incontig<0 && newreadisleft==false)){
 
1973
            if(CON_readpool->getRead(newid).getTemplateSegment()==1){
 
1974
              templateguess.splace_seen=ReadGroupLib::SPLACE_SF;
 
1975
            }else{
 
1976
              templateguess.splace_seen=ReadGroupLib::SPLACE_SB;
 
1977
            }
 
1978
          }else{
 
1979
            if(CON_readpool->getRead(newid).getTemplateSegment()==1){
 
1980
              templateguess.splace_seen=ReadGroupLib::SPLACE_SB;
 
1981
            }else{
 
1982
              templateguess.splace_seen=ReadGroupLib::SPLACE_SF;
 
1983
            }
 
1984
          }
 
1985
        }
 
1986
      }
 
1987
#ifndef PUBLICQUIET
 
1988
      cout << " nril " << newreadisleft
 
1989
           << " " << leftrx << "-" << leftry
 
1990
           << " " << rightrx << "-" << rightry
 
1991
           << " " << templateguess;
 
1992
      cout.flush();
 
1993
#endif
 
1994
 
 
1995
 
1902
1996
      // ** Check for segment placement **
1903
 
      if(CON_readpool->getRead(newid).getReadGroupID().getSegmentPlacementCode() == ReadGroupLib::SPLACE_SF
1904
 
         || CON_readpool->getRead(newid).getReadGroupID().getSegmentPlacementCode() == ReadGroupLib::SPLACE_SB){
1905
 
#ifndef PUBLICQUIET
1906
 
        cout << " spl";
1907
 
        cout.flush();
1908
 
#endif
 
1997
      // unknown does obviously need no check, "samedir unknown" already checked implicitly by direction check above
 
1998
      {
 
1999
        auto spc=CON_readpool->getRead(newid).getReadGroupID().getSegmentPlacementCode();
 
2000
        if(spc != ReadGroupLib::SPLACE_UNKNOWN
 
2001
           && spc != ReadGroupLib::SPLACE_SU
 
2002
           && templateguess.splace_seen != ReadGroupLib::SPLACE_UNKNOWN
 
2003
           && templateguess.splace_seen != ReadGroupLib::SPLACE_SU){
 
2004
#ifndef PUBLICQUIET
 
2005
          cout << " spl";
 
2006
          cout.flush();
 
2007
#endif
 
2008
          // if reads overlap, stay cool and don't check as it could be that sequencing errors and clipping
 
2009
          //  muddy the picture
 
2010
          // That is: the +10 bp criterion of above is not taken into account here, only the non-overlapping part
 
2011
 
 
2012
          bool hasplacementerror=false;
 
2013
          if(leftry<rightrx){
 
2014
#ifndef PUBLICQUIET
 
2015
            cout << " chk";
 
2016
            cout.flush();
 
2017
#endif
 
2018
            if(spc!=templateguess.splace_seen) hasplacementerror=true;
 
2019
          }
 
2020
 
 
2021
          if(hasplacementerror){
 
2022
            errstat.code=ESEGMENTPLACEMENT;
 
2023
 
 
2024
            if(xcut<0) xcut=0;
 
2025
            if(static_cast<uint32>(ycut)>CON_counts.size()){
 
2026
              ycut=CON_counts.size();
 
2027
            }
 
2028
            gettimeofday(&us_start,nullptr);
 
2029
            if(CON_readpool->getRead(refid).isRail())  {
 
2030
              getRailsAsReadsAffected(refid, errstat.reads_affected, xcut, ycut);
 
2031
            }else{
 
2032
              getReadORPIDsAtContigPosition(errstat.reads_affected, xcut, ycut);
 
2033
            }
 
2034
            CON_us_steps[USCLO_GRACP]+=diffsuseconds(us_start);
 
2035
            FUNCEND();
 
2036
            return;
 
2037
          }
 
2038
        }
1909
2039
      }
1910
2040
 
1911
2041
      // ** Check for template size **
1916
2046
        cout.flush();
1917
2047
#endif
1918
2048
 
1919
 
        // Note: play it safe while adding right extend ... take only
1920
 
        //  part of it
1921
 
        int32 leftrx=xcut;
1922
 
        int32 leftry=ycut;
1923
 
        if(direction_newid_incontig>0) {
1924
 
          leftry+=CON_readpool->getRead(newid).getRightExtend()*3/2;
1925
 
        }else{
1926
 
          leftrx-=CON_readpool->getRead(newid).getRightExtend()*3/2;
1927
 
        }
1928
 
 
1929
 
        int32 rightrx=pcrI.getReadStartOffset();
1930
 
        int32 rightry=pcrI.getReadStartOffset()+pcrI->getLenClippedSeq();
1931
 
        if(direction_refid>0) {
1932
 
          rightry+=pcrI->getRightExtend()*3/2;
1933
 
        }else{
1934
 
          rightrx-=pcrI->getRightExtend()*3/2;
1935
 
        }
1936
 
 
1937
 
        if(rightrx<leftrx) {
1938
 
          swap(rightrx,leftrx);
1939
 
          swap(rightry,leftry);
1940
 
        }
1941
 
 
1942
 
        int32 actinsertsize=rightry - leftrx;
1943
 
        if(rightry<leftry) actinsertsize=leftry - leftrx;
1944
 
 
1945
2049
        // allow 15% error in calculated insert size
1946
2050
        int32 aisp10=actinsertsize+actinsertsize*15/100;
1947
2051
        int32 aism10=actinsertsize-actinsertsize*15/100;
1948
2052
 
1949
 
        int32 tif=pcrI->getInsizeFrom();
1950
 
        int32 tit=pcrI->getInsizeTo();
 
2053
        int32 tif=tppcrI->getInsizeFrom();
 
2054
        int32 tit=tppcrI->getInsizeTo();
1951
2055
 
1952
 
        if(pcrI->getInsizeFrom() >=0 && aisp10 < tif){
 
2056
        if(tppcrI->getInsizeFrom() >=0 && aisp10 < tif){
1953
2057
        // distance too small
1954
2058
#ifndef PUBLICQUIET
1955
 
          cout << "templ too small: " << tif << " min allowed, got " << aisp10;
 
2059
          cout << "templ too small: " << tif << " min allowed, got " << aism10 << "-" << aisp10;
1956
2060
#endif
1957
2061
          havematchingtemplatepartner=false;
1958
2062
          if(!CON_readpool->getRead(newid).getReadGroupID().getTSInfoOnly()){
1959
 
            errstat.code=ETEMPLATESIZE;
 
2063
            errstat.code=ETEMPLATESIZELT;
1960
2064
          }
1961
2065
        }
1962
 
        if(errstat.code != ETEMPLATESIZE && pcrI->getInsizeTo() >=0 && aism10 > tit){
 
2066
        if(errstat.code == ENOERROR && tppcrI->getInsizeTo() >=0 && aism10 > tit){
1963
2067
          // distance too big
1964
2068
#ifndef PUBLICQUIET
1965
 
          cout << "templ too big: " << tit << " max allowed, got " << aism10;
 
2069
          cout << "templ too big: " << tit << " max allowed, got " << aism10 << "-" << aisp10;
1966
2070
#endif
1967
2071
          havematchingtemplatepartner=false;
1968
2072
          if(!CON_readpool->getRead(newid).getReadGroupID().getTSInfoOnly()){
1969
 
            errstat.code=ETEMPLATESIZE;
 
2073
            errstat.code=ETEMPLATESIZEGT;
1970
2074
          }
1971
2075
        }
1972
 
        if(errstat.code == ETEMPLATESIZE){
 
2076
        if(errstat.code != ENOERROR){
1973
2077
          if(xcut<0) xcut=0;
1974
2078
          if(static_cast<uint32>(ycut)>CON_counts.size()){
1975
2079
            ycut=CON_counts.size();
2771
2875
//      canmap=false;
2772
2876
//      }
2773
2877
 
 
2878
      //cout << "\ncanmap2: " << canmap << endl;
 
2879
 
2774
2880
      // maybe we'd like to have reads at contig ends not mapped (e.g.
2775
2881
      //  for scaffolding)
2776
2882
      if(canmap
2788
2894
          dist=CON_readpool->getRead(newid).getInsizeTo();
2789
2895
        }
2790
2896
 
2791
 
//      cout << "mpl: " << CON_markerpositions[CON_mpindex_msrkceu_left]
2792
 
//         << "\tmpr: " << CON_markerpositions[CON_mpindex_msrkceu_right]
2793
 
//         << "\nxcut: " << xcut
2794
 
//         << "\tycut: " << ycut
2795
 
//         << "\tdist: " << dist
2796
 
//         << "\nl: " << dist+CON_markerpositions[CON_mpindex_msrkceu_left]
2797
 
//         << "\tl: " << CON_markerpositions[CON_mpindex_msrkceu_right]-dist
2798
 
//      ;
 
2897
        //cout << "mpl: " << CON_markerpositions[CON_mpindex_msrkceu_left]
 
2898
        //     << "\tmpr: " << CON_markerpositions[CON_mpindex_msrkceu_right]
 
2899
        //     << "\nxcut: " << xcut
 
2900
        //     << "\tycut: " << ycut
 
2901
        //     << "\tdist: " << dist
 
2902
        //     << "\nl: " << dist+CON_markerpositions[CON_mpindex_msrkceu_left]
 
2903
        //     << "\tl: " << CON_markerpositions[CON_mpindex_msrkceu_right]-dist
 
2904
        //  ;
2799
2905
 
2800
2906
        // if no distance or when read has no paired-end partner, we could map
2801
2907
        //  anyway and would not need to check further
2802
2908
        if(dist>0){
2803
2909
          // check whether we are in boundaries, if yes, do not map
2804
2910
          if(xcut <= dist+CON_markerpositions[CON_mpindex_msrkceu_left]
2805
 
             || ycut >= CON_markerpositions[CON_mpindex_msrkceu_left]-dist){
 
2911
             || ycut >= CON_markerpositions[CON_mpindex_msrkceu_right]-dist){
2806
2912
            canmap=false;
2807
2913
          }
2808
2914
        }
2809
2915
 
2810
 
        //cout << "\ncanmap: " << canmap << endl;
 
2916
        //cout << "\ncanmap2: " << canmap << endl;
2811
2917
      }
2812
2918
 
2813
2919
      // map if possible, else the read will be normally added
4515
4621
 
4516
4622
/*************************************************************************
4517
4623
 *
4518
 
 *
4519
 
 *
 
4624
 * "doshiftreads==false" is a hack for trimMapOverhang() which needs
 
4625
 *  chompFront() functionality WITHOUT shifting reads
4520
4626
 *
4521
4627
 *************************************************************************/
4522
4628
 
4523
4629
//#define CEBUG(bla)   {cout << bla; cout.flush();}
4524
 
int32 Contig::chompFront(int32 maxchecklen)
 
4630
int32 Contig::chompFront(int32 maxchecklen, bool doshiftreads)
4525
4631
{
4526
 
  FUNCSTART("int32 Contig::chompFront(int32 maxchecklen)");
 
4632
  FUNCSTART("int32 Contig::chompFront(int32 maxchecklen, bool doshiftreads)");
4527
4633
 
4528
4634
#ifdef CLOCK_STEPS_DRFC
4529
4635
  timeval us_tv;
4569
4675
    }
4570
4676
    if(CON_us_steps_drfc.size()) CON_us_steps_drfc[USCLODRFC_SDT]+=diffsuseconds(us_tv);
4571
4677
    // push down the offsets of the other reads
4572
 
    gettimeofday(&us_tv,nullptr);
4573
 
    CON_reads.shiftReads(1,-frontdeletions);
4574
 
    if(CON_us_steps_drfc.size()) CON_us_steps_drfc[USCLODRFC_SR]+=diffsuseconds(us_tv);
 
4678
    if(doshiftreads){
 
4679
      gettimeofday(&us_tv,nullptr);
 
4680
      CON_reads.shiftReads(1,-frontdeletions);
 
4681
      if(CON_us_steps_drfc.size()) CON_us_steps_drfc[USCLODRFC_SR]+=diffsuseconds(us_tv);
 
4682
    }
4575
4683
    // same thing for contig marker positions
4576
4684
    gettimeofday(&us_tv,nullptr);
4577
4685
    shiftMarkerPositions(1,-frontdeletions);
4653
4761
#define CEBUG(bla)
4654
4762
 
4655
4763
 
 
4764
/*************************************************************************
 
4765
 *
 
4766
 *
 
4767
 *
 
4768
 *
 
4769
 *************************************************************************/
 
4770
//#define CEBUG(bla)   {cout << bla; cout.flush();}
 
4771
void Contig::trimMapOverhang()
 
4772
{
 
4773
  FUNCSTART("void Contig::trimMapOverhang()");
 
4774
 
 
4775
  //cout << "MMMMMMMMMM\n";
 
4776
  //CON_reads.debugDump(false);
 
4777
 
 
4778
  if(getNumBackbones()==0) return;
 
4779
 
 
4780
  int32 firstknownbbpos=-1;
 
4781
  int32 lastknownbbpos=-1;
 
4782
  for(auto crI= CON_reads.begin(); crI!=CON_reads.end(); ++crI){
 
4783
    if(crI->isBackbone()){
 
4784
      if(static_cast<int32>(crI.getReadStartOffset()+crI->getLenClippedSeq()) > lastknownbbpos){
 
4785
        lastknownbbpos=crI.getReadStartOffset()+crI->getLenClippedSeq();
 
4786
      }
 
4787
      if(firstknownbbpos<0){
 
4788
        firstknownbbpos=crI.getReadStartOffset();
 
4789
      }
 
4790
    }
 
4791
  }
 
4792
 
 
4793
  CEBUG("fkbb: " << firstknownbbpos << endl);
 
4794
  CEBUG("lkbb: " << lastknownbbpos << endl);
 
4795
  CEBUG("gcl: " << getContigLength() << endl);
 
4796
 
 
4797
  BUGIFTHROW(firstknownbbpos<0,"firstknownbbpos " << firstknownbbpos << " <0 ???");
 
4798
  BUGIFTHROW(lastknownbbpos<0,"lastknownbbpos " << lastknownbbpos << " <0 ???");
 
4799
 
 
4800
  // check: are there reads completely outside the backbone? if yes: no-can-do
 
4801
  for(auto crI= CON_reads.begin(); crI!=CON_reads.end(); ++crI){
 
4802
    if(!crI->isBackbone()){
 
4803
      if(crI.getReadStartOffset() < firstknownbbpos
 
4804
         && crI.getReadStartOffset()+crI->getLenClippedSeq() <= firstknownbbpos){
 
4805
        // oooops, impossible at front
 
4806
        return;
 
4807
      }
 
4808
      if(crI.getReadStartOffset() >= lastknownbbpos){
 
4809
        // oooops, impossible at end
 
4810
        return;
 
4811
      }
 
4812
    }
 
4813
  }
 
4814
 
 
4815
  bool docorrect=false;
 
4816
  // we're fine, now get things rolling
 
4817
  //
 
4818
  // first, the back of the contig if needed
 
4819
  if(lastknownbbpos < getContigLength()){
 
4820
    docorrect=true;
 
4821
    // zero out CON_counts
 
4822
    auto ccI=CON_counts.end(); // cannot use rbegin, not implemented in HDeque
 
4823
    CEBUG("numdelsteps: " << getContigLength()-lastknownbbpos << endl);
 
4824
    for(auto numdelsteps=getContigLength()-lastknownbbpos; numdelsteps!=0; --numdelsteps){
 
4825
      --ccI;
 
4826
      *ccI=CON_concounts_zero;
 
4827
    }
 
4828
 
 
4829
    // shorten overhanging reads
 
4830
    for(auto crI= CON_reads.begin(); crI!=CON_reads.end(); ++crI){
 
4831
      int32 shorten=static_cast<int32>(crI.getReadStartOffset()+crI->getLenClippedSeq())-lastknownbbpos;
 
4832
      CEBUG(crI->getName() << " shortenb " << shorten << endl);
 
4833
      if(shorten>0){
 
4834
        // very hacky, but we're forcing our hand on the reads of the PCR container here
 
4835
        Read & ncr=const_cast<Read &>(*crI);
 
4836
        if(crI.getReadDirection()>0){
 
4837
          ncr.setRQClipoff(crI->getRightClipoff()-shorten);
 
4838
        }else{
 
4839
          ncr.setLQClipoff(crI->getLeftClipoff()+shorten);
 
4840
        }
 
4841
      }
 
4842
    }
 
4843
  }
 
4844
 
 
4845
  // now the front of the contig if needed
 
4846
  if(firstknownbbpos>0){
 
4847
    docorrect=true;
 
4848
    // zero out CON_counts
 
4849
    auto ccI=CON_counts.begin();
 
4850
    for(auto numdelsteps=firstknownbbpos; numdelsteps!=0; --numdelsteps, ++ccI){
 
4851
      *ccI=CON_concounts_zero;
 
4852
    }
 
4853
 
 
4854
    // shorten overhanging reads
 
4855
    for(auto crI= CON_reads.begin(); crI!=CON_reads.end(); ++crI){
 
4856
      int32 shorten=firstknownbbpos-static_cast<int32>(crI.getReadStartOffset());
 
4857
      CEBUG(crI->getName() << "\tfkbbp: " << firstknownbbpos << "\trso: " << crI.getReadStartOffset() <<"\tshortenf " << shorten << endl);
 
4858
      if(shorten>0){
 
4859
        // very hacky, but we're forcing our hand on the reads of the PCR container here
 
4860
        Read & ncr=const_cast<Read &>(*crI);
 
4861
        if(crI.getReadDirection()>0){
 
4862
          ncr.setLQClipoff(crI->getLeftClipoff()+shorten);
 
4863
        }else{
 
4864
          ncr.setRQClipoff(crI->getRightClipoff()-shorten);
 
4865
        }
 
4866
      }
 
4867
    }
 
4868
  }
 
4869
 
 
4870
  if(docorrect){
 
4871
    chompBack(-1);
 
4872
 
 
4873
    if(firstknownbbpos>0){
 
4874
      // hacky thing
 
4875
      chompFront(-1,false);
 
4876
      // hacky thing
 
4877
      CON_reads.shiftReadsBounceZero(0,-firstknownbbpos);
 
4878
    }
 
4879
  }
 
4880
}
 
4881
//#define CEBUG(bla)
 
4882
 
4656
4883
 
4657
4884
/*************************************************************************
4658
4885
 *
4774
5001
    CEBUG("New element:\t");
4775
5002
    CEBUG("Offset: " << rle.offset_in_contig << "\tID: " << rle.id);
4776
5003
    CEBUG("\tDir: " << rle.direction << endl);
4777
 
    CEBUG("Read:\n" << CON_readpool->getRead(rle.id).getName() << endl);
 
5004
    CEBUG("Read: " << CON_readpool->getRead(rle.id).getName() << endl);
4778
5005
 
4779
5006
    auto pcrI=CON_reads.placeRead(CON_readpool->getRead(rle.id),
4780
5007
                                  rle.id,
4782
5009
                                  rle.direction);
4783
5010
    computedclen=max(computedclen,static_cast<uint32>(pcrI.getReadStartOffset())+static_cast<uint32>(pcrI->getLenClippedSeq()));
4784
5011
 
4785
 
    auto coveragemultiplier=CON_readpool->getRead(pcrI->getStrainID()).getDigiNormMultiplier();
 
5012
    auto coveragemultiplier=CON_readpool->getRead(rle.id).getDigiNormMultiplier();
4786
5013
    CON_readsperstrain[pcrI->getStrainID()]+=coveragemultiplier;
4787
5014
    CON_readsperreadgroup[pcrI->getReadGroupID().getLibId()]+=coveragemultiplier;
4788
5015
    if(pcrI->getLenClippedSeq() > CON_longestreadseen) CON_longestreadseen=pcrI->getLenClippedSeq();
4795
5022
    }
4796
5023
  }
4797
5024
 
 
5025
  CEBUG("end of list" << endl);
 
5026
 
4798
5027
  // CAF as written by gap5 may start at another offset than 0
4799
5028
  // MIRA expects the first read to have an offset of 0, therefore correct for that
4800
5029
  int32 shiftoffset=-CON_reads.begin().getReadStartOffset();
4854
5083
    }
4855
5084
  }
4856
5085
 
 
5086
  CEBUG("CON_counts before trimming: " << CON_counts.size() << endl);
 
5087
 
4857
5088
  // trim end of contig if needed
4858
5089
  uint64 numpops=0;
4859
5090
  while(!CON_counts.empty() &&
4880
5111
    CON_fixedconsqual.erase(CON_fixedconsqual.begin(),CON_fixedconsqual.begin()+numpops);
4881
5112
  }
4882
5113
 
 
5114
  BUGIFTHROW(CON_counts.empty(),"After trimming, contig has zero length?");
4883
5115
 
4884
5116
  // must copy each tag as we have to "upgrade" a tag_t to a consensustag_t
4885
5117
  {
6886
7118
    }
6887
7119
    break;
6888
7120
  }
 
7121
  case Contig::ENOTCALLED : {
 
7122
    if(longmsg){
 
7123
      cout << "\t-\tnot called\n";
 
7124
    }else{
 
7125
      cout << ' ';
 
7126
    }
 
7127
    break;
 
7128
  }
6889
7129
  case Contig::ENOALIGN : {
6890
7130
    if(longmsg){
6891
7131
      cout << "\t-\tno align found\n";
6906
7146
    if(longmsg){
6907
7147
      cout << "\t-\ttemplate in wrong direction\n";
6908
7148
    }else{
6909
 
      cout << 't';
6910
 
    }
6911
 
    break;
6912
 
  }
6913
 
  case Contig::ETEMPLATESIZE : {
6914
 
    if(longmsg){
6915
 
      cout << "\t-\tmismatch in template size\n";
6916
 
    }else{
6917
 
      cout << 's';
 
7149
      cout << 'T';
 
7150
    }
 
7151
    break;
 
7152
  }
 
7153
  case Contig::ESEGMENTPLACEMENT : {
 
7154
    if(longmsg){
 
7155
      cout << "\t-\tplacement of segment wrong\n";
 
7156
    }else{
 
7157
      cout << 'P';
 
7158
    }
 
7159
    break;
 
7160
  }
 
7161
  case Contig::ETEMPLATESIZELT : {
 
7162
    if(longmsg){
 
7163
      cout << "\t-\tmismatch in template size (<)\n";
 
7164
    }else{
 
7165
      cout << '<';
 
7166
    }
 
7167
    break;
 
7168
  }
 
7169
  case Contig::ETEMPLATESIZEGT : {
 
7170
    if(longmsg){
 
7171
      cout << "\t-\tmismatch in template size (>)\n";
 
7172
    }else{
 
7173
      cout << '>';
6918
7174
    }
6919
7175
    break;
6920
7176
  }
6922
7178
    if(longmsg){
6923
7179
      cout << "\t-\tmismatch in SRMB zone\n";
6924
7180
    }else{
6925
 
      cout << 'p';
 
7181
      cout << 'R';
6926
7182
    }
6927
7183
    break;
6928
7184
  }