RetroSearch Browse

192 const string

) :

194

x_CalcWhat(file_or_accn, level,

)),

195

m_file_or_accn_basename(

196

x_CalcFileBaseNameOrAccn(file_or_accn)),

203

~SOutMessage() throw() { }

205 const string

m_file_or_accn_basename;

206 const string

m_level;

214 static string

x_CalcFileBaseNameOrAccn(

215 const string

& file_or_accn);

217 static string

x_CalcWhat(

218 const string

& file_or_accn,

219 const string

& level,

220 const string

221 const string

);

224 const string

SOutMessage::kErrorStr(

"ERROR"

);

225 const string

SOutMessage::kFatalStr(

"FATAL"

);

227 void

SOutMessage::WriteAsXML(

& out_strm)

const 230 xml::node

& expn_root_node = expn_doc.get_root_node();

234 if

( ! m_file_or_accn_basename.empty() ) {

235

expn_root_attribs.

(

236 "file_or_accn"

, m_file_or_accn_basename.c_str());

238

expn_root_attribs.

(

"severity"

, m_level.c_str());

239

expn_root_attribs.

(

"code"

, m_code.c_str());

246 void

SOutMessage::WriteAsText(

& out_strm)

const 248

out_strm << what() << endl;

251 string

SOutMessage::x_CalcFileBaseNameOrAccn(

252 const string

& file_or_accn)

254 CFile

maybe_file(file_or_accn);

255 if

( maybe_file.IsFile() ) {

257 return

maybe_file.GetName();

264 string

SOutMessage::x_CalcWhat(

265 const string

& file_or_accn,

266 const string

& level,

267 const string

268 const string

)

273 const string

file_or_accn_basename =

274

x_CalcFileBaseNameOrAccn(file_or_accn);

276

answer_str_strm << level <<

": ["

"] "

;

277 if

( ! file_or_accn_basename.empty() ) {

278

answer_str_strm << file_or_accn_basename <<

": "

;

280

answer_str_strm <<

;

297 virtual void Init

(

void

);

298 virtual int Run

(

void

);

340 const

SOutMessage &out_message,

& out_strm)

const

;

357

m_eOutFormat(eOutFormat_ASCIITable)

373

arg_desc->SetUsageContext(

GetArguments

().GetProgramBasename(),

374 "Gap analysis program"

false

);

378

arg_desc->SetCurrentGroup(

"BASIC"

);

382 "The files or accessions to do gap analysis on. " 383 "ASN.1, XML, and FASTA are some of the supported file formats."

386

vector<string> gap_types;

387 ITERATE

(TAddGapTypeMap, add_gap_type_it, sc_addgaptypename) {

388

gap_types.push_back(add_gap_type_it->first);

391

arg_desc->AddOptionalKey(

394 "A comma-separated list of types of gaps we look at. " 395 "If none are specified, all types will be shown. Possibilities: "

396 NStr::Join

(gap_types,

", "

) +

" but note that 'all' means adding " 397 "one column that combines all types. Also, there is " 398 "the gap-type 'both' as a shortcut for showing a column " 399 "for every possible gap type " 400 "(except 'all', unless you also add that to -gap-types)"

403

arg_desc->AddDefaultKey(

404 "out-format"

"Format"

405 "Specifies how the results should be printed."

409 "ascii-table"

"xml"

));

411

arg_desc->AddDefaultKey(

"sort-on"

"HowToSortResults"

412 "Specify the order of the summary. length sorts on the gap length. " 413 "num_seqs sorts on the number of sequences " 414 "each gap length appears on. " 415 "num_gaps sorts on the total number of times a gap " 416 "of each gap length appears. "

420 "length"

"num_seqs"

"num_gaps"

));

422

arg_desc->AddFlag(

"rev-sort"

423 "Set this to reverse the sorting order."

);

425

arg_desc->AddFlag(

"show-seqs-for-gap-lengths"

426 "This will show the sequences that each gap size has. It " 427 "is not affected by the sorting options."

);

429

arg_desc->AddDefaultKey(

"mol"

"MolTypesToLookAt"

430 "Specify whether you just want to look at sequences which are " 431 "nucleic acids (na), amino acids (aa), or any."

435 "na"

"aa"

"any"

));

437

arg_desc->AddDefaultKey(

"assume-mol"

"AssumedMolType"

438 "If unable to determine mol of sequence from ASN.1, " 439 "FASTA mods, etc. this specifies what mol we assume it is."

445

arg_desc->SetCurrentGroup(

"RESOLVING SEQUENCES"

);

448 "no-gbload"

"If set, avoids using the GenBank data loader"

);

450

arg_desc->AddDefaultKey(

451 "max-resolve-count"

"MAX_COUNT"

452 "How deep to resolve when following far references. Note that being " 453 "unable to resolve a far reference whether due to reaching max " 454 "resolve count or any other reason will result in an error since " 455 "unable to calculate 'gaps' if cannot get sequence. Zero disables " 456 "all far reference resolution."

460

arg_desc->SetConstraint(

463

arg_desc->SetCurrentGroup(

"HISTOGRAM"

);

465

arg_desc->AddFlag(

"show-hist"

466 "Set this flag to see a histogram of gap data"

);

468

arg_desc->AddDefaultKey(

"hist-bins"

"NumBins"

469 "Set the number of histogram bins to aim for (not guaranteed " 470 "to be that exact number). Default is 0, which means to " 471 "automatically pick a reasonable number of bins"

474

arg_desc->SetConstraint(

"hist-bins"

476

arg_desc->SetDependency(

479

arg_desc->AddDefaultKey(

"hist-algo"

"HistogramAlgorithm"

480 "Set this if you want the histogram binner to try to " 481 "use a different binning algorithm. The default should " 482 "be fine for most people."

486 "cluster"

"even_bins"

));

487

arg_desc->SetDependency(

514 const string

& sOutFormat = args[

"out-format"

].AsString();

515 if

(

"ascii-table"

== sOutFormat ) {

517

}

else if

(

"xml"

== sOutFormat ) {

521

pResultBeginEndStr.

reset

(

522 new

CBeginEndStrToCoutRAII(

"<result>"

"</result>"

));

533

}

catch

(

const

std::exception & ex ) {

534

p_out_message.

reset

(

537 "INTERNAL_ERROR"

, ex.what()));

539

p_out_message.

reset

(

542 "INTERNAL_ERROR"

"(---UNKNOWN INTERNAL ERROR TYPE---)"

));

570 const string

& sMol = args[

"mol"

].AsString();

571 if

( sMol ==

"na"

) {

573

}

else if

( sMol ==

"aa"

) {

575

}

else if

( sMol ==

"any"

) {

584 const string

& sAssumeMol = args[

"assume-mol"

].AsString();

585 if

( sAssumeMol ==

"na"

) {

587

}

else if

( sAssumeMol ==

"aa"

) {

595 if

( args[

"gap-types"

HasValue

() ) {

596 string

gap_type_arg = args[

"gap-types"

].AsString();

602

vector< string > raw_gap_types;

603 NStr::Split

(args[

"gap-types"

].AsString(),

","

, raw_gap_types, 0);

604 if

( raw_gap_types.empty() ) {

607 "NO_GAP_TYPES_GIVEN"

608 "-gap-types must be given at least one gap"

);

612 ITERATE

(vector< string >, raw_gap_types_it, raw_gap_types)

614 const string

& raw_gap_type_str = *raw_gap_types_it;

620

vector<string> gap_types_from_this_arg;

621 if

( raw_gap_type_str ==

"both"

) {

624 const size_t

old_num_gap_types =

625

gap_types_from_this_arg.size();

626 ITERATE

(TAddGapTypeMap, add_gap_type_it, sc_addgaptypename) {

627 const string

& gap_type_str = add_gap_type_it->first;

628 if

( gap_type_str !=

"all"

) {

629

gap_types_from_this_arg.push_back(gap_type_str);

632 const size_t

new_num_gap_types =

633

gap_types_from_this_arg.size();

635 if

( (old_num_gap_types + 2) != new_num_gap_types ) {

643

gap_types_from_this_arg.push_back(raw_gap_type_str);

646 ITERATE

(vector<string>, gap_type_str_it, gap_types_from_this_arg)

648 const string

& gap_type_str = *gap_type_str_it;

649

TAddGapTypeMap::const_iterator find_it =

650

sc_addgaptypename.find(gap_type_str.c_str());

651 if

( find_it == sc_addgaptypename.end() ) {

654

SOutMessage::kFatalStr,

656 "This gap type is not recognized: '"

657

gap_type_str +

"'"

);

660 const

SGapRelatedInfo & gap_related_flags =

672 m_eSort

= GA::eSortGapLength_Length;

673 const string

sSortOn = args[

"sort-on"

].AsString();

674 if

(

"length"

== sSortOn ) {

675 m_eSort

= GA::eSortGapLength_Length;

676

}

else if

(

"num_seqs"

== sSortOn ) {

677 m_eSort

= GA::eSortGapLength_NumSeqs;

678

}

else if

(

"num_gaps"

== sSortOn ) {

679 m_eSort

= GA::eSortGapLength_NumGaps;

687

GA::eSortDir_Descending :

688

GA::eSortDir_Ascending );

692 for

(

size_t

ii = 1; ii <= args.

GetNExtra

(); ++ii ) {

693 const string

sFileOrAccn = args[ii].AsString();

697 catch

(

const

ncbi::objects::CObjMgrException& ex) {

698 if

(ex.GetErrCode() == ncbi::objects::CObjMgrException::eAddDataError

699

&& ex.GetMsg().find(

"duplicate Bioseq id"

) == 0) {

700

SOutMessage out_message(

701

sFileOrAccn, SOutMessage::kErrorStr,

702

ex.GetErrCodeString(), ex.GetMsg());

707 catch

(

const

SOutMessage & out_message) {

714 catch

(

const

ncbi::CException& ex) {

715

SOutMessage out_message(

716

sFileOrAccn, SOutMessage::kFatalStr,

717

ex.GetErrCodeString(), ex.GetMsg());

728

SOutMessage out_message(

729

sFileOrAccn, SOutMessage::kFatalStr,

730 "INTERNAL_ERROR"

"Unexpected exception"

);

740 if

( args[

"show-seqs-for-gap-lengths"

] ) {

744 if

( args[

"show-hist"

] ) {

745 const Uint8

num_bins = args[

"hist-bins"

].AsInt8();

746 const string

& sHistAlgo = args[

"hist-algo"

].AsString();

749 if

(

"cluster"

== sHistAlgo ) {

751

}

else if

(

"even_bins"

== sHistAlgo ) {

778 if

( !

GetArgs

()[

"no-gbload"

] ) {

791 string

answer = gap_name;

818 ITERATE

( GA::TVectorGapLengthSummary,

819

summary_unit_it, *pGapLenSummary )

821

all_gap_lengths_list->push_back((*summary_unit_it)->gap_length);

827

== all_gap_lengths_list->end());

830 return

all_gap_lengths_list;

841 if

(

CDirEntry

(sFileOrAccn).Exists() ) {

843 if

( !

CDirEntry

(sFileOrAccn).IsFile() ) {

846

SOutMessage::kErrorStr,

848 FORMAT

(

"This exists but is not a plain file: " 883

pSeqEntry = fasta_reader.

ReadSet

();

890

SOutMessage::kErrorStr,

891 "UNSUPPORTED_FORMAT"

892 FORMAT

(

"This format is not yet supported: " 906 if

( ! pSeqSubmit->IsEntrys() ||

907

pSeqSubmit->GetData().GetEntrys().size() != 1 )

911

SOutMessage::kErrorStr,

912 "SEQ_SUBMIT_MULTIPLE_SEQ_ENTRIES"

914 "Only Seq-submits with exactly " 916 "inside are supported."

));

918

pSeqEntry = *pSeqSubmit->SetData().SetEntrys().begin();

931

pSeqEntry = pNewSeqEntry;

945

pSeqEntry->

SetSeq

( *pBioseq );

953

SOutMessage::kErrorStr,

954 "INVALID_FORMAT_OR_BAD_OBJ_TYPE"

955 FORMAT

(

"Invalid ASN.1 or unsupported object type"

));

969

SOutMessage::kErrorStr,

978

SOutMessage::kErrorStr,

979 "ACCESSION_NOT_FOUND"

980 FORMAT

(

"Accession could not be found"

));

1010 xml::node

gap_len_infos_node(

"gap_len_infos"

);

1014 typedef

pair<GA::TGapLength, GA::EGapType>

1017

TGapLenTypeToSummaryMap;

1018

TGapLenTypeToSummaryMap gap_len_type_to_summary_map;

1022 const

GA::EGapType eGapType = gap_type_name_it->first;

1030 ITERATE

(GA::TVectorGapLengthSummary, gap_summary_it,

1031

*p_gap_len_summary )

1035

TGapLenTypeKey gap_len_type(p_one_summary->gap_length, eGapType);

1036

pair<TGapLenTypeToSummaryMap::iterator, bool> insert_ret =

1037

gap_len_type_to_summary_map.insert(

1042 if

( ! insert_ret.second ) {

1052 TGapLengthVec

& all_gap_lengths_list = *p_all_gap_lengths_list;

1058 const

GA::TGapLength gap_len = *gap_length_it;

1060 xml::node

one_gap_len_info(

"one_gap_len_info"

);

1064

one_gap_len_attributes.

(

1070 const

GA::EGapType eGapType = gap_type_name_it->first;

1071 const CTempString

pchGapName = gap_type_name_it->second;

1080

TGapLenTypeToSummaryMap::const_iterator find_it =

1081

gap_len_type_to_summary_map.find(

1082

make_pair(gap_len, eGapType));

1083 if

( find_it != gap_len_type_to_summary_map.end() ) {

1084

p_one_summary = find_it->second;

1086

p_one_summary.

Reset

(

1087 new

GA::SOneGapLengthSummary(

1098 "one_gap_len_type_info"

);

1102

one_gap_len_type_info_attributes.

(

1105

one_gap_len_type_info_attributes.

(

1108

one_gap_len_type_info_attributes.

(

1112

one_gap_len_info.

(one_gap_len_type_info);

1115

gap_len_infos_node.

(one_gap_len_info);

1117

gap_info_root_node.

(gap_len_infos_node);

1125

cout <<

"SUMMARY:"

<< endl;

1127 bool

bAnyGapOfLenZero =

false

;

1129 const size_t

kDigitsInUint8 = numeric_limits<Uint8>::digits10;

1131

vecColInfos.

AddCol

(

"Gap Length"

, kDigitsInUint8,

1135 const

GA::EGapType eGapType = gap_type_name_it->first;

1136 string

pchGapName = gap_type_name_it->second;

1143 "#Seqs with "

+ pchGapName, kDigitsInUint8,

1146 "# of "

+ pchGapName, kDigitsInUint8,

1153 const xml::node

& one_gap_len_node = *one_gap_len_it;

1157 const

GA::TGapLength gap_len =

1158

to_uint8(find_attrib_attr_or_die(one_gap_len_info,

"len"

));

1159

table_printer << gap_len << CellEnd();

1160 if

( 0 == gap_len ) {

1161

bAnyGapOfLenZero =

true

;

1167 const xml::node

& gap_len_summary_node = *child_it;

1168 const Uint8

num_seqs = to_uint8(

1169

find_node_attr_or_die(

1170

gap_len_summary_node,

"num_seqs"

));

1171

table_printer << num_seqs << CellEnd();

1173 const Uint8

num_gaps = to_uint8(find_node_attr_or_die(

1174

gap_len_summary_node,

"num_gaps"

));

1175

table_printer << num_gaps << CellEnd();

1182 if

( bAnyGapOfLenZero ) {

1183

cout <<

"* Note: Gap of length zero means " 1184

"'completely unknown length'."

<< endl;

1205 const

GA::TGapLength gap_len = *all_gap_lens_it;

1207 xml::node

gap_seqs_one_len_node(

"gap_length_info"

);

1212 const

GA::EGapType eGapType = gap_type_name_it->first;

1213 const char

* pchGapName = gap_type_name_it->second;

1219 xml::node

gap_seqs_one_len_and_gap_type(

"gap_type_seq_ids"

);

1221 "gap_type"

, pchGapName);

1223 const

GA::TMapGapLengthToSeqIds & map_len_to_seq_ids =

1226

GA::TMapGapLengthToSeqIds::const_iterator find_seq_ids_it =

1227

map_len_to_seq_ids.find(gap_len);

1230 if

( find_seq_ids_it != map_len_to_seq_ids.end() ) {

1231 const

GA::TSetSeqIdConstRef & set_seq_id_const_ref =

1232

find_seq_ids_it->second;

1234

GA::TSetSeqIdConstRef, seq_id_ref_it, set_seq_id_const_ref)

1238 "seq_id"

, (*seq_id_ref_it)->AsFastaString().c_str());

1240

gap_seqs_one_len_and_gap_type.

push_back

(one_seq_node);

1244

gap_seqs_one_len_node.

push_back

(gap_seqs_one_len_and_gap_type);

1247

gap_seqs_root_node.

push_back

(gap_seqs_one_len_node);

1283

cout <<

"SEQ-IDS FOR EACH GAP-LENGTH:"

<< endl;

1286 const

GA::TGapLength iGapLength =

1287

to_uint8(find_node_attr_or_die(

1288

*gap_len_node_it,

"len"

));

1289

cout <<

"\tSeq-ids with a gap of length " 1290

<< iGapLength <<

':'

<< endl;

1293 const CTempString

pchGapName = find_node_attr_or_die(

1294

*gap_type_seq_ids_it,

"gap_type"

);

1295

cout <<

"\t\t"

<< pchGapName <<

":"

<< endl;

1297 if

( gap_type_seq_ids_it->size() ) {

1300

<< find_node_attr_or_die(*seq_info_it,

"seq_id"

)

1304

cout <<

"\t\t\t(NONE)"

<< endl;

1328 const

GA::EGapType eGapType = gap_type_it->first;

1329 const char

* pchGapName = gap_type_it->second;

1335 xml::node

histogram_node(

"histogram"

);

1338

histogram_node_attrs.

(

"gap_type"

, pchGapName);

1341 m_gapAnalysis

.GetGapHistogram(eGapType, num_bins, eHistAlgo));

1360

histogram_node.

(bin_node);

1363

hist_root_node.

(histogram_node);

1373 const size_t

kDigitsInUint8 = numeric_limits<Uint8>::digits10;

1378 const xml::node

& hist_node = *hist_node_it;

1381

hist_node_it->get_attributes();

1383 const char

* pchGapName = find_attrib_attr_or_die(

1384

hist_node_attrs,

"gap_type"

).data();

1387

vecColInfos.

AddCol

(

"Range"

, 1 + 2*kDigitsInUint8);

1388

vecColInfos.

AddCol

(

"Number in Range"

, kDigitsInUint8,

1392

cout <<

"HISTOGRAM FOR "

<< pchGapName <<

":"

<< endl;

1397

bin_node_it->get_attributes();

1399 const Uint8

start = to_uint8(find_attrib_attr_or_die(

1400

bin_node_attrs,

"start_inclusive"

));

1401 const Uint8

end = to_uint8(find_attrib_attr_or_die(

1402

bin_node_attrs,

"end_inclusive"

));

1403 const Uint8

num_appearances = to_uint8(find_attrib_attr_or_die(

1404

bin_node_attrs,

"num_appearances"

));

1406

table_printer << start <<

'-'

<< end << CellEnd();

1407

table_printer << num_appearances << CellEnd();

1417 const

SOutMessage &out_message,