gnGenomeSpec.cpp Source File

00001 
00002 // File:            gnGenomeSpec.cpp
00003 // Purpose:         implements gnMultiSpec for genomes
00004 // Description:     
00005 // Changes:        
00006 // Version:         libGenome 0.1.0 
00007 // Author:          Aaron Darling 
00008 // Last Edited:     April 15, 2001, 11:13:00pm 
00009 // Modified by:     
00010 // Copyright:       (c) Aaron Darling 
00011 // Licenses:        Proprietary 
00013 
00014 #include "gn/gnGenomeSpec.h"
00015 #include <string>
00016 
00017 gnGenomeSpec::gnGenomeSpec()
00018 {
00019         Clear();
00020 }
00021 
00022 gnGenomeSpec::gnGenomeSpec( const gnGenomeSpec& s )
00023 {
00024         m_sourceName = s.m_sourceName;
00025         m_name = s.m_name;
00026         m_reverseComplement = s.m_reverseComplement;
00027         m_circular = s.m_circular;
00028         //copy the header list.
00029         uint32 list_size =  s.m_headerList.size();
00030         m_headerList.reserve(list_size);
00031         for(uint32 i=0; i < list_size; i++)
00032                 m_headerList.push_back(s.m_headerList[i]->Clone());
00033         //copy the fragment list.
00034         list_size =  s.m_SpecList.size();
00035         m_SpecList.reserve(list_size);
00036         for(uint32 i=0; i < list_size; i++)
00037                 m_SpecList.push_back(s.m_SpecList[i]->Clone());
00038 }
00039 gnGenomeSpec::~gnGenomeSpec()
00040 {
00041         Clear();
00042 }
00043 // Clone
00044 void gnGenomeSpec::Clear()
00045 {
00046         uint32 list_size =  m_SpecList.size();
00047         for(uint32 i=0; i < list_size; i++)
00048                 delete m_SpecList[i];
00049         m_SpecList.clear();
00050         gnMultiSpec::Clear();
00051 }
00052 
00053 void gnGenomeSpec::SetReverseComplement( const boolean value )
00054 {
00055         if(value == m_reverseComplement)
00056                 return;
00057         //reverse the spec list entries
00058         vector<gnFragmentSpec*> tmp_spec_list;
00059         for(uint32 i=0; i < GetSpecListLength(); i++){
00060                 //transmit rev_comp down to each fragment spec
00061                 GetSpec(i)->SetReverseComplement(!GetSpec(i)->IsReverseComplement());
00062                 tmp_spec_list.insert(tmp_spec_list.begin(), GetSpec(i));
00063         }
00064         m_SpecList = tmp_spec_list;
00065         m_reverseComplement = value;
00066 }
00067 
00068 gnGenomeSpec* gnGenomeSpec::CloneRange( const gnSeqI startI, const gnSeqI len ) const{
00069         if(len == 0)
00070                 return new gnGenomeSpec();
00071 
00072         //find the valid range of specs to copy
00073         uint32 firstSpec = GetSpecIndexByBase(startI);
00074         gnSeqI total_copylen = len;
00075         uint32 endSpec;
00076         if(len != GNSEQI_END){
00077                 endSpec = GetSpecIndexByBase(startI + len - 1);
00078         }else{
00079                 endSpec = GetSpecListLength() - 1;
00080                 total_copylen = GetLength() - startI;
00081         }
00082 
00083         //find their starting and ending bases
00084         gnSeqI firstBase = startI - GetSpecStartBase(firstSpec);
00085         gnSeqI firstSpecLen = GetSpec(firstSpec)->GetLength();
00086         boolean spans_specs = true;
00087         gnSeqI firstCopyLen = firstSpecLen - firstBase;
00088         if(firstCopyLen >= total_copylen){
00089                 spans_specs = false;
00090                 firstCopyLen = total_copylen;
00091         }
00092         
00093         gnGenomeSpec* destSpec = new gnGenomeSpec();
00094         gnFragmentSpec* newSpec = m_SpecList[firstSpec]->CloneRange(firstBase, firstCopyLen);
00095         destSpec->AddSpec( newSpec );
00096 
00097         gnSeqI cur_copylen = firstCopyLen;
00098         //add all the completely covered specs in the middle
00099         for(uint32 specI = firstSpec + 2; specI <= endSpec; specI++){
00100                 destSpec->AddSpec(GetSpec(specI-1)->Clone());
00101                 cur_copylen += GetSpec(specI-1)->GetLength();
00102         }
00103         //add the last spec if necessary
00104         if(spans_specs){
00105                 newSpec = m_SpecList[endSpec]->CloneRange( 0, total_copylen - cur_copylen);
00106                 destSpec->AddSpec(newSpec);
00107         }
00108         return destSpec;
00109 }
00110 
00111 
00112 //IMPLEMENT ME: What to do with headers?
00113 void gnGenomeSpec::MergeFragments(const uint32 startF, const uint32 endF){
00114         if(startF > m_SpecList.size() || endF > m_SpecList.size())
00115                 Throw_gnEx(FragmentIndexOutOfBounds());
00116         if(startF > endF)
00117                 Throw_gnEx(FragmentIndexOutOfBounds());
00118         
00119         for(uint32 i = startF + 1; i < endF; i++){
00120                 gnFragmentSpec* delspec = m_SpecList[startF + 1];
00121                 m_SpecList.erase(m_SpecList.begin() + startF + 1);
00122 
00123                 for(uint32 j = 0; j < delspec->GetSpecListLength(); j++)
00124                         m_SpecList[startF]->AddSpec(delspec->GetSpec(j));
00125                 delete delspec;
00126         }
00127 }
00128 
00129 uint32 gnGenomeSpec::AddFeature( gnBaseFeature* feat ){
00130         uint32 count = 0;
00131         uint32 len = 0;
00132         uint32 featureI = 0;
00133         uint32 specListLen = GetSpecListLength();
00134 
00135         for(uint32 specI = 0; specI < specListLen; specI++){
00136                 len = GetSpec(specI)->GetLength();
00137                 gnLocation lt(count, count+len);
00138                 if(feat->IsContainedBy(lt)){
00139                         return featureI + GetSpec(specI)->AddFeature(feat);
00140                 }
00141                 count += len;
00142                 featureI += GetSpec(specI)->GetFeatureListLength();
00143         }
00144         //if we get this far then the feature has invalid coordinates
00145         Throw_gnEx(SeqIndexOutOfBounds());
00146 }
00147 
00148 uint32 gnGenomeSpec::GetFeatureListLength() const
00149 {
00150         uint32 len = 0;
00151         for(uint32 i=0; i < GetSpecListLength(); i++)
00152                 len += GetSpec(i)->GetFeatureListLength();
00153         return len;
00154 }
00155 
00156 gnBaseFeature* gnGenomeSpec::GetFeature(const uint32 i ) const
00157 {
00158         uint32 count = 0;
00159         uint32 len = 0;
00160         for(uint32 specI=0; specI < GetSpecListLength(); specI++){
00161                 len = GetSpec(specI)->GetFeatureListLength();
00162                 if(count <= i && i < count + len){
00163                         gnBaseFeature* feat = GetSpec(specI)->GetFeature(i - count);
00164                         feat->MovePositive(GetSpecStartBase(specI));
00165                         return feat;
00166                 }
00167                 count += len;
00168         }
00169         Throw_gnEx(FeatureIndexOutOfBounds());
00170 }
00171 
00172 void gnGenomeSpec::GetContainedFeatures(const gnLocation& lt, vector<gnBaseFeature*>& feature_vector, vector<uint32>& index_vector) const{
00173         uint32 ss_size = GetSpecListLength();
00174         uint32 fl_size = 0;
00175         gnSeqI start_base = 0;
00176         for(uint32 i=0; i < ss_size; i++){
00177                 gnLocation sub_lt = lt;
00178                 gnSeqI sub_len = GetSpec(i)->GetLength();
00179                 sub_lt.MoveNegative(start_base);
00180                 sub_lt.CropEnd(sub_len);
00181                 GetSpec(i)->GetContainedFeatures(sub_lt, feature_vector, index_vector);
00182                 uint32 fvs = feature_vector.size();
00183                 for(uint32 j = 0; j < fvs; j++){
00184                         feature_vector[j]->MovePositive(start_base);
00185                         index_vector[j]+= fl_size;
00186                 }
00187                 if(fvs > 0)
00188                         return;
00189                 start_base += sub_len;
00190                 fl_size += GetSpec(i)->GetFeatureListLength();
00191         }
00192 }
00193 
00194 void gnGenomeSpec::GetIntersectingFeatures(const gnLocation& lt, vector<gnBaseFeature*>& feature_vector, vector<uint32>& index_vector) const{
00195         uint32 ss_size = GetSpecListLength();
00196         uint32 fl_size = 0;
00197         gnSeqI start_base = 0;
00198         for(uint32 i=0; i < ss_size; i++){
00199                 gnLocation sub_lt = lt;
00200                 gnSeqI sub_len = GetSpec(i)->GetLength();
00201                 sub_lt.MoveNegative(start_base);
00202                 sub_lt.CropEnd(sub_len);
00203                 GetSpec(i)->GetIntersectingFeatures(sub_lt, feature_vector, index_vector);
00204                 uint32 fvs = feature_vector.size();
00205                 for(uint32 j = 0; j < fvs; j++){
00206                         feature_vector[j]->MovePositive(start_base);
00207                         index_vector[j]+= fl_size;
00208                 }
00209                 if(fvs > 0)
00210                         return;
00211                 start_base += sub_len;
00212                 fl_size += GetSpec(i)->GetFeatureListLength();
00213         }
00214 }
00215 void gnGenomeSpec::GetBrokenFeatures(const gnLocation& lt, vector<gnBaseFeature*>& feature_vector) const{
00216         uint32 ss_size = GetSpecListLength();
00217         uint32 fl_size = 0;
00218         gnSeqI start_base = 0;
00219         for(uint32 i=0; i < ss_size; i++){
00220                 gnLocation sub_lt = lt;
00221                 gnSeqI sub_len = GetSpec(i)->GetLength();
00222                 sub_lt.MoveNegative(start_base);
00223                 sub_lt.CropEnd(sub_len);
00224                 GetSpec(i)->GetBrokenFeatures(sub_lt, feature_vector);
00225                 uint32 fvs = feature_vector.size();
00226                 for(uint32 j = 0; j < fvs; j++)
00227                         feature_vector[j]->MovePositive(start_base);
00228                 if(fvs > 0)
00229                         return;
00230                 start_base += sub_len;
00231                 fl_size += GetSpec(i)->GetFeatureListLength();
00232         }
00233 }
00234 
00235 void gnGenomeSpec::RemoveFeature( const uint32 i ){
00236         uint32 count = 0;
00237         uint32 len = 0;
00238         uint32 specI=0;
00239         for(; specI < GetSpecListLength(); specI++){
00240                 len = GetSpec(specI)->GetFeatureListLength();
00241                 if(count <= i && i < count + len){
00242                         GetSpec(specI)->RemoveFeature( i - count);
00243                 }
00244                 count += len;
00245         }
00246         Throw_gnEx(FeatureIndexOutOfBounds());
00247 }