home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / DMegaSequence.cpp < prev    next >
Encoding:
Text File  |  1996-07-05  |  28.1 KB  |  1,237 lines  |  [TEXT/R*ch]
  1. // DSequence.cp 
  2. // d.g.gilbert
  3.  
  4.  
  5.  
  6. #include "DMegaSequence.h"
  7. #include "DSeqFile.h"
  8. #include "DREnzyme.h"
  9. #include <ncbi.h>
  10. #include <dgg.h>
  11. #include <DUtil.h>
  12. #include <DFile.h>
  13. #include "ureadseq.h"
  14.  
  15.  
  16.  
  17.  
  18.  
  19.  
  20.  
  21.  
  22.  
  23.  
  24. // DMegaSequence -----------------
  25.  
  26.  
  27. DMegaSequence::DMegaSequence()
  28. {
  29.     fIsMega        = TRUE;
  30. }
  31.  
  32.  
  33. DMegaSequence::~DMegaSequence()  
  34. {
  35. }
  36.  
  37.  
  38. DObject* DMegaSequence::Clone() // override 
  39. {    
  40.     DMegaSequence* aSeq= (DMegaSequence*) DSequence::Clone();
  41. #if 0
  42.     aSeq->fBases= NULL; aSeq->SetBases( fBases, true);  
  43.     aSeq->fMasks= NULL; aSeq->SetMasks( fMasks, true); 
  44.     aSeq->fInfo=  StrDup( fInfo);
  45.     aSeq->fName=  StrDup( fName);
  46.     aSeq->fSearchRec.targ= StrDup( fSearchRec.targ);
  47. #endif
  48.     return aSeq;
  49. }
  50.  
  51. void DMegaSequence::CopyContents( DSequence* fromSeq)
  52. {
  53.     DSequence::CopyContents( fromSeq);
  54. }
  55.  
  56.  
  57. char* DMegaSequence::Bases() const {
  58.     return fBases; // fix later
  59.     }
  60.  
  61. char* DMegaSequence::Masks() const {
  62.     return fMasks; // fix later
  63.     }
  64.     
  65.     
  66. enum { kHugePartSize = 3000000 }; // = kMegaPartSize 
  67.  
  68. void DMegaSequence::countparts( char* seq, long& len, long& nparts, long& partsize)
  69. {
  70.     char* cp= seq;
  71.     partsize= kHugePartSize;  
  72.     len= 0;
  73.     nparts= 0;
  74.     if (!seq) return;
  75.     
  76.     while ( *cp != '\0' ) {
  77.         long  n;
  78.         for (n=0; *cp && n<kHugePartSize; n++) cp++;
  79.         len += n;
  80.         nparts++;    
  81.         }    
  82.     if (nparts == 1) partsize= len;
  83. }
  84.  
  85.  
  86. void DMegaSequence::SetBases( char*& theBases, Boolean duplicate)
  87. {
  88.     long  newlen= (theBases) ? StrLen(theBases) : 0;
  89.     if (fBases) { 
  90.         if (fLength>0 && 
  91.             (fLength != newlen || MemCmp( theBases, fBases, newlen)!= 0))
  92.                 fChanged= TRUE;
  93.         MemFree(fBases); 
  94.         fBases= NULL; fLength= 0; 
  95.         }
  96.  
  97.     if (theBases) {
  98.         countparts( theBases, fLength, fNparts, fPartsize);
  99.         
  100.         if (fNparts==1) {
  101.             if (duplicate) fBases = StrDup(theBases);  
  102.             else { fBases= theBases; theBases= NULL; }
  103.             }
  104.             
  105.         else if (fNparts>1) {
  106.                 // temp for testing !
  107.             fLength= fPartsize;
  108.             fBases= (char*)MemNew( fLength+1);
  109.             MemCpy( fBases, theBases, fLength);
  110.             fBases[fLength]= 0;
  111.             }
  112.         
  113.         }
  114.     fModTime= gUtil->time();
  115.     fMasksOk= false; // assume SetBases allways called before SetMask...
  116. }
  117.  
  118.  
  119. void DMegaSequence::UpdateBases( char* theBases, long theLength)
  120. {
  121.     fBases =  theBases;  
  122.     fLength= theLength;
  123.     fMasksOk= false; // assume this allways called before SetMask...
  124. }
  125.  
  126.  
  127. void DMegaSequence::SetMasks( char*& theMasks, Boolean duplicate)
  128. {
  129.     if (fMasks) MemFree(fMasks);
  130.     if (theMasks) {
  131.         if (duplicate) fMasks = StrDup(theMasks);  
  132.         else { fMasks= theMasks; theMasks= NULL; }
  133.         FixMasks();
  134.         }
  135.     else {
  136.         fMasks= NULL;
  137.         fMasksOk= false;
  138.         }
  139. }
  140.  
  141. void DMegaSequence::FixMasks()
  142. {
  143.     ulong i, mlen= 0;
  144.     if (!fMasks) { 
  145.         mlen= 1+fLength;
  146.         fMasks= (char*) MemGet(mlen,true);  
  147.         for (i=0; i<mlen; i++) fMasks[i] |= 0x80; // set high bit so Strxxx sees it as non-null
  148.         fMasks[mlen-1]= 0;
  149.         }
  150.     else { 
  151.         mlen= 1+StrLen(fMasks);
  152.         if (fLength+1 != mlen) {
  153.             fMasks= (char*) Nlm_MemExtend( fMasks, fLength+1, mlen);
  154.             mlen= fLength+1; 
  155.             for (i=0; i<mlen; i++) fMasks[i] |= 0x80; 
  156.             fMasks[mlen-1]= 0;
  157.             } 
  158.         }
  159.     fMasksOk= true;
  160. }
  161.  
  162. short DMegaSequence::MaskAt(long baseindex, short masklevel)
  163. {
  164.     short b = kMaskEmpty;
  165.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength) {
  166.         b= (unsigned char) fMasks[baseindex];
  167.         switch (masklevel) {
  168.             case 0: b &= 0x7f; break; // return full mask - top (0x80) bit
  169.             case 1: b &= 1;  break;  
  170.             case 2: b &= 2;  break;
  171.             case 3: b &= 4;  break;
  172.             case 4: b &= 8;  break;
  173.             default: b = kMaskEmpty; break;
  174.             }
  175.         }
  176.     return b;
  177. }
  178.  
  179. Boolean DMegaSequence::IsMasked(unsigned char maskbyte, short masklevel)
  180. {
  181.     short b= (unsigned char) maskbyte;
  182.     switch (masklevel) {
  183.         case 0: b &= 0x7f; break; // return full mask - top (0x80) bit
  184.         case 1: b &= 1;  break;  
  185.         case 2: b &= 2;  break;
  186.         case 3: b &= 4;  break;
  187.         case 4: b &= 8;  break;
  188.         default: b = kMaskEmpty; break;
  189.         }
  190.     return (b > 0);
  191. }
  192.  
  193.  
  194. inline void SetMaskByte( char& b, short masklevel, short maskval)
  195. {
  196.     switch (masklevel) {
  197.         case 1: if (maskval) b |= 1; else b &= ~1; break;
  198.         case 2: if (maskval) b |= 2; else b &= ~2; break;
  199.         case 3: if (maskval) b |= 4; else b &= ~4; break;
  200.         case 4: if (maskval) b |= 8; else b &= ~8; break;
  201.         // reserve bits 5, 6, 7 & 8 for now 
  202.         // -- 7 & 8 may go unused to store data in printchar form
  203.         default: break;
  204.         }
  205.     //if (b) masklevel--; // debugging b val
  206. }
  207.  
  208. void DMegaSequence::SetMaskAt(long baseindex, short masklevel, short maskval)
  209. {
  210.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  211.         ::SetMaskByte(fMasks[baseindex],masklevel, maskval);
  212. }
  213.  
  214.  
  215. void DMegaSequence::SetMask(short masklevel, short maskval)
  216. {
  217.     long  baseindex;
  218.     if (fMasks && fMasksOk) 
  219.      for (baseindex=0; baseindex<fLength; baseindex++) 
  220.         ::SetMaskByte(fMasks[baseindex], masklevel, maskval);
  221. }
  222.  
  223.  
  224. inline void FlipMaskByte( char& b, short masklevel)
  225. {
  226.     switch (masklevel) {
  227.         case 1: b ^= 1;  break;  
  228.         case 2: b ^= 2;  break;
  229.         case 3: b ^= 4;  break;
  230.         case 4: b ^= 8;  break;
  231.         default:  break;
  232.         }
  233.     //if (b) masklevel--; // debugging b val
  234. }
  235.  
  236. void DMegaSequence::FlipMaskAt(long baseindex, short masklevel)
  237. {
  238.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  239.         ::FlipMaskByte(fMasks[baseindex],masklevel);
  240. }
  241.  
  242. void DMegaSequence::FlipMask(short masklevel)
  243. {
  244.     long  baseindex;
  245.     if (fMasks && fMasksOk) 
  246.      for (baseindex=0; baseindex<fLength; baseindex++)  
  247.         FlipMaskByte(fMasks[baseindex],masklevel);
  248. }
  249.  
  250.  
  251.  
  252. void DMegaSequence::UpdateFlds()
  253. {
  254.     long         seqlen;
  255.     short     kind;
  256.     unsigned long check;
  257.     fLength= StrLen(fBases);
  258.     GetSeqStats(fBases, fLength, &kind, &check, &seqlen);
  259.     //fValidLength= seqlen; //?? is fLength == StrLen(fBases) or is it <= that, only valid bases??
  260.     fKind            = kind;
  261.     fChecksum    = check;
  262.     fChanged    = FALSE;
  263. }
  264.  
  265.  
  266.  
  267.  
  268.  
  269. void DMegaSequence::SetSelection( long start, long nbases)
  270. {
  271.     fSelStart= start;
  272.     fSelBases= nbases;
  273. }
  274.  
  275. void DMegaSequence::ClearSelection()
  276. {
  277.     fSelStart= 0;
  278.     fSelBases= 0;
  279. }
  280.  
  281. void    DMegaSequence::GetSelection( long& start, long& nbases)
  282. {
  283.     if (fSelBases==0) {
  284.         start= 0;
  285.         nbases= fLength;
  286.         }    
  287.     else if (fLength>0) {
  288.         start= fSelStart;
  289.         nbases= Min( fSelBases, fLength - fSelStart);
  290.         }    
  291.     else {
  292.         start= 0;
  293.         nbases= 0;
  294.         }
  295. }
  296.  
  297.  
  298. DSequence* DMegaSequence::CopyRange()
  299. {  
  300.             //note: start == 0 == 1st position 
  301.     DMegaSequence* aSeq= (DMegaSequence*) Clone();
  302.     if (fSelStart > 0 || fSelBases > 0) {
  303.         char* h= aSeq->fBases;
  304.         char* m= aSeq->fMasks;
  305.         Boolean domask= (fMasks && fMasksOk);
  306.  
  307.         long len= aSeq->fLength;
  308.         if (fSelStart > 0) {
  309.             len -= fSelStart;
  310.             MemMove( h, h + fSelStart, len);
  311.             if (domask) MemMove( m, m + fSelStart, len);
  312.             }
  313.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  314.         h= (char*) MemMore( h, len+1);
  315.         h[len]= '\0'; //((char*)h+len) = '\0';
  316.         aSeq->fBases= h;
  317.         if (domask) {
  318.             m= (char*) MemMore( m, len+1);
  319.             m[len]= '\0';    
  320.             aSeq->fMasks= m;
  321.             }
  322.         }
  323.     aSeq->UpdateFlds();
  324.     aSeq->Modified();
  325.     return aSeq;
  326. }
  327.  
  328.  
  329. DSequence* DMegaSequence::CutRange()
  330. {
  331. /* aSeq.Cutrange means:  
  332.         (a) selected range (fSelStart..fSelStart+fSelBases) is copied into Result; &&
  333.         (b) selected range is removed from this
  334. */
  335.  
  336.     DMegaSequence* aSeq= (DMegaSequence*) CopyRange();
  337.     if (fSelStart > 0 || fSelBases > 0) {
  338.         if (fSelStart<0) fSelStart= 0;
  339.         char* h= this->fBases;
  340.         char* m= this->fMasks;
  341.         Boolean domask= (fMasks && fMasksOk);
  342.         long len= this->fLength;
  343.         long top= len - fSelStart - fSelBases;
  344.         if (fSelBases > 0 && fSelBases < len) {
  345.             MemMove( h+fSelStart, h+fSelStart+fSelBases, top);
  346.             if (domask) MemMove( m+fSelStart, m+fSelStart+fSelBases, top);
  347.             }
  348.         len= fSelStart + top;
  349.         h= (char*) MemMore( h, len+1);
  350.         *((char*)h+len)= '\0';       //h[len]= '\0';
  351.         this->fBases= h;
  352.         if (domask) {
  353.             m= (char*) MemMore( m, len+1);
  354.             *((char*)m+len)= '\0';    
  355.             this->fMasks= m;
  356.             }
  357.         }
  358.     UpdateFlds();    
  359.     Modified();
  360.     return aSeq;
  361. }
  362.  
  363. void DMegaSequence::InsertSpacers( long insertPoint, long insertCount, char spacer)
  364. { 
  365.     long i;
  366.     char *cp;
  367.     char* h= this->fBases;
  368.     char* m= this->fMasks;
  369.     Boolean domask= (fMasks && fMasksOk);
  370.     long len= this->fLength;
  371.     if (insertPoint<0) insertPoint= 0;
  372.     else if (insertPoint>len-1) insertPoint= len-1;
  373.     h= (char*) MemMore( h, len+insertCount+1);
  374.     cp = h+insertPoint;
  375.     MemMove( cp+insertCount, cp, len-insertPoint);
  376.     for ( i= 0; i<insertCount; i++)  *cp++= spacer;
  377.     h[len+insertCount]= '\0';
  378.     this->fBases= h;
  379.     if (domask) {
  380.         m= (char*) MemMore( m, len+insertCount+1);
  381.         cp = m+insertPoint;
  382.         char cmask= *cp;
  383.         MemMove( cp+insertCount, cp, len-insertPoint);
  384.         for ( i= 0; i<insertCount; i++)  *cp++= cmask;
  385.         m[len+insertCount]= '\0';
  386.         this->fMasks= m;
  387.         }
  388.     this->UpdateFlds();    
  389.     this->Modified();
  390. }
  391.  
  392.  
  393. DSequence* DMegaSequence::PasteBases( long insertPoint, char* fromBases, char* fromMasks)
  394. { 
  395.     long add;
  396.     DMegaSequence* aSeq= (DMegaSequence*) Clone();
  397.     if (fromBases) add= StrLen( fromBases);  else add= 0;
  398.     if (add>0) {
  399.         char* h= aSeq->fBases;
  400.         char* m= aSeq->fMasks;
  401.         Boolean domask= (fMasks && fMasksOk);
  402.         long len = aSeq->fLength;
  403.         if (insertPoint<=0) {
  404.             fromBases= StrDup( fromBases);
  405.             StrExtendCat( &fromBases, h);
  406.             MemFree( h);
  407.             h= fromBases;
  408.             if (domask && fromMasks) {
  409.                 fromMasks= StrDup( fromMasks);
  410.                 StrExtendCat( &fromMasks, m);
  411.                 MemFree( m);
  412.                 m= fromMasks;
  413.                 }
  414.             }        
  415.         else if (insertPoint>=len) { 
  416.             StrExtendCat( &h, fromBases);
  417.             if (domask) StrExtendCat( &m, fromMasks);
  418.             }
  419.         else { 
  420.             h= (char*) MemMore( h, len+add);
  421.             MemMove( h+insertPoint+add, h+insertPoint, len-insertPoint);
  422.             MemMove( h+insertPoint+add, fromBases, add);
  423.             *((char*)h+len+add)= '\0';
  424.             if (domask) {
  425.                 m= (char*) MemMore( m, len+add);
  426.                 MemMove( m+insertPoint+add, m+insertPoint, len-insertPoint);
  427.                 MemMove( m+insertPoint+add, fromMasks, add);
  428.                 *((char*)m+len+add)= '\0';
  429.                 }
  430.             }
  431.         aSeq->fBases= h;
  432.         aSeq->fMasks= m;
  433.         }
  434.     aSeq->UpdateFlds();    
  435.     aSeq->Modified();
  436.     return aSeq;
  437. }
  438.  
  439.  
  440. DSequence* DMegaSequence::Reverse()
  441. {
  442.     long i;
  443.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  444.     if (fSelStart < 0) fSelStart= 0;
  445.     long len= aSeq->fLength - fSelStart;
  446.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  447.     char* pc= this->fBases+fSelStart + len;
  448.     char* pr= aSeq->fBases+fSelStart;
  449.     for (i= 0; i<len; i++) *pr++= *--pc;
  450.     Boolean domask= (fMasks && fMasksOk);
  451.     if (domask) {
  452.         pc= this->fMasks+fSelStart + len;
  453.         pr= aSeq->fMasks+fSelStart;
  454.         for (i= 0; i<len; i++) *pr++= *--pc;
  455.         }
  456.     aSeq->UpdateFlds();
  457.     aSeq->Modified();
  458.     return aSeq;
  459. }
  460.  
  461.  
  462. DSequence* DMegaSequence::Complement()
  463. // A->T, T->A, G->C, C->G  
  464. {
  465.     Boolean        isrna;
  466.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  467.     
  468.     if (fKind == kRNA) isrna= TRUE;
  469.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  470.     else return aSeq;    
  471.         
  472.     if (fSelStart < 0) fSelStart= 0;
  473.     long len= aSeq->fLength - fSelStart;
  474.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  475.     char* pc= this->fBases+fSelStart;
  476.     char* pr= aSeq->fBases+fSelStart;
  477.     NucleicComplement( isrna, pc, pr, len);
  478.         
  479.     aSeq->UpdateFlds();
  480.     aSeq->Modified();
  481.     return aSeq;
  482. }
  483.  
  484.  
  485.  
  486. DSequence* DMegaSequence::ChangeCase(Boolean toUpper)
  487. // g->G or G->g
  488. {
  489.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  490.     
  491.     if (fSelStart < 0) fSelStart= 0;
  492.     char* h= aSeq->fBases;
  493.     long len= aSeq->fLength - fSelStart;
  494.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  495.     char* pc= this->fBases+fSelStart;
  496.     char* pr= aSeq->fBases+fSelStart;
  497.     
  498.     long i;
  499.     if (toUpper) for (i= 0; i<len; i++) {
  500.         *pr++= toupper(*pc); pc++;
  501.         }
  502.     else for (i= 0; i<len; i++) {
  503.         *pr++= tolower(*pc); pc++;
  504.         }    
  505.         
  506.     aSeq->UpdateFlds();
  507.     aSeq->Modified();
  508.     return aSeq;
  509. }
  510.  
  511. DSequence* DMegaSequence::Dna2Rna(Boolean toRna)
  512. // T->U or U->T 
  513. {
  514.     Boolean        isrna;
  515.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  516.     
  517.     if (fKind==kRNA) isrna= TRUE;
  518.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  519.     else return aSeq;    
  520.      
  521.     if (fSelStart < 0) fSelStart= 0;
  522.     long len= aSeq->fLength - fSelStart;
  523.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  524.     char* pc= this->fBases+fSelStart;
  525.     char* pr= aSeq->fBases+fSelStart;
  526.     
  527.     long i;
  528.     char b, c;
  529.     if (toRna) for (i= 0; i<len; i++) {
  530.         b= *pc++;
  531.         if (b=='T') c= 'U';
  532.         else if (b=='t') c= 'u';
  533.         else c= b;
  534.         *pr++= c;
  535.         }
  536.     else for (i= 0; i<len; i++) {
  537.         b= *pc++;
  538.         if (b=='U') c= 'T';
  539.         else if (b=='u') c= 't';
  540.         else c= b;
  541.         *pr++= c;
  542.         }    
  543.         
  544.     aSeq->UpdateFlds();
  545.     aSeq->Modified();
  546.     return aSeq;
  547. }
  548.  
  549.  
  550. DSequence* DMegaSequence::OnlyORF(char nonorf)
  551. // convert sequence to ORF part only 
  552. {
  553.     return DSequence::OnlyORF(nonorf);    
  554. }
  555.  
  556. DSequence* DMegaSequence::LockIndels( Boolean doLock)
  557. // indelSoft -> indelHard or vice versa 
  558. {    
  559.     char b, c, fromc, toc;
  560.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  561.     
  562.     if (fSelStart < 0) fSelStart= 0;
  563.     long  len= aSeq->fLength - fSelStart;
  564.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  565.     char* pc= this->fBases+fSelStart;
  566.     char* pr= aSeq->fBases+fSelStart;
  567.     
  568.     if (doLock) { fromc= indelSoft; toc= indelHard; }
  569.     else { fromc= indelHard; toc= indelSoft; }
  570.     
  571.     for (long i= 0; i<len; i++) {
  572.         b= *pc++;
  573.         if (b == fromc) c= toc;
  574.         else c= b;
  575.         *pr++= c;
  576.         }    
  577.         
  578.     aSeq->UpdateFlds();
  579.     aSeq->Modified();
  580.     return aSeq;
  581. }
  582.  
  583.  
  584. DSequence* DMegaSequence::Slide( long slideDist)
  585. /* insert "-" below (+dist) or above (-dist) of seq range,
  586.     squeeze out "-" above (below) of range
  587. */
  588. {     
  589.     long         slideSave, newlen, fulllen, len, i, j= 0;
  590.     char         b, cutchar, cutmask;
  591.     char        *pc, *pr;
  592.     
  593.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  594.     if (fSelStart < 0) fSelStart= 0;
  595.     fulllen= aSeq->fLength;
  596.     len= fulllen - fSelStart;
  597.     Boolean domask= (fMasks && fMasksOk);
  598.     Boolean dosqueeze= true; //(indelSoft & 0x80 != 0); // quick fix flag
  599.     indelSoft &= 0x7f;
  600.     
  601.     if (slideDist < 0) {
  602.         // tough... do if below works...
  603.         // really need ability to extend seq in 5' (lower) direction 
  604.         //! build pr in reverse, then flip it....
  605.  
  606.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  607.         slideDist= -slideDist;
  608.         slideSave= slideDist;
  609.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fulllen+slideDist);
  610.         pr= aSeq->fBases;
  611.         pc= this->fBases + fulllen;
  612.         for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  613.              
  614.         if (pr == aSeq->fBases || *(pr-1) == indelEdge) cutchar= indelEdge;
  615.         else cutchar= indelSoft; //??
  616.         for (i=0; i<slideDist; i++) *pr++= cutchar;
  617.  
  618.         pc= this->fBases+fSelStart+len;
  619.         for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  620.             
  621.         pc= this->fBases+fSelStart;
  622.         for (i= fSelStart-1; i>=0; i--) {
  623.             b= *--pc;
  624.             if (dosqueeze && slideDist>0 && (b==indelSoft || b==indelEdge)) slideDist--;
  625.             else {
  626.                 *pr++= b;
  627.                 }
  628.             }
  629.  
  630.         newlen= pr - aSeq->fBases;
  631.         char* htmp= (char*) MemNew(newlen+1);
  632.         pc= htmp;
  633.         for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  634.         pc= '\0';
  635.         MemFree(aSeq->fBases);
  636.         aSeq->fBases = htmp;  
  637.         aSeq->fLength= newlen;
  638.         aSeq->fOrigin= fulllen - newlen; //== this->fLength - aSeq->fLength
  639.  
  640.         if (domask) {
  641.             slideDist= slideSave;
  642.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fulllen+slideDist);
  643.             pr= aSeq->fMasks;
  644.             pc= this->fMasks + fulllen;
  645.             for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  646.                  
  647.             cutmask= 0x80; // required empty mask  
  648.             for (i=0; i<slideDist; i++) *pr++= cutmask;
  649.     
  650.             pc= this->fMasks+fSelStart+len;
  651. #ifdef WIN_MSWIN
  652.             // mswin is crashing here for no obvious reason
  653.             //pc--; len--;
  654.             for (i= fSelStart + len - 1; i >= fSelStart; i--) {
  655.                  register short amask;
  656.                  --pc;
  657.                  amask= *pc;    //mswin bomb is here !
  658.                  *pr= amask;
  659.                  pr++;
  660.          }
  661. #else
  662.             for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  663. #endif
  664.             pc= this->fMasks+fSelStart;
  665.             for (i= fSelStart-1; i>=0; i--) {
  666.                 b= *--pc;
  667.                 if (slideDist>0 && (b==cutmask)) slideDist--;
  668.                 else {
  669.                     *pr++= b;
  670.                     }
  671.                 }
  672.                 
  673.             newlen= pr - aSeq->fMasks;
  674.             htmp= (char*) MemNew(newlen+1);
  675.             pc= htmp;
  676.             for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  677.             pc= '\0';
  678.             MemFree(aSeq->fMasks);
  679.             aSeq->fMasks = htmp;  
  680.             }
  681.             
  682.         }    
  683.         
  684.     else if (slideDist > 0) {
  685.         fulllen= len;
  686.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  687.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fSelStart+fulllen+slideDist); //make it large enough
  688.  
  689.         pr= aSeq->fBases;
  690.       if (fSelStart==0 || pr[fSelStart-1]==indelEdge) cutchar= indelEdge; 
  691.         else cutchar= indelSoft; 
  692.  
  693.         pc= this->fBases+fSelStart;
  694.         pr= aSeq->fBases+fSelStart;
  695.  
  696.         for (j= 0; j<slideDist; j++) *pr++= cutchar;
  697.         for (i= 0; i<len; i++) *pr++= *pc++; 
  698.     
  699.         for (i= len; i<fulllen; i++) {
  700.             b= *pc++;
  701.             if (dosqueeze && slideDist>0 && (b==indelSoft || b==indelEdge)) 
  702.                 slideDist--;
  703.             else {
  704.                 *pr++= b;
  705.                 }
  706.             }
  707.         
  708.         *pr= 0;
  709.         newlen= pr - aSeq->fBases; 
  710.         aSeq->fBases= (char*)MemMore(aSeq->fBases, newlen+1);  
  711.         
  712.         if (domask) {
  713.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fSelStart+fulllen+slideDist);
  714.           cutmask= 0x80; 
  715.             pc= this->fMasks+fSelStart;
  716.             pr= aSeq->fMasks+fSelStart;
  717.     
  718.             for (j= 0; j<slideDist; j++) *pr++= cutmask;
  719.             for (i= 0; i<len; i++) *pr++= *pc++; 
  720.  
  721.             for (i= len; i<fulllen; i++) {
  722.                 b= *pc++;
  723.                 if (slideDist>0 && (b==cutmask)) // this is BAD? must be done w/ base slide !?
  724.                     slideDist--;
  725.                 else {
  726.                     *pr++= b;
  727.                     }
  728.                 }
  729.             
  730.             *pr= 0;
  731.             newlen= pr - aSeq->fMasks; 
  732.             aSeq->fMasks= (char*)MemMore(aSeq->fMasks, newlen+1);  
  733.             }
  734.         }
  735.             
  736.     aSeq->UpdateFlds();
  737.     aSeq->Modified();
  738.     return aSeq;
  739. }
  740.  
  741.  
  742. DSequence* DMegaSequence::Compress()
  743. // pull out "-", ".", ? anything but nuc/aa codes ? 
  744. {
  745.     long  len, fulllen, i;
  746.     char     b, c;
  747.     char * mc, * mr;
  748.     Boolean domask= (fMasks && fMasksOk);
  749.  
  750.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  751.     if (fSelStart < 0) fSelStart= 0;
  752.     len= aSeq->fLength - fSelStart;
  753.     fulllen= len;
  754.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  755.     char* pc= this->fBases+fSelStart;
  756.     char* pr= aSeq->fBases+fSelStart;
  757.     if (domask) {
  758.         mc= this->fMasks+fSelStart;
  759.         mr= aSeq->fMasks+fSelStart;
  760.         }
  761.         
  762.     if (domask) for (i= 0; i<len; i++) {
  763.         b= *pc++;
  764.         c= *mc++;
  765.         if (b!=indelSoft && b!=indelEdge) { 
  766.             *pr++= b;
  767.             *mr++= c;
  768.             }
  769.         }
  770.     else for (i= 0; i<len; i++) {
  771.         b= *pc++;
  772.         if (b!=indelSoft && b!=indelEdge) 
  773.             *pr++= b;
  774.         }
  775.     for (i=len; i<fulllen; i++) *pr++= *pc++;
  776.     if (domask) for (i=len; i<fulllen; i++) *mr++= *mc++;
  777.  
  778.     long newlen= pr - aSeq->fBases;
  779.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen+1);
  780.     aSeq->fBases[newlen]= 0;
  781.     if (domask) {
  782.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen+1);
  783.         aSeq->fMasks[newlen]= 0;
  784.         }
  785.         
  786.     aSeq->UpdateFlds();
  787.     aSeq->Modified();
  788.     return aSeq;
  789. }
  790.  
  791.  
  792.  
  793.  
  794. DSequence* DMegaSequence::CompressFromMask(short masklevel)
  795. {
  796.     long  len, fulllen, i;
  797.     char     b, c;
  798.     char * mc, * mr;
  799.     Boolean domask= (fMasks && fMasksOk && masklevel>0);
  800.  
  801.     if (!domask) return NULL; // or aSeq ???
  802.     DMegaSequence* aSeq= (DMegaSequence*) Clone();    
  803.     //if (!domask) return aSeq;  
  804.  
  805.     if (fSelStart < 0) fSelStart= 0;
  806.     len= aSeq->fLength - fSelStart;
  807.     fulllen= len;
  808.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  809.     char* pc= this->fBases+fSelStart;
  810.     char* pr= aSeq->fBases+fSelStart;
  811.     mc= this->fMasks+fSelStart;
  812.     mr= aSeq->fMasks+fSelStart;
  813.          
  814.     for (i= 0; i<len; i++) {
  815.         b= *pc++;
  816.         c= *mc++;
  817.         if (IsMasked(c, masklevel)) { 
  818.             *pr++= b;
  819.             *mr++= c;
  820.             }
  821.         }
  822.     for (i=len; i<fulllen; i++) *pr++= *pc++;
  823.     for (i=len; i<fulllen; i++) *mr++= *mc++;
  824.  
  825.     long newlen= pr - aSeq->fBases;
  826.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen+1);
  827.     aSeq->fBases[newlen]= 0;
  828.     aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen+1);
  829.     aSeq->fMasks[newlen]= 0;
  830.         
  831.     aSeq->UpdateFlds();
  832.     aSeq->Modified();
  833.     return aSeq;
  834. }
  835.  
  836.  
  837.  
  838. DSequence* DMegaSequence::Translate(Boolean keepUnselected)
  839. {    
  840.     char  c;
  841.     char    acodon[3];
  842.     Boolean isrna, nocodon;
  843.     char    *pr, *pc, *mr, *mc;
  844.     long    len, fulllen, newlen, i, j, k;
  845.     
  846.     DMegaSequence* aSeq= (DMegaSequence*) Clone();        
  847.     if (DCodons::NotAvailable()) return aSeq; // safer than returning NULL !?
  848.     // !? what do we do w/ fMasks !?
  849.     Boolean domask= (fMasks && fMasksOk);
  850.     
  851.     if (fSelStart < 0) fSelStart= 0;
  852.     fulllen= this->fLength - fSelStart;
  853.     len= fulllen;
  854.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  855.     
  856.     switch (fKind) {
  857.     
  858.         case kRNA:
  859.         case kDNA:
  860.         case kNucleic: 
  861.             {
  862.             if (keepUnselected) {
  863.                 pr= aSeq->fBases + fSelStart; // start at selection, leaving untranslated 5'
  864.                 if (domask) mr= aSeq->fMasks + fSelStart;
  865.                 }
  866.             else {
  867.                 pr= aSeq->fBases;
  868.                  if (domask) mr= mr= aSeq->fMasks;
  869.                  }
  870.                  
  871.             pc= this->fBases + fSelStart;
  872.             if (domask) mc= this->fMasks + fSelStart;
  873.             isrna= (fKind == kRNA);
  874.             i= fSelStart;
  875.             while (i < len) {
  876.                 for (j=0; j<3; j++) {
  877.                     c= toupper( *pc); pc++; // toupper is sometimes a MACRO ! no (*pc++) !
  878.                     if (isrna && c=='U') c= 'T';
  879.                     acodon[j]= c;
  880.                     }
  881.                 i += 3;
  882.                 for (k= 0, nocodon= true; k<64; k++) 
  883.                     if ( MemCmp(DCodons::codon(k),acodon,3) == 0) { 
  884.                         *pr++= DCodons::amino(k); //gCodonTable[k].amino;
  885.                         nocodon= false;
  886.                         break;
  887.                         }
  888.                 if (nocodon) *pr++= '?';
  889.                 if (domask) {
  890.                     *mr++= *mc;
  891.                     mc += 3;
  892.                     }
  893.                 }
  894.                  // copy over untranslated portion
  895.             if (keepUnselected) {
  896.                 long savei= i;
  897.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  898.                 if (domask) for ( i= savei; i < fulllen; i++)  *mr++= *mc++; 
  899.                 }
  900.             }
  901.             break;
  902.             
  903.         case kAmino: 
  904.             {
  905.             char *bestcodon;
  906.             long  selend;
  907.             
  908.             if (len < fulllen) {
  909.                 newlen= fulllen + len*3;
  910.                 selend= fSelStart + len;
  911.                 }
  912.             else {
  913.                 newlen= fulllen * 3;
  914.                 selend= len;
  915.                 }
  916.             pc= this->fBases + fSelStart;
  917.             aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  918.             if (keepUnselected)
  919.                 pr= aSeq->fBases + fSelStart;
  920.             else
  921.                 pr= aSeq->fBases;
  922.                 
  923.             if (domask) {
  924.                 mc= this->fMasks + fSelStart;
  925.                 aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  926.                 if (keepUnselected)
  927.                     mr= aSeq->fMasks + fSelStart;
  928.                 else
  929.                     mr= aSeq->fMasks;
  930.                 }
  931.                 
  932.             for (i= fSelStart; i<selend; i++) {
  933.                 c= toupper( *pc); pc++;// toupper is sometimes a MACRO ! no (*pc++) !
  934.                 bestcodon= (char*) DCodons::FindBestCodon(c);
  935.                 MemCpy( pr, bestcodon, 3);  
  936.                 pr += 3;
  937.                 if (domask) {
  938.                     for (j=0; j<3; j++) *mr++= *mc;
  939.                     mc++;
  940.                     }
  941.                 }
  942.                  // copy over untranslated portion
  943.             if (keepUnselected) {
  944.                 long savei= i;
  945.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  946.                 if (domask) for ( i=savei; i < fulllen; i++)  *mr++= *mc++; 
  947.                 }
  948.             }
  949.             break;
  950.         
  951.         default:
  952.             return NULL;
  953.          
  954.       }
  955.         
  956.     *pr= '\0';
  957.     newlen= pr - aSeq->fBases;
  958.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  959.     if (domask) {
  960.         *mr= '\0';
  961.         newlen= mr - aSeq->fMasks;
  962.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  963.         }
  964.     aSeq->UpdateFlds();
  965.     aSeq->Modified();
  966.     return aSeq;    
  967. }
  968.  
  969.  
  970.  
  971.  
  972. void DMegaSequence::DoWrite(DFile* aFile, short format)
  973. {  
  974.     if (fLength>0) {
  975.         DSeqFile::WriteSeqWrapper( aFile, fBases, fLength, format, fName);
  976.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks) 
  977.             DSeqFile::WriteMaskWrapper( aFile, fMasks, fLength, format, fName);
  978.         }
  979. }
  980.  
  981.  
  982. void DMegaSequence::DoWriteSelection(DFile* aFile, short format)
  983. {
  984.     if (fSelBases==0) 
  985.         DoWrite( aFile, format);
  986.     else if (fLength>0) {
  987.         long aStart= Min( fSelStart, fLength);
  988.         long aLength= Min( fSelBases, fLength - fSelStart);
  989.         DSeqFile::WriteSeqWrapper(  aFile, fBases+aStart, aLength, format, fName);
  990.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks)  
  991.             DSeqFile::WriteMaskWrapper( aFile, fMasks+aStart, aLength, format, fName);
  992.         }
  993. }
  994.  
  995.  
  996.  
  997.  
  998.  
  999.  
  1000.  
  1001.  
  1002.  
  1003. void DMegaSequence::SearchORF( long& start, long& stop) 
  1004. {
  1005.     char  startc[2];
  1006.     long  i, j, k;
  1007.     short    ns;
  1008.     CodonStat  stops[10];
  1009.     
  1010.     start= kSearchNotFound;
  1011.     stop = kSearchNotFound;
  1012.     if (DCodons::NotAvailable()) return;
  1013.     
  1014.     for (ns= 0, i= 0; i<64; i++) {
  1015.         if (ns<10 && DCodons::amino(i) == '*') 
  1016.             stops[ns++] = DCodons::fCodons[i];
  1017.         }
  1018.     if (ns == 0) return;
  1019.         
  1020.     switch (fKind) {  
  1021.         case kRNA:
  1022.         case kDNA:
  1023.         case kNucleic: 
  1024.              start= Search( DCodons::startcodon(), FALSE);
  1025.              if (start == kSearchNotFound) return;
  1026.              for (j= start+3; j<fLength; j += 3) {
  1027.                  for (k=0; k<ns; k++) {
  1028.                      if (StrNICmp( fBases+j, stops[k].codon, 3) == 0) {
  1029.                          stop= j+3; // +3 to include this codon
  1030.                          return;
  1031.                          }
  1032.                      }
  1033.                  }
  1034.             break;
  1035.             
  1036.         case kAmino:
  1037.             startc[0]= DCodons::startamino();
  1038.             startc[1]= 0;
  1039.              start= Search( startc, FALSE);
  1040.              if (start == kSearchNotFound) return;
  1041.              for (j= start+1; j<fLength; j++) {
  1042.                  for (k=0; k<ns; k++) {
  1043.                      if (toupper(fBases[j]) == toupper(stops[k].amino)) {
  1044.                          stop= j+1; // +1 to include this codon
  1045.                          return;
  1046.                          }
  1047.                      }
  1048.                  }
  1049.             break;
  1050.             
  1051.         default:
  1052.             break;
  1053.         }
  1054. }
  1055.  
  1056.  
  1057.  
  1058. void DMegaSequence::NucleicSearch(char* source, char* target, SearchRec& aSearchRec)
  1059. {
  1060.     long k, j, tBits, lens, ti;
  1061.     Boolean done;
  1062.  
  1063.                     //find target[ti] where tBits == kMaskA,C,G,T for fastest search
  1064.                     //This part is wasted time for SearchAgain... store tBits/ti in saverec...
  1065.     lens= StrLen(target);
  1066.     for (j= 0, ti= 0, tBits= 0; j<lens && !tBits; j++) {
  1067.         ti= j;
  1068.         tBits= NucleicBits( target[ti]);
  1069.         tBits &= kMaskNucs; //    bClr( tBits, kBitExtra); //drop RNA bit
  1070.         }
  1071.  
  1072.     do {
  1073. reloop:
  1074.             aSearchRec.indx += aSearchRec.step; 
  1075.             aSearchRec.lim  -= aSearchRec.step;    
  1076.             
  1077.             if (aSearchRec.step>0) {
  1078.                 for (j= 0; j<aSearchRec.lim; j++)
  1079.                  // if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx]))) 
  1080.                  if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1081.                  {
  1082.                     k= j; 
  1083.                     goto nextindx;
  1084.                  }
  1085.                 k= aSearchRec.lim;
  1086.                 }            
  1087.             else {  //step<0 IS BAD .. see UTextDoc version...
  1088.                 for (j= 0; j > -aSearchRec.lim; j--) 
  1089.                   //if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx])))
  1090.                   if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1091.                 {
  1092.                     k= j; 
  1093.                     goto nextindx;
  1094.                     }
  1095.                 k= -aSearchRec.lim;
  1096.                 }
  1097.                 
  1098. nextindx:  
  1099.             aSearchRec.indx += k; 
  1100.             aSearchRec.lim  -= k;
  1101.                 //? not enough source left to match full target 
  1102.             if (aSearchRec.step>0) {
  1103.                 if (aSearchRec.lim < lens-ti) aSearchRec.lim= 0;
  1104.                 }            
  1105.             else {
  1106.                 if (aSearchRec.lim < ti) aSearchRec.lim= 0; //? or lim < ti-1 
  1107.                 }
  1108.             done= (aSearchRec.lim==0);
  1109.             
  1110.             if (!done) {
  1111.                 long srcstart= aSearchRec.indx - ti;
  1112.                 for (j= 0; j<lens; j++) 
  1113.                  if (j != ti) {
  1114.                     //if (!IsNucleicMatch( NucleicBits( target[j]), NucleicBits( source[j+srcstart]))) 
  1115.                   if (!(kMaskNucs & NucleicBits( target[j]) & NucleicBits( source[j+srcstart]))) 
  1116.                       goto reloop;
  1117.                     }
  1118.                 done= true;
  1119.                 }
  1120.         } while (!done);
  1121. }  
  1122.     
  1123.  
  1124. void DMegaSequence::ProteinSearch(char* source, char* target, SearchRec& aSearchRec)
  1125. {
  1126.     //just toupper match of chars 
  1127.     char      tChar; 
  1128.     Boolean  done;
  1129.     long        lens, k, j;
  1130.  
  1131.     lens= StrLen( target);
  1132.     tChar= toupper( target[0]);
  1133.     
  1134.     do {
  1135. reloop:
  1136.             aSearchRec.indx += aSearchRec.step; 
  1137.             aSearchRec.lim  -= aSearchRec.step;    
  1138.             
  1139.             if (aSearchRec.step>0) {
  1140.                 for (j= 0; j<aSearchRec.lim; j++) 
  1141.                  if (tChar == toupper( source[j+aSearchRec.indx])) {
  1142.                     k= j;
  1143.                     goto nextindx;
  1144.                     }
  1145.                 k= aSearchRec.lim;
  1146.                 }            
  1147.             else {
  1148.                 for (j= 0; j > -aSearchRec.lim; j--)
  1149.                  if (tChar == toupper(source[j+aSearchRec.indx])) {
  1150.                     k= j; 
  1151.                     goto nextindx;
  1152.                     }
  1153.                 k= -aSearchRec.lim;
  1154.                 }
  1155.                 
  1156. nextindx:  
  1157.             aSearchRec.indx += k; 
  1158.             aSearchRec.lim  -= k;
  1159.                 //? not enough source left to match full target 
  1160.             if (aSearchRec.step>0) {
  1161.                 if (aSearchRec.lim < lens-1) aSearchRec.lim= 0;
  1162.                 }            
  1163.             else {
  1164.                 if (aSearchRec.lim < 1) aSearchRec.lim= 0;  
  1165.                 }
  1166.             done= (aSearchRec.lim==0);
  1167.             
  1168.             if (!done) {
  1169.                 for (j= 1; j<lens; j++) {
  1170.                     if (toupper( target[j]) != toupper( source[j+aSearchRec.indx])) 
  1171.                         goto reloop;
  1172.                     }
  1173.                 done= true;
  1174.                 }
  1175.         } while (!done);
  1176. }  
  1177.  
  1178.  
  1179.     
  1180.  
  1181. long DMegaSequence::Search( char* target, Boolean backwards)
  1182. /* search for target in sequence;
  1183.     search == abs index into fBases[0..fLength] or kSearchNotFound
  1184.   if (target='') then find next  
  1185.   if (backwards) then backwards search
  1186. */
  1187. {
  1188.     long        limit, aStart;
  1189.     SearchRec  aSearchRec;
  1190.  
  1191.     aStart= fSelStart;
  1192.     if (!target || !*target) {  //search again
  1193.         if (fSearchRec.targ && (fSearchRec.step == -1 || fSearchRec.step == 1)) {  
  1194.             aSearchRec= fSearchRec;
  1195.             target= aSearchRec.targ;
  1196.             }        
  1197.         else {
  1198.             return kSearchNotFound;
  1199.             }
  1200.         }        
  1201.     else {
  1202.         limit= fLength; //StrLen(fBases);
  1203.         if (aStart < 0) aStart= 0;
  1204.         if (backwards) {
  1205.             if (!aStart) aStart= limit-1;
  1206.             else limit= aStart+1;  //? off-by-one ??
  1207.             }        
  1208.         else {
  1209.             if (aStart > limit) aStart= 0; 
  1210.             else limit -= aStart;
  1211.             }
  1212.         if (fSelBases > 0 && fSelBases < limit) limit= fSelBases;
  1213.         if (backwards) aSearchRec.step= -1;
  1214.         else aSearchRec.step= 1;
  1215.         aSearchRec.indx= aStart - aSearchRec.step;  //indx always in [0..handsize-1], except 1st is -1 
  1216.         aSearchRec.lim = limit + aSearchRec.step;  //lim always in [1..handsize] or 0=done, except 1st is hand+1 
  1217.         aSearchRec.targ= StrDup(target);
  1218.         }
  1219.  
  1220.     switch (fKind) {  
  1221.         case kRNA:
  1222.         case kDNA:
  1223.         case kNucleic: 
  1224.             NucleicSearch(fBases,target,aSearchRec); 
  1225.             break;
  1226.         default:
  1227.             ProteinSearch(fBases,target,aSearchRec);
  1228.             break;
  1229.         }
  1230.     
  1231.     fSearchRec= aSearchRec;
  1232.     if (aSearchRec.lim == 0) return kSearchNotFound;
  1233.     else return aSearchRec.indx;
  1234. }  
  1235.  
  1236.  
  1237.