seqsummarycommand.cpp

   1 /*
   2  *  seqcoordcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 5/30/09.
   6  *  Copyright 2009 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10 #include "seqsummarycommand.h"
  11
  12 //***************************************************************************************************************
  13
  14 SeqSummaryCommand::SeqSummaryCommand(){
  15         try {
  16                 globaldata = GlobalData::getInstance();
  17
  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);   }
  22
  23                 readSeqs->read();
  24                 db = readSeqs->getDB();
  25                 numSeqs = db->size();
  26         }
  27         catch(exception& e) {
  28                 cout << "Standard Error: " << e.what() << " has occurred in the SeqCoordCommand class Function SeqCoordCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  29                 exit(1);
  30         }
  31         catch(...) {
  32                 cout << "An unknown error has occurred in the SeqCoordCommand class function SeqCoordCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  33                 exit(1);
  34         }
  35 }
  36
  37 //***************************************************************************************************************
  38
  39 SeqSummaryCommand::~SeqSummaryCommand(){
  40 }
  41
  42 //***************************************************************************************************************
  43
  44 int SeqSummaryCommand::execute(){
  45         try{
  46
  47                 ofstream outfile;
  48                 string summaryFile = getRootName(globaldata->inputFileName) + "summary";
  49                 openOutputFile(summaryFile, outfile);
  50
  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);
  56
  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;
  67                         }
  68                 }
  69                 else{
  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;
  77                         }
  78                 }
  79
  80                 sort(seqLength.begin(), seqLength.end());
  81                 sort(ambigBases.begin(), ambigBases.end());
  82                 sort(longHomoPolymer.begin(), longHomoPolymer.end());
  83
  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;
  90
  91                 cout << endl;
  92                 if(db->get(0).getIsAligned() == 1){
  93                         sort(startPosition.begin(), startPosition.end());
  94                         sort(endPosition.begin(), endPosition.end());
  95
  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;
 104                 }
 105                 else{
 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;
 114                 }
 115                 cout << "# of Seqs:\t" << numSeqs << endl;
 116
 117                 return 0;
 118         }
 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";
 121                 exit(1);
 122         }
 123         catch(...) {
 124                 cout << "An unknown error has occurred in the FilterSeqsCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 125                 exit(1);
 126         }
 127
 128 }
 129
 130 //***************************************************************************************************************
 131
 132