secondarystructurecommand.cpp

   1 /*
   2  *  secondarystructurecommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/18/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "secondarystructurecommand.h"
  11 #include "sequence.hpp"
  12
  13 //**********************************************************************************************************************
  14
  15 AlignCheckCommand::AlignCheckCommand(string option){
  16         try {
  17                 abort = false;
  18
  19                 //allow user to run help
  20                 if(option == "help") { help(); abort = true; }
  21
  22                 else {
  23                         //valid paramters for this command
  24                         string Array[] =  {"fasta","map"};
  25                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  26
  27                         OptionParser parser(option);
  28                         map<string,string> parameters = parser.getParameters();
  29
  30                         ValidParameters validParameter;
  31
  32                         //check to make sure all parameters are valid for command
  33                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  34                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  35                         }
  36
  37                         //check for required parameters
  38                         mapfile = validParameter.validFile(parameters, "map", true);
  39                         if (mapfile == "not open") { abort = true; }
  40                         else if (mapfile == "not found") {  mapfile = "";  mothurOut("You must provide an map file."); mothurOutEndLine(); abort = true; }
  41
  42                         fastafile = validParameter.validFile(parameters, "fasta", true);
  43                         if (fastafile == "not open") { abort = true; }
  44                         else if (fastafile == "not found") {  fastafile = "";  mothurOut("You must provide an fasta file."); mothurOutEndLine(); abort = true;  }
  45
  46                 }
  47
  48         }
  49         catch(exception& e) {
  50                 errorOut(e, "AlignCheckCommand", "RemoveSeqsCommand");
  51                 exit(1);
  52         }
  53 }
  54 //**********************************************************************************************************************
  55
  56 void AlignCheckCommand::help(){
  57         try {
  58                 mothurOut("The align.check command reads a fasta file and map file.\n");
  59                 mothurOut("It outputs a file containing the secondary structure matches in the .align.check file.\n");
  60                 mothurOut("The align.check command parameters are fasta and map, both are required.\n");
  61                 mothurOut("The align.check command should be in the following format: align.check(fasta=yourFasta, map=yourMap).\n");
  62                 mothurOut("Example align.check(map=silva.ss.map, fasta=amazon.fasta).\n");
  63                 mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
  64         }
  65         catch(exception& e) {
  66                 errorOut(e, "AlignCheckCommand", "help");
  67                 exit(1);
  68         }
  69 }
  70
  71 //**********************************************************************************************************************
  72
  73 int AlignCheckCommand::execute(){
  74         try {
  75
  76                 if (abort == true) { return 0; }
  77
  78                 //get secondary structure info.
  79                 readMap();
  80
  81                 ifstream in;
  82                 openInputFile(fastafile, in);
  83
  84                 ofstream out;
  85                 string outfile = getRootName(fastafile) + "align.check";
  86                 openOutputFile(outfile, out);
  87
  88                 out << "name" << '\t' << "pound" << '\t' << "dash" << '\t' << "plus" << '\t' << "equal" << '\t';
  89                 out << "loop" << '\t' << "tilde" << '\t' << "total" << endl;
  90
  91
  92                 while(!in.eof()){
  93
  94                         Sequence seq(in);
  95                         statData data = getStats(seq.getAligned());
  96
  97                         out << seq.getName() << '\t' << data.pound << '\t' << data.dash << '\t' << data.plus << '\t' << data.equal << '\t';
  98                         out << data.loop << '\t' << data.tilde << '\t' << data.total << endl;
  99                 }
 100
 101                 in.close();
 102                 out.close();
 103
 104                 return 0;
 105         }
 106
 107         catch(exception& e) {
 108                 errorOut(e, "AlignCheckCommand", "execute");
 109                 exit(1);
 110         }
 111 }
 112 //**********************************************************************************************************************
 113 void AlignCheckCommand::readMap(){
 114         try {
 115
 116                 structMap.resize(1, 0);
 117                 ifstream in;
 118
 119                 openInputFile(mapfile, in);
 120
 121                 while(!in.eof()){
 122                         int position;
 123                         in >> position;
 124                         structMap.push_back(position);
 125                         gobble(in);
 126                 }
 127                 in.close();
 128
 129                 seqLength = structMap.size();
 130
 131
 132                 //check you make sure is structMap[10] = 380 then structMap[380] = 10.
 133                 for(int i=0;i<seqLength;i++){
 134                         if(structMap[i] != 0){
 135                                 if(structMap[structMap[i]] != i){
 136                                         mothurOut("Your map file contains an error:  line " + toString(i) + " does not match line " + toString(structMap[i]) + "."); mothurOutEndLine();
 137                                 }
 138                         }
 139                 }
 140
 141
 142         }
 143         catch(exception& e) {
 144                 errorOut(e, "AlignCheckCommand", "readMap");
 145                 exit(1);
 146         }
 147 }
 148 /**************************************************************************************************/
 149
 150 statData AlignCheckCommand::getStats(string sequence){
 151         try {
 152
 153                 statData data;
 154                 sequence = "*" + sequence; // need to pad the sequence so we can index it by 1
 155
 156                 int seqLength = sequence.length();
 157                 for(int i=1;i<seqLength;i++){
 158                         if(structMap[i] != 0){
 159                                 if(sequence[i] == 'A'){
 160                                         if(sequence[structMap[i]] == 'T')               {       data.tilde++;   }
 161                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   }
 162                                         else if(sequence[structMap[i]] == 'G')  {       data.equal++;   }
 163                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 164                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 165                                         data.total++;
 166                                 }
 167                                 else if(sequence[i] == 'T'){
 168                                         if(sequence[structMap[i]] == 'T')               {       data.plus++;    }
 169                                         else if(sequence[structMap[i]] == 'A')  {       data.tilde++;   }
 170                                         else if(sequence[structMap[i]] == 'G')  {       data.dash++;    }
 171                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 172                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 173                                         data.total++;
 174                                 }
 175                                 else if(sequence[i] == 'G'){
 176                                         if(sequence[structMap[i]] == 'T')               {       data.dash++;    }
 177                                         else if(sequence[structMap[i]] == 'A')  {       data.equal++;   }
 178                                         else if(sequence[structMap[i]] == 'G')  {       data.pound++;   }
 179                                         else if(sequence[structMap[i]] == 'C')  {       data.tilde++;   }
 180                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 181                                         data.total++;
 182                                 }
 183                                 else if(sequence[i] == 'C'){
 184                                         if(sequence[structMap[i]] == 'T')               {       data.pound++;   }
 185                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   }
 186                                         else if(sequence[structMap[i]] == 'G')  {       data.tilde++;   }
 187                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 188                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 189                                         data.total++;
 190                                 }
 191                                 else if(sequence[i] == '-'){
 192                                         if(sequence[structMap[i]] == 'T')               {       data.pound++;   data.total++;   }
 193                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   data.total++;   }
 194                                         else if(sequence[structMap[i]] == 'G')  {       data.pound++;   data.total++;   }
 195                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   data.total++;   }
 196                                         else if(sequence[structMap[i]] == '-')  {               /*donothing*/                           }
 197                                 }
 198                         }
 199                         else if(isalnum(sequence[i])){
 200                                 data.loop++;
 201                                 data.total++;
 202                         }
 203                 }
 204                 return data;
 205
 206         }
 207         catch(exception& e) {
 208                 errorOut(e, "AlignCheckCommand", "getStats");
 209                 exit(1);
 210         }
 211 }
 212
 213
 214 //**********************************************************************************************************************