Enter your search terms GNU-Darwin Web Submit search form

testTests.cpp Go to the documentation of this file. #include "gn/gnSequence.h" #include "gn/gnGBKSource.h" #include "gn/gnFASSource.h" #include <algorithm> static boolean LocationLessthan(const gnLocation& a, const gnLocation& b); static boolean LocationLessthan(const gnLocation& a, const gnLocation& b){ return a.GetStart() < b.GetStart(); } static boolean LocationEndLessthan(const gnLocation& a, const gnLocation& b); static boolean LocationEndLessthan(const gnLocation& a, const gnLocation& b){ return a.GetEnd() < b.GetStart(); } static boolean LocationSizeLessthan(const gnLocation& a, const gnLocation& b); static boolean LocationSizeLessthan(const gnLocation& a, const gnLocation& b){ return a.GetEnd() - a.GetStart() < b.GetEnd() - b.GetStart(); } // create a location list which contains every intervening region // in l_list void ComplementLocation(vector<gnLocation>& l_list){ if(l_list.size() < 2) return; //need 2 or more locations to get the in between regions vector<gnLocation> comp_list; sort(l_list.begin(), l_list.end(), &LocationLessthan); for(uint32 locationI = 1; locationI < l_list.size(); locationI++){ //if they don't intersect there is an intervening region if( !l_list[locationI].Intersects( l_list[locationI-1] ) ){ //create the region between locationI and locationI - 1 gnLocation between_loc; between_loc.SetStart(l_list[locationI-1].GetEnd()); between_loc.SetEnd(l_list[locationI].GetStart()); //add it to the complement list comp_list.push_back(between_loc); } } l_list = comp_list; } void DeleteRepeats(list<gnLocation>& loc_list){ loc_list.sort(&LocationLessthan); list<gnLocation>::iterator loc_iter = loc_list.begin(); list<gnLocation>::iterator cur_loc; for(; loc_iter != loc_list.end(); loc_iter++){ cur_loc = loc_iter; cur_loc++; while(cur_loc != loc_list.end() && !LocationEndLessthan(*loc_iter, *cur_loc)){ if(loc_iter->Intersects(*cur_loc)){ *loc_iter = loc_iter->GetUnion(*cur_loc); list<gnLocation>::iterator to_del = cur_loc; cur_loc--; loc_list.erase(to_del); } cur_loc++; } } } void WriteData(gnSequence& file_sequence, list<gnLocation>& loc_list, string& filename){ gnSequence loc_seq; cout << "Gathering sequence data to write...\n"; list<gnLocation>::iterator loc_iter = loc_list.begin(); uint32 loc_count = loc_list.size(); uint32 notice_interval = loc_count / 10; uint32 cur_locI = 0; for(;loc_iter != loc_list.end(); loc_iter++){ loc_seq += file_sequence.ToString( loc_iter->GetEnd() - loc_iter->GetStart(), loc_iter->GetStart()); if(cur_locI % notice_interval == 0){ cout << (cur_locI / loc_count) * 10 << "%.."; } cur_locI++; } cout << "\nWriting sequence data...\n"; gnFASSource::Write(loc_seq, filename); } void WriteList(list<gnLocation>& loc_list, string& filename){ ofstream list_out(filename.c_str()); list<gnLocation>::iterator loc_iter = loc_list.begin(); for(;loc_iter != loc_list.end(); loc_iter++){ list_out << loc_iter->GetStart() << '\t' << loc_iter->GetEnd() << '\n'; } list_out.close(); } void print_usage(char* pname){ cout << "Usage: " << pname << " <genbank file> <exon_list> <intron_list> <exon_seq> <intron_seq>\n"; } int main( int argc, char* argv[] ) { boolean run_interactive = false; string seq_filename; string exon_list_filename; string intron_list_filename; string exon_seq_filename; string intron_seq_filename; // define a gnSequence to store the sequence gnSequence file_sequence; if(!run_interactive){ // check for correct calling semantics if(argc != 6){ print_usage(argv[0]); return -1; } seq_filename = argv[1]; exon_list_filename = argv[2]; intron_list_filename = argv[3]; exon_seq_filename = argv[4]; intron_seq_filename = argv[5]; }else{ // Get the name of the sequence to load cout << "Enter the name of the sequence file to load:\n" ; cin >> seq_filename; } // Load the sequence and tell the user if it loaded successfully if(file_sequence.LoadSource(seq_filename)) { cout << seq_filename << " loaded successfully. " << file_sequence.length() << " base pairs.\n"; }else{ cout << "Error loading file.\n"; return -1; } uint32 feature_count = file_sequence.getFeatureListLength(); uint32 cds_count = 0; //define a sequence for each type of data gnSequence exon_seq; gnSequence intron_seq; list<gnLocation> exon_list; list<gnLocation> intron_list; cout << "There are " << feature_count << " known features.\n"; for(uint32 featureI = 0; featureI < feature_count; featureI++){ gnBaseFeature* cur_feat = file_sequence.getFeature(featureI); if(cur_feat == NULL){ cout << "cur_feat is NULL, trace me!\n"; cur_feat = file_sequence.getFeature(featureI); continue; } if(cur_feat->GetName() == "CDS"){ cds_count++; uint32 loc_count = cur_feat->GetLocationListLength(); vector<gnLocation> cur_exon_list, cur_intron_list; for(uint32 locationI = 0; locationI < loc_count; locationI++){ cur_exon_list.push_back(cur_feat->GetLocation(locationI)); exon_list.push_back(cur_exon_list[locationI]); } if(loc_count >= 2){ cur_intron_list = cur_exon_list; ComplementLocation(cur_intron_list); for(uint32 locationI = 0; locationI < cur_intron_list.size(); locationI++) intron_list.push_back(cur_intron_list[locationI]); } /* Code for exon size tracking */ /* sort(cur_exon_list.begin(), cur_exon_list.end(), &LocationSizeLessthan); gnSeqI cur_size = cur_exon_list[loc_count - 1].GetEnd() - cur_exon_list[loc_count - 1].GetStart(); size_out << cur_size << '\n'; if(cur_size > 5000){ cout << "psycho bob is at it again\n"; cout << "Found " << cur_size << " length exon in feature: " << cur_feat->GetName() << "\n"; cout << "Feature Index: " << featureI << "\n"; for(uint32 i=0; i < cur_feat->GetQualifierListLength(); i++){ cout << cur_feat->GetQualifierName(i) << "\t" << cur_feat->GetQualifierValue(i) << "\n"; } cout << "\n"; } */ delete cur_feat; } } cout << "There are " << cds_count << " cds features.\n"; //now sort the lists and delete repeats and fix overlaps uint32 exon_count = exon_list.size(); cout << "Deleting exon repeats... "; DeleteRepeats(exon_list); cout << "Eliminated " << exon_count - exon_list.size() << " repeats for a total of "; cout << exon_list.size() << " repeats\n"; uint32 intron_count = intron_list.size(); cout << "Deleting intron repeats... "; DeleteRepeats(intron_list); cout << "Eliminated " << intron_count - intron_list.size() << " repeats for a total of "; cout << intron_list.size() << " repeats\n"; WriteList(exon_list, exon_list_filename); WriteList(intron_list, intron_list_filename); WriteData(file_sequence, exon_list, exon_seq_filename); WriteData(file_sequence, intron_list, intron_seq_filename); if(run_interactive) cin >> seq_filename; return 0; };