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);  gobble(in);
  95                         if (seq.getName() != "") {
  96                                 statData data = getStats(seq.getAligned());
  97
  98                                 out << seq.getName() << '\t' << data.pound << '\t' << data.dash << '\t' << data.plus << '\t' << data.equal << '\t';
  99                                 out << data.loop << '\t' << data.tilde << '\t' << data.total << endl;
 100                         }
 101                 }
 102
 103                 in.close();
 104                 out.close();
 105
 106                 return 0;
 107         }
 108
 109         catch(exception& e) {
 110                 errorOut(e, "AlignCheckCommand", "execute");
 111                 exit(1);
 112         }
 113 }
 114 //**********************************************************************************************************************
 115 void AlignCheckCommand::readMap(){
 116         try {
 117
 118                 structMap.resize(1, 0);
 119                 ifstream in;
 120
 121                 openInputFile(mapfile, in);
 122
 123                 while(!in.eof()){
 124                         int position;
 125                         in >> position;
 126                         structMap.push_back(position);
 127                         gobble(in);
 128                 }
 129                 in.close();
 130
 131                 seqLength = structMap.size();
 132
 133
 134                 //check you make sure is structMap[10] = 380 then structMap[380] = 10.
 135                 for(int i=0;i<seqLength;i++){
 136                         if(structMap[i] != 0){
 137                                 if(structMap[structMap[i]] != i){
 138                                         mothurOut("Your map file contains an error:  line " + toString(i) + " does not match line " + toString(structMap[i]) + "."); mothurOutEndLine();
 139                                 }
 140                         }
 141                 }
 142
 143
 144         }
 145         catch(exception& e) {
 146                 errorOut(e, "AlignCheckCommand", "readMap");
 147                 exit(1);
 148         }
 149 }
 150 /**************************************************************************************************/
 151
 152 statData AlignCheckCommand::getStats(string sequence){
 153         try {
 154
 155                 statData data;
 156                 sequence = "*" + sequence; // need to pad the sequence so we can index it by 1
 157
 158                 int seqLength = sequence.length();
 159                 for(int i=1;i<seqLength;i++){
 160                         if(structMap[i] != 0){
 161                                 if(sequence[i] == 'A'){
 162                                         if(sequence[structMap[i]] == 'T')               {       data.tilde++;   }
 163                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   }
 164                                         else if(sequence[structMap[i]] == 'G')  {       data.equal++;   }
 165                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 166                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 167                                         data.total++;
 168                                 }
 169                                 else if(sequence[i] == 'T'){
 170                                         if(sequence[structMap[i]] == 'T')               {       data.plus++;    }
 171                                         else if(sequence[structMap[i]] == 'A')  {       data.tilde++;   }
 172                                         else if(sequence[structMap[i]] == 'G')  {       data.dash++;    }
 173                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 174                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 175                                         data.total++;
 176                                 }
 177                                 else if(sequence[i] == 'G'){
 178                                         if(sequence[structMap[i]] == 'T')               {       data.dash++;    }
 179                                         else if(sequence[structMap[i]] == 'A')  {       data.equal++;   }
 180                                         else if(sequence[structMap[i]] == 'G')  {       data.pound++;   }
 181                                         else if(sequence[structMap[i]] == 'C')  {       data.tilde++;   }
 182                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 183                                         data.total++;
 184                                 }
 185                                 else if(sequence[i] == 'C'){
 186                                         if(sequence[structMap[i]] == 'T')               {       data.pound++;   }
 187                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   }
 188                                         else if(sequence[structMap[i]] == 'G')  {       data.tilde++;   }
 189                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   }
 190                                         else if(sequence[structMap[i]] == '-')  {       data.pound++;   }
 191                                         data.total++;
 192                                 }
 193                                 else if(sequence[i] == '-'){
 194                                         if(sequence[structMap[i]] == 'T')               {       data.pound++;   data.total++;   }
 195                                         else if(sequence[structMap[i]] == 'A')  {       data.pound++;   data.total++;   }
 196                                         else if(sequence[structMap[i]] == 'G')  {       data.pound++;   data.total++;   }
 197                                         else if(sequence[structMap[i]] == 'C')  {       data.pound++;   data.total++;   }
 198                                         else if(sequence[structMap[i]] == '-')  {               /*donothing*/                           }
 199                                 }
 200                         }
 201                         else if(isalnum(sequence[i])){
 202                                 data.loop++;
 203                                 data.total++;
 204                         }
 205                 }
 206                 return data;
 207
 208         }
 209         catch(exception& e) {
 210                 errorOut(e, "AlignCheckCommand", "getStats");
 211                 exit(1);
 212         }
 213 }
 214
 215
 216 //**********************************************************************************************************************