home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / DSeqList.cpp < prev    next >
Encoding:
Text File  |  1996-07-05  |  14.1 KB  |  621 lines  |  [TEXT/R*ch]
  1. // DSeqList.cp
  2. // d.g.gilbert
  3.  
  4.  
  5. #include "DSeqList.h"
  6. #include "DSequence.h"
  7. #include "DSeqFile.h"
  8. #include <ncbi.h>
  9. #include <dgg.h>
  10. #include <DList.h>
  11. #include <DFile.h>
  12. #include <ureadseq.h>
  13.  
  14.  
  15.  
  16. Global short gLinesPerSeqWritten = 0;
  17. short DSeqList::gMinCommonPercent = 70;
  18. short DSeqList::gMinORFsize = 20;
  19.  
  20.  
  21. // DSeqList  -------------------------------------
  22.  
  23. DSeqList::DSeqList() :
  24.     DList(),
  25.     fSortOrder(kSortByItem)
  26. {
  27. }
  28.  
  29. DSeqList::~DSeqList()
  30. {
  31.         // PLUG MEMORY LEAK !! ~DList doesn't call object destructors !!!!
  32.     if (fDeleteObjects) {
  33.         long i, n= GetSize();
  34.         for (i=0; i<n; i++) {
  35.             DSequence* aseq= this->SeqAt(i);
  36.             delete aseq;
  37.             }
  38.         fDeleteObjects= false;
  39.         }
  40. }
  41.  
  42. Boolean DSeqList::IsEmpty()
  43. {
  44.     if (GetSize() > 1) 
  45.         return FALSE;
  46.     else if (GetSize() == 1) {
  47.         DSequence* aSeq= (DSequence*) this->First();
  48.         return (aSeq && aSeq->LengthF() > 0);
  49.         }    
  50.     else
  51.       return TRUE;
  52. }
  53.  
  54.  
  55.  
  56.  
  57.  
  58. class CFileIndex 
  59. {
  60. public:
  61.     long*    fIndices;
  62.     long    fMax, fNum;
  63.  
  64.     CFileIndex();
  65.     ~CFileIndex();
  66.      long* Indices() { return fIndices; }
  67.      long    IndexCount() { return fNum; }
  68.     void Indexit( long index);
  69.     void Indexit( DFile* aFile);
  70.     void IndexEOF( DFile* aFile);
  71. };
  72.  
  73. CFileIndex::CFileIndex() 
  74. {
  75.     fNum= 0;
  76.     fMax= 0;
  77.     fIndices= (long*) MemGet(sizeof(long),true);
  78. }
  79.  
  80. CFileIndex::~CFileIndex()
  81. {
  82.     if (fIndices) MemFree( fIndices);
  83. }
  84.      
  85. void CFileIndex::Indexit( long index)
  86. {
  87.     if (fNum >= fMax){
  88.         fMax += 20; 
  89.         fIndices= (long*) MemMore( fIndices, sizeof(long)*fMax); 
  90.         }
  91.     if (fIndices) fIndices[fNum++]= index;
  92. }
  93.  
  94. void CFileIndex::Indexit( DFile* aFile)
  95. {
  96.     Indexit(aFile->Tell());
  97. }
  98.  
  99. void CFileIndex::IndexEOF( DFile* aFile)
  100. {
  101.     Indexit(aFile->Tell());
  102.     fNum--;
  103. }
  104.     
  105.  
  106.  
  107.     
  108.  
  109. void DSeqList::DoWrite( DFile* aFile, short format)  // change to iostreams !!
  110. {
  111.     Boolean        needSameSize = false, sizesDiffer = false, isInterleaved = false;
  112.     short         iseq, nseqs;
  113.     short            seqtype = DSequence::kOtherSeq;
  114.     long          minbases = -1;
  115.     long             start, nbases;
  116.     DSequence * aSeq;
  117.     DFile         * outfile;
  118.     DFile            * tmpfile = NULL;
  119.     CFileIndex aFileIndex;
  120.  
  121.     
  122.     nseqs= this->GetSize();
  123.     outfile= aFile;
  124.     
  125.     if (nseqs>1) // deal w/ one-per-file formats 
  126.         switch (format) {
  127.             case DSeqFile::kGCG            : format= DSeqFile::kMSF; break;
  128.             case DSeqFile::kStrider    : format= DSeqFile::kIG; break;
  129.             case DSeqFile::kNoformat:
  130.             case DSeqFile::kPlain:
  131.             case DSeqFile::kUnknown : format= DSeqFile::kGenBank; break;
  132.             default                :    break;
  133.             }
  134.         
  135.     //this->Each( findMinSize);    
  136.     for (iseq= 0; iseq<nseqs; iseq++) {
  137.         aSeq= this->SeqAt(iseq);
  138.         aSeq->GetSelection( start, nbases);
  139.         seqtype= aSeq->Kind();
  140.         //if (minbases<0) minbases= nbases;
  141.         if (minbases<=0) minbases= nbases; // ?? temp fix problem w/ empty seq
  142.         else {
  143.             if (nbases!=minbases) sizesDiffer= true;
  144.             minbases= Min( nbases, minbases);
  145.             }
  146.         }
  147.         
  148.     if (!gOutputName) gOutputName= (char*)aFile->GetShortname();
  149.     DSeqFile::WriteSeqHeader( aFile, format, nseqs, minbases, gOutputName, needSameSize, isInterleaved);
  150.     gOutputName= NULL;
  151.     
  152.     if (isInterleaved) {
  153.         tmpfile= new DTempFile();
  154.         outfile= tmpfile;
  155.         }
  156.     
  157.     if ((needSameSize || isInterleaved) && sizesDiffer) {
  158.         //this->Each( SetSameSize);
  159.         for (iseq= 0; iseq<nseqs; iseq++) {
  160.             aSeq= this->SeqAt(iseq);
  161.             aSeq->GetSelection( start, nbases);
  162.             aSeq->SetSelection( start, minbases); 
  163.             }
  164.         }
  165.  
  166.     if (gPretty.isactive && gPretty.numtop) {
  167.         aSeq= this->SeqAt(nseqs-1);
  168.         gPretty.numline = 1;
  169.         if (isInterleaved) aFileIndex.Indexit(outfile);
  170.         aSeq->DoWriteSelection( outfile, format);
  171.         gPretty.numline = 2;
  172.         if (isInterleaved) aFileIndex.Indexit(outfile);
  173.         aSeq->DoWriteSelection( outfile, format);
  174.         gPretty.numline = 0;
  175.         }
  176.  
  177.     for (iseq= 0; iseq<nseqs; iseq++) {
  178.         aSeq= this->SeqAt(iseq);
  179.         if (isInterleaved) aFileIndex.Indexit(outfile);
  180.         aSeq->DoWriteSelection( outfile, format);
  181.         }
  182.                 
  183.     if (gPretty.isactive && gPretty.numbot) {
  184.         aSeq= this->SeqAt(nseqs-1);
  185.         gPretty.numline = 2;
  186.         if (isInterleaved) aFileIndex.Indexit(outfile);
  187.         aSeq->DoWriteSelection( outfile, format);
  188.         gPretty.numline = 1;
  189.         if (isInterleaved) aFileIndex.Indexit(outfile);
  190.         aSeq->DoWriteSelection( outfile, format);
  191.         gPretty.numline = 0;
  192.         }
  193.  
  194.     if (isInterleaved) {
  195.         aFileIndex.IndexEOF( outfile);
  196.         //tmpfile->Close();   // don't need
  197.         //tmpfile->Open("r"); // don't need
  198.         DSeqFile::DoInterleave( aFile, tmpfile, format, seqtype, nseqs,
  199.                                     minbases, gLinesPerSeqWritten, 
  200.                                     aFileIndex.IndexCount(), aFileIndex.Indices());
  201.         outfile= aFile;
  202.         //tmpfile->Close(); // don't need
  203.         tmpfile->Delete();  // debug tmp
  204.         }
  205.         
  206.     DSeqFile::WriteSeqTrailer( aFile, format, this->GetSize(), gLinesPerSeqWritten, minbases);
  207. }
  208.  
  209.  
  210.  
  211. void DSeqList::DoWrite( char* aFileName, short format)  
  212. {
  213.     gOutputName= aFileName;
  214.     DFile* aFile= new DFile( aFileName, "w"); //, "TEXT", fgSire);
  215.     aFile->Open();
  216.     DoWrite( aFile, format);
  217.     aFile->Close();
  218. }
  219.  
  220.  
  221.  
  222. void DSeqList::AddNewSeq()
  223. {
  224.     DSequence* aSeq = new DSequence();
  225.     InsertLast( aSeq);
  226.     //aSeq->fIndex= this->GetSize();
  227. }
  228.  
  229.  
  230. void DSeqList::ClearSelections()
  231. {
  232.     long i, nseq= GetSize();
  233.     for (i= 0; i<nseq; i++)  
  234.         SeqAt(i)->ClearSelection();
  235. }
  236.  
  237. short DSeqList::ZeroOrigin()
  238. {
  239.     DSequence* aSeq;
  240.     short shiftall= 0;
  241.     long i, nseq= GetSize();
  242.     
  243.     for (i= 0; i<nseq; i++) {
  244.         aSeq= SeqAt(i);
  245.         if (aSeq->Origin()<0) shiftall= Min(shiftall, aSeq->Origin());
  246.         }
  247.         
  248.     for (i= 0; i<nseq; i++) {
  249.         aSeq= SeqAt(i);
  250.         aSeq->SetOrigin(aSeq->Origin() - shiftall); // -= shiftall;
  251.         if (aSeq->Origin() > 0) {
  252.             aSeq->InsertSpacers( 0, aSeq->Origin(), DSequence::indelEdge);
  253.             aSeq->SetOrigin(0);
  254.             }
  255.         }
  256.     return shiftall;
  257. }
  258.  
  259.  
  260.  
  261. short dslCompareByName( void* item1, void* item2)
  262. {
  263.     char*    name1= (char*) ((DSequence*)item1)->Name();
  264.     char*    name2= (char*) ((DSequence*)item2)->Name();
  265.     return StrICmp( name1, name2);
  266. }
  267.  
  268. short dslCompareBySize( void* item1, void* item2)
  269. {
  270.         // inverted sort, largest 1st
  271.     long name1= ((DSequence*)item1)->LengthF();
  272.     long name2= ((DSequence*)item2)->LengthF();
  273.     if (name1 > name2) return -1;
  274.     else if (name1 < name2) return 1;
  275.     else  return 0;
  276. }
  277.  
  278. short dslCompareByKind( void* item1, void* item2)
  279. {
  280.     short name1= ((DSequence*)item1)->Kind();
  281.     short name2= ((DSequence*)item2)->Kind();
  282.     if (name1 > name2) return 1;
  283.     else if (name1 < name2) return -1;
  284.     else  return 0;
  285. }
  286.  
  287. short dslCompareByDate( void* item1, void* item2)
  288. {
  289.         // inverted sort, latest 1st
  290.     unsigned long name1= ((DSequence*)item1)->ModTime();
  291.     unsigned long name2= ((DSequence*)item2)->ModTime();
  292.     if (name1 > name2) return -1;
  293.     else if (name1 < name2) return 1;
  294.     else  return 0;
  295. }
  296.  
  297.  
  298. Boolean DSeqList::SortList(Sorts sortorder)
  299. {
  300.     if (sortorder == fSortOrder)
  301.         return false;
  302.     else {
  303.         fSortOrder= sortorder;
  304.         switch( sortorder) {
  305.             case kSortByKind :
  306.                 this->SortBy( dslCompareByKind);
  307.                 break;
  308.             case kSortBySize : 
  309.                 this->SortBy( dslCompareBySize);
  310.                 break;
  311.             case kSortByDate :
  312.                 this->SortBy( dslCompareByDate);
  313.                 break;
  314.             case kSortByName : 
  315.                 this->SortBy( dslCompareByName);
  316.                 break;
  317.             case kSortByItem: 
  318.             default:
  319.                 this->Sort();
  320.                 break;
  321.             }
  322.         return true;
  323.         }
  324. }
  325.  
  326. DSequence* DSeqList::FindNamedSeq(char* name, Boolean respectCase)
  327. {    
  328.     if (name && *name) {
  329.         DSequence* aSeq;
  330.         long i, nseq= GetSize();
  331.         if (respectCase) for (i= 0; i<nseq; i++) {
  332.             aSeq= SeqAt(i);
  333.           if (Nlm_StringCmp(aSeq->Name(), name)==0) return aSeq;
  334.             }
  335.         else for (i= 0; i<nseq; i++) {
  336.             aSeq= SeqAt(i);
  337.             if (Nlm_StringICmp(aSeq->Name(), name)==0) return aSeq;
  338.             }
  339.         }
  340.     return NULL;    
  341. }
  342.  
  343. DSequence* DSeqList::Consensus(void)
  344. {    
  345.     long i, nseq= GetSize();
  346.     for (i= 0; i<nseq; i++) {
  347.         DSequence* aSeq= SeqAt(i);
  348.         if (aSeq->IsConsensus()) return aSeq;
  349.         }
  350.     return NULL;    
  351. }
  352.  
  353.  
  354. short DSeqList::ConsensusRow()
  355. {
  356.     DSequence* cons= this->Consensus();        
  357.     if (cons) 
  358.         return this->GetIdentityItemNo( cons);
  359.     else
  360.         return -1;
  361. }
  362.  
  363.  
  364.  
  365.     
  366. void DSeqList::MakeConsensus()
  367. {        
  368.     DSequence* cons= this->Consensus();        
  369.     if (!cons) {
  370.         cons = new DSequence(); 
  371.         cons->SetName(DSequence::kConsensus);
  372.         InsertLast( cons); 
  373.         }
  374.         
  375.     if (cons) {
  376.         //short arow= this->GetIdentityItemNo( cons);
  377.         char* hCon= cons->Bases();   
  378.         long conlen= cons->LengthF(); //StrLen( hCon);
  379.         long count= 0;
  380.         long i, iseq, nseq= GetSize();
  381.         for (iseq= 0; iseq<nseq; iseq++) {
  382.             DSequence* aSeq= this->SeqAt(iseq);
  383.             if (!aSeq->IsConsensus() 
  384.             //&& aSeq->Kind() != DSequence::kOtherSeq
  385.             && (aSeq->Kind() == DSequence::kDNA 
  386.                || aSeq->Kind() == DSequence::kRNA
  387.                || aSeq->Kind() == DSequence::kNucleic)
  388.              ) {
  389.                 count++;
  390.                 char* hSeq= aSeq->Bases();
  391.                 long  len = aSeq->LengthF();
  392.                 if (len > conlen) {
  393.                     conlen= len;
  394.                     hCon= (char*) MemMore( hCon, conlen);
  395.                     cons->UpdateBases( hCon, conlen);
  396.                     }
  397.                 if (count==1) {
  398.                     for (i= 0; i<len; i++) hCon[i]= hSeq[i];
  399.                     }            
  400.                 else {
  401.                     for (i= 0; i<len; i++) 
  402. // !! THis is only good for Nucleic bases, not for Aminos !!
  403. // change to table of matches !!
  404.                         hCon[i]= DSequence::NucleicConsensus(
  405.                             DSequence::NucleicBits(hCon[i]), 
  406.                             DSequence::NucleicBits(hSeq[i]));  
  407.                     }
  408.                 }
  409.             }
  410.         }
  411. }
  412.  
  413.  
  414.  
  415.  
  416.  
  417. char* DSeqList::Distances( short correct, Boolean doSimilarity)
  418. {  
  419.     struct NucCounts {
  420.         float apct, cpct, gpct, tpct;
  421.         };
  422.     long iseq, nseq= GetSize();
  423.     NucCounts* nuccounts = NULL;
  424.     char    linebuf[256], * buf = NULL;
  425.     char    name0[128];
  426.     
  427.     if (!nseq) return NULL;
  428.     if (correct == DColsen) nuccounts= new NucCounts[ nseq];
  429.     buf= (char*) MemNew(1);
  430.     *buf= 0;
  431.     *name0= 0;
  432.     
  433.     if (correct == DColsen) {
  434.         for (iseq=0; iseq<nseq; iseq++) {
  435.             long    na, nc, ng, nt, ntotal;
  436.             DSequence* aSeq= this->SeqAt(iseq);
  437.             char* bases= aSeq->Bases();
  438.             long  ibase, nbases= aSeq->LengthF();
  439.             na= nc= ng= nt= ntotal=0;
  440.                     // create Olsen correction array
  441.             for (ibase=0; ibase<nbases; ibase++) {
  442.                 long baseval= DSequence::NucleicBits(bases[ibase]);
  443.                 if (baseval & DSequence::kMaskA) na++;
  444.                 if (baseval & DSequence::kMaskC) nc++;
  445.                 if (baseval & DSequence::kMaskG) ng++;
  446.                 if (baseval & DSequence::kMaskT) nt++;
  447.                 if (baseval) ntotal++;
  448.                 }
  449.             nuccounts[iseq].apct = (float) na / (float) ntotal;
  450.             nuccounts[iseq].cpct = (float) nc / (float) ntotal;
  451.             nuccounts[iseq].gpct = (float) ng / (float) ntotal;
  452.             nuccounts[iseq].tpct = (float) nt / (float) ntotal;
  453.             }
  454.         }
  455.     
  456.     for (long j=1; j<nseq; j++)
  457.         for (long k=0; k<j; k++) {
  458.             long jk;
  459.             long numer= 0, denom= 0;
  460.             DSequence* jseq= this->SeqAt( j);
  461.             DSequence* kseq= this->SeqAt( k);
  462.             long len= MIN( jseq->LengthF(), kseq->LengthF());
  463.             for (jk= 0; jk < len; jk++) {
  464.                 long jval= jseq->NucleicBits(jseq->Bases()[jk]);
  465.                 long kval= kseq->NucleicBits(kseq->Bases()[jk]);
  466.                 if ( jval & kval) numer++;
  467.                 if ( jval | kval) denom++;
  468.                 }
  469.                 
  470.             float distv, cor;
  471.             switch (correct) {
  472.               case DColsen: 
  473.                   cor = nuccounts[j].apct * nuccounts[k].apct
  474.                       + nuccounts[j].cpct * nuccounts[k].cpct
  475.                       + nuccounts[j].gpct * nuccounts[k].gpct
  476.                       + nuccounts[j].tpct * nuccounts[k].tpct;
  477.                     break;
  478.                 case DCjukes: cor = 0.25; break;
  479.                 case DCnone: 
  480.                 default:             cor = 0.0;  break;
  481.                 }
  482.         
  483.             if (correct == DCnone)
  484.                 distv= 1.0 - (float)numer / (float)denom;
  485.             else
  486.                 distv= -(1.0-cor) * log( 1.0 / (1.0 - cor)
  487.                     * ((float)numer / (float)denom - cor) );
  488.             
  489.             if (k == 0) {
  490.                 if (j == 1) {  
  491.                     sprintf( name0, "%2d=%-15s",k+1,kseq->Name());
  492.                     //StrExtendCat( &buf, linebuf);
  493.                     }
  494.                 sprintf( linebuf, "%2d:%-15s|",j+1,jseq->Name());
  495.                 StrExtendCat( &buf, linebuf);
  496.                 }
  497.                 
  498.             if (doSimilarity)
  499.                 sprintf(linebuf, "%6.1f",100 - distv*100.0);
  500.             else
  501.                 sprintf(linebuf, "%6.1f",distv*100.0);
  502.             StrExtendCat( &buf, linebuf);
  503.         
  504.             if (k == j-1) {
  505.                 StrExtendCat( &buf, " \n");
  506.                 }
  507.                 
  508.             }
  509.  
  510.     sprintf( linebuf, "   %-15s ","-------------");
  511.     StrExtendCat( &buf, linebuf);
  512.     for (iseq=1; iseq<nseq; iseq++) {
  513.         sprintf(linebuf, " %5s", "----");
  514.         StrExtendCat( &buf, linebuf);
  515.         }
  516.     sprintf( linebuf, "\n%-18s ",name0);
  517.     StrExtendCat( &buf, linebuf);
  518.     for (iseq=1; iseq<nseq; iseq++) {
  519.         sprintf(linebuf, "%5d ",iseq);
  520.         StrExtendCat( &buf, linebuf);
  521.         }
  522.         
  523.     char * cors, *skind;
  524.     switch (correct) {
  525.         case DCjukes: cors= "(Jukes correction)"; break;
  526.         case DColsen: cors= "(Olsen correction)"; break;
  527.         default:
  528.         case DCnone:  cors= ""; break;
  529.         }
  530.     if (doSimilarity) skind= "Similarity"; else skind="Distance";
  531.     sprintf( linebuf, "\n  %s matrix %s\n", skind, cors);
  532.     StrExtendCat( &buf, linebuf);
  533.             
  534.     delete[] nuccounts;
  535.     return buf;
  536. }
  537.  
  538.  
  539.  
  540.  
  541.  
  542. char* DSeqList::FindCommonBases(short minCommonPerCent, char*& firstCommon)
  543. {  
  544.     enum { kLetrange = 26, kLetlast = 255 };
  545.     char    * hSeq,  * hMaxbase;
  546.     char        ch;
  547.     long        len, maxlen, ibase;
  548.     short        spccount, nilcount;
  549.     short        maxlet, let, nseq, maxcnt;
  550.     short        letcount[kLetrange], letfirst[kLetrange];
  551.     float     rmaxc;
  552.  
  553.     DSequence* cons= this->Consensus();        
  554.     if (!cons) cons= (DSequence*) this->First(); //any seq will do...
  555.     
  556.     hMaxbase= NULL;
  557.     firstCommon= NULL;
  558.     if (cons) {        
  559.         maxlen= cons->LengthF();
  560.         hMaxbase= (char*) Nlm_MemGet(maxlen+1,true); 
  561.         firstCommon= (char*) Nlm_MemGet(maxlen+1,true); 
  562.         
  563.         for (ibase= 0; ibase<maxlen; ibase++) {
  564.             for (let= 0; let < kLetrange; let++) { letcount[let]= 0; letfirst[let]= kLetlast; }
  565.             nilcount= 0; 
  566.             spccount= 0;
  567.             nseq= 0;
  568.             
  569.             long iseq, maxseq= GetSize();
  570.             for (iseq=0; iseq<maxseq; iseq++) {
  571.                 DSequence* aSeq= this->SeqAt(iseq);
  572.                 if (!aSeq->IsConsensus() && aSeq->Kind() != DSequence::kOtherSeq) {
  573.                     hSeq= aSeq->Bases();
  574.                     len= aSeq->LengthF();
  575.                     if (ibase<len) ch= toupper( hSeq[ibase]);
  576.                     else ch= 0;
  577.                     if (ch >= 'A' && ch <= 'Z') {
  578.                         let = ch - 'A';
  579.                         letcount[let]++;
  580.                         if (letfirst[let] == kLetlast) letfirst[let]= Min(kLetlast, iseq);
  581.                         }
  582.                     else if (
  583.                          ch == DSequence::indelSoft 
  584.                         || ch == DSequence::indelHard 
  585.                         || ch == DSequence::indelEdge 
  586.                         || ch == ' ' || ch == '_' )  
  587.                         spccount++;
  588.                     else 
  589.                         nilcount++;
  590.                     nseq++;
  591.                     }
  592.                 }
  593.                 
  594.             for (let= 0, maxlet= 0, maxcnt= -1; let < kLetrange; let++) 
  595.                 if (letcount[let] > maxcnt) { 
  596.                     maxcnt= letcount[let]; 
  597.                     maxlet= let; 
  598.                     }
  599.                     
  600.             if (maxcnt>0) {
  601.                 firstCommon[ibase]= letfirst[maxlet];
  602.                 rmaxc= (100.0 * maxcnt) / nseq;
  603.                 if (rmaxc >= minCommonPerCent)   //mark maxlet's
  604.                     hMaxbase[ibase]= maxlet + 'A';
  605.                 else
  606.                     hMaxbase[ibase]= '!';
  607.                 }
  608.             else {
  609.                 firstCommon[ibase]= kLetlast;
  610.                 hMaxbase[ibase]= '!';
  611.                 }
  612.             }
  613.             
  614.         hMaxbase[maxlen]= 0;
  615.         firstCommon[maxlen]= 0;
  616.         }
  617.     
  618.     return hMaxbase;
  619. }
  620.  
  621.