5 * Created by Pat Schloss on 5/30/09.
6 * Copyright 2009 Patrick D. Schloss. All rights reserved.
10 #include "seqsummarycommand.h"
12 //***************************************************************************************************************
14 SeqSummaryCommand::SeqSummaryCommand(){
16 globaldata = GlobalData::getInstance();
18 if(globaldata->getFastaFile() != "") { readSeqs = new ReadFasta(globaldata->inputFileName); }
19 else if(globaldata->getNexusFile() != "") { readSeqs = new ReadNexus(globaldata->inputFileName); }
20 else if(globaldata->getClustalFile() != "") { readSeqs = new ReadClustal(globaldata->inputFileName); }
21 else if(globaldata->getPhylipFile() != "") { readSeqs = new ReadPhylip(globaldata->inputFileName); }
24 db = readSeqs->getDB();
28 cout << "Standard Error: " << e.what() << " has occurred in the SeqCoordCommand class Function SeqCoordCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
32 cout << "An unknown error has occurred in the SeqCoordCommand class function SeqCoordCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
37 //***************************************************************************************************************
39 SeqSummaryCommand::~SeqSummaryCommand(){
42 //***************************************************************************************************************
44 int SeqSummaryCommand::execute(){
48 string summaryFile = getRootName(globaldata->inputFileName) + "summary";
49 openOutputFile(summaryFile, outfile);
51 vector<int> startPosition(numSeqs, 0);
52 vector<int> endPosition(numSeqs, 0);
53 vector<int> seqLength(numSeqs, 0);
54 vector<int> ambigBases(numSeqs, 0);
55 vector<int> longHomoPolymer(numSeqs, 0);
57 if(db->get(0).getIsAligned() == 1){
58 outfile << "seqname\tstart\tend\tlength\tambiguities\tlonghomopolymer" << endl;
59 for(int i = 0; i < numSeqs; i++) {
60 Sequence current = db->get(i);
61 startPosition[i] = current.getStartPos();
62 endPosition[i] = current.getEndPos();
63 seqLength[i] = current.getNumBases();
64 ambigBases[i] = current.getAmbigBases();
65 longHomoPolymer[i] = current.getLongHomoPolymer();
66 outfile << current.getName() << '\t' << startPosition[i] << '\t' << endPosition[i] << '\t' << seqLength[i] << '\t' << ambigBases[i] << '\t' << longHomoPolymer[i] << endl;
70 outfile << "seqname\tlength\tambiguities\tlonghomopolymer" << endl;
71 for(int i=0;i<numSeqs;i++){
72 Sequence current = db->get(i);
73 seqLength[i] = current.getNumBases();
74 ambigBases[i] = current.getAmbigBases();
75 longHomoPolymer[i] = current.getLongHomoPolymer();
76 outfile << current.getName() << '\t' << seqLength[i] << '\t' << ambigBases[i] << '\t' << longHomoPolymer[i] << endl;
80 sort(seqLength.begin(), seqLength.end());
81 sort(ambigBases.begin(), ambigBases.end());
82 sort(longHomoPolymer.begin(), longHomoPolymer.end());
84 int median = int(numSeqs * 0.500);
85 int lowestPtile = int(numSeqs * 0.025);
86 int lowPtile = int(numSeqs * 0.250);
87 int highPtile = int(numSeqs * 0.750);
88 int highestPtile = int(numSeqs * 0.975);
89 int max = numSeqs - 1;
92 if(db->get(0).getIsAligned() == 1){
93 sort(startPosition.begin(), startPosition.end());
94 sort(endPosition.begin(), endPosition.end());
96 cout << "\t\tStart\tEnd\tLength\tN's\tPolymer" << endl;
97 cout << "Minimum:\t" << startPosition[0] << '\t' << endPosition[0] << '\t' << seqLength[0] << '\t' << ambigBases[0] << '\t' << longHomoPolymer[0] << endl;
98 cout << "2.5%-tile:\t" << startPosition[lowestPtile] << '\t' << endPosition[lowestPtile] << '\t' << seqLength[lowestPtile] << '\t' << ambigBases[lowestPtile] << '\t' << longHomoPolymer[lowestPtile] << endl;
99 cout << "25%-tile:\t" << startPosition[lowPtile] << '\t' << endPosition[lowPtile] << '\t' << seqLength[lowPtile] << '\t' << ambigBases[lowPtile] << '\t' << longHomoPolymer[lowPtile] << endl;
100 cout << "Median: \t" << startPosition[median] << '\t' << endPosition[median] << '\t' << seqLength[median] << '\t' << ambigBases[median] << '\t' << longHomoPolymer[median] << endl;
101 cout << "75%-tile:\t" << startPosition[highPtile] << '\t' << endPosition[highPtile] << '\t' << seqLength[highPtile] << '\t' << ambigBases[highPtile] << '\t' << longHomoPolymer[highPtile] << endl;
102 cout << "97.5%-tile:\t" << startPosition[highestPtile] << '\t' << endPosition[highestPtile] << '\t' << seqLength[highestPtile] << '\t' << ambigBases[highestPtile] << '\t' << longHomoPolymer[highestPtile] << endl;
103 cout << "Maximum:\t" << startPosition[max] << '\t' << endPosition[max] << '\t' << seqLength[max] << '\t' << ambigBases[max] << '\t' << longHomoPolymer[max] << endl;
106 cout << "\t\tLength\tN's\tPolymer" << endl;
107 cout << "Minimum:\t" << seqLength[0] << '\t' << ambigBases[0] << '\t' << longHomoPolymer[0] << endl;
108 cout << "2.5%-tile:\t" << seqLength[lowestPtile] << '\t' << ambigBases[lowestPtile] << '\t' << longHomoPolymer[lowestPtile] << endl;
109 cout << "25%-tile:\t" << seqLength[lowPtile] << '\t' << ambigBases[lowPtile] << '\t' << longHomoPolymer[lowPtile] << endl;
110 cout << "Median: \t" << seqLength[median] << '\t' << ambigBases[median] << '\t' << longHomoPolymer[median] << endl;
111 cout << "75%-tile:\t"<< seqLength[highPtile] << '\t' << ambigBases[highPtile] << '\t' << longHomoPolymer[highPtile] << endl;
112 cout << "97.5%-tile:\t"<< seqLength[highestPtile] << '\t' << ambigBases[highestPtile] << '\t' << longHomoPolymer[highestPtile] << endl;
113 cout << "Maximum:\t" << seqLength[max] << '\t' << ambigBases[max] << '\t' << longHomoPolymer[max] << endl;
115 cout << "# of Seqs:\t" << numSeqs << endl;
119 catch(exception& e) {
120 cout << "Standard Error: " << e.what() << " has occurred in the FilterSeqsCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
124 cout << "An unknown error has occurred in the FilterSeqsCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
130 //***************************************************************************************************************