chimera.cpp

   1 /*
   2  *  chimera.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 8/11/09.
   6  *  Copyright 2009 Schloss Lab Umass Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "chimera.h"
  11
  12 //***************************************************************************************************************
  13 //this is a vertical soft filter
  14 string Chimera::createFilter(vector<Sequence*> seqs, float t) {
  15         try {
  16                 filterString = "";
  17                 int threshold = int (t * seqs.size());
  18 //cout << "threshhold = " << threshold << endl;
  19
  20                 vector<int> gaps;       gaps.resize(seqs[0]->getAligned().length(), 0);
  21                 vector<int> a;          a.resize(seqs[0]->getAligned().length(), 0);
  22                 vector<int> t;          t.resize(seqs[0]->getAligned().length(), 0);
  23                 vector<int> g;          g.resize(seqs[0]->getAligned().length(), 0);
  24                 vector<int> c;          c.resize(seqs[0]->getAligned().length(), 0);
  25
  26                 filterString = (string(seqs[0]->getAligned().length(), '1'));
  27
  28                 //for each sequence
  29                 for (int i = 0; i < seqs.size(); i++) {
  30
  31                         if (m->control_pressed) { return filterString; }
  32
  33                         string seqAligned = seqs[i]->getAligned();
  34
  35                         for (int j = 0; j < seqAligned.length(); j++) {
  36                                 //if this spot is a gap
  37                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
  38                                 else if (toupper(seqAligned[j]) == 'A')                                 {       a[j]++;         }
  39                                 else if (toupper(seqAligned[j]) == 'T')                                 {       t[j]++;         }
  40                                 else if (toupper(seqAligned[j]) == 'G')                                 {       g[j]++;         }
  41                                 else if (toupper(seqAligned[j]) == 'C')                                 {       c[j]++;         }
  42                         }
  43                 }
  44
  45                 //zero out spot where all sequences have blanks
  46                 //zero out spot where all sequences have blanks
  47                 int numColRemoved = 0;
  48                 for(int i = 0;i < seqs[0]->getAligned().length(); i++){
  49
  50                         if (m->control_pressed) { return filterString; }
  51
  52                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
  53
  54                         else if (((a[i] < threshold) && (t[i] < threshold) && (g[i] < threshold) && (c[i] < threshold))) {      filterString[i] = '0';  numColRemoved++;  }
  55                         //cout << "a = " << a[i] <<  " t = " << t[i] <<  " g = " << g[i] <<  " c = " << c[i] << endl;
  56                 }
  57
  58                 m->mothurOut("Filter removed " + toString(numColRemoved) + " columns.");  m->mothurOutEndLine();
  59                 return filterString;
  60         }
  61         catch(exception& e) {
  62                 m->errorOut(e, "Chimera", "createFilter");
  63                 exit(1);
  64         }
  65 }
  66 //***************************************************************************************************************
  67 map<int, int> Chimera::runFilter(Sequence* seq) {
  68         try {
  69                 map<int, int> maskMap;
  70                 string seqAligned = seq->getAligned();
  71                 string newAligned = "";
  72                 int count = 0;
  73
  74                 for (int j = 0; j < seqAligned.length(); j++) {
  75                         //if this spot is a gap
  76                         if (filterString[j] == '1') {
  77                                 newAligned += seqAligned[j];
  78                                 maskMap[count] = j;
  79                                 count++;
  80                         }
  81                 }
  82
  83                 seq->setAligned(newAligned);
  84
  85                 return maskMap;
  86         }
  87         catch(exception& e) {
  88                 m->errorOut(e, "Chimera", "runFilter");
  89                 exit(1);
  90         }
  91 }
  92 //***************************************************************************************************************
  93 vector<Sequence*> Chimera::readSeqs(string file) {
  94         try {
  95
  96                 m->mothurOut("Reading sequences... "); cout.flush();
  97                 ifstream in;
  98                 openInputFile(file, in);
  99
 100                 vector<Sequence*> container;
 101                 int count = 0;
 102                 length = 0;
 103                 unaligned = false;
 104
 105                 //read in seqs and store in vector
 106                 while(!in.eof()){
 107
 108                         if (m->control_pressed) { return container; }
 109
 110                         Sequence* current = new Sequence(in);  gobble(in);
 111
 112                         if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
 113                         else if (length != current->getAligned().length()) { //seqs are unaligned
 114                                 unaligned = true;
 115                         }
 116
 117                         if (current->getName() != "") {  container.push_back(current);  }
 118                 }
 119
 120                 in.close();
 121                 m->mothurOut("Done."); m->mothurOutEndLine();
 122
 123                 return container;
 124         }
 125         catch(exception& e) {
 126                 m->errorOut(e, "Chimera", "readSeqs");
 127                 exit(1);
 128         }
 129 }
 130 //***************************************************************************************************************
 131 void Chimera::setMask(string filename) {
 132         try {
 133
 134                 if (filename == "default") {
 135                         //default is from wigeon  236627 EU009184.1 Shigella dysenteriae str. FBD013
 136                         seqMask = ".....................................................................................................AAATTGAAGAGTTT-GA--T-CA-T-G-GCTC-AG-AT-TGAA-C-GC--TGG-C--G-GC-A-GG--C----C-T--AACACA-T-GC-A-AGT-CGA-A-CG----------G-TAA-CA-G----------------------------GAAG-A-AG----------------------------------------------------CTT-G----------------------------------------------------------------------------------CT-TCTTT----------------G-CT--G--AC--G--AG-T-GG-C-GG-A--C-------------GGG-TGAGT-A--AT-GT-C-T-G-GG---A-A--A-CT-G--C-C-TGA--TG-G------------------------------------------------------------------A-GG----GGG-AT-AA-CTA-------------------------C-T-G-----------------------GAA-A---CGG-TAG-CTAA-TA---CC-G--C-AT-A----------A--------------------C-------------------------------------GT-C-----------------------------------------------------------------------------------------------------------------------G-CA-A--------------------------------------------------------------------------------------------------------------------------------------G-A-C---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CAAA--G-A-G-GG-----G--GA-C-CT--------------------------------------------------------------------------------------------------------------------TCG-G----------------------------------------------------------------------------------------------------------------------G----CC-TC--T---T-G--------------C----C-A---T-CG-G---AT---G-T-----G-CCC-AGA--T-GGG--A------TT--A--G-CT-A----G---TAGG-T-G-GG-G-T----AAC-GG-C-T-C-ACCT--A-GG-C-G--A-CG-A------------TCC-C-T------AG-CT-G-G-TCT-G-AG----A--GG-AT--G-AC-C-AG-CCAC-A-CTGGA--A-C-TG-A-GA-C-AC-G-G-TCCAGA-CTCC-TAC-G--G-G-A-G-GC-A-GC-A-G-TG---GG-G-A-ATA-TTGCA-C-AA-T-GG--GC-GC-A----A-G-CC-T-GA-TG-CA-GCCA-TGCC-G-CG-T---G-T-A--T--GA-A-G--A--A-G-G-CC-----TT-CG---------G-G-T-T-G-T--A---AA-G-TAC--------TT-TC-A-G--C-GGG----GA-G--G---AA-GGGA---GTAA-AG----T--T--AA-T---A----C-----CT-T-TGC-TCA-TT-GA-CG-TT-A-C-CC-G-CA-G---------AA-----------GAAGC-ACC-GG-C-TAA---C--T-CCGT--GCCA--G-C---A--GCCG---C-GG--TA-AT--AC---GG-AG-GGT-GCA-A-G-CG-TTAA-T-CGG-AA-TT-A--C-T--GGGC-GTA----AA-GCGC-AC--G-CA-G-G-C-G------------G--T-TT-G-T-T-AA----G-T-C-A---G-ATG-TG-A-AA-TC--CC-CGG-G--------------------------------------------------------------------CT-C-AA-------------------------------------------------------------------------CC-T-G-GG-AA-C----T-G-C-A-T-C--------T--GA-T-A-C-T-G-GCA--A-G-C---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-T-G-A-G-T-C-----T-CG--TA-G-A------------G-GG-G-GG-T----AG--AATT-CCA-G-GT--GT-A-GCG-GTGAAA-TG-CGT-AGAG-A-TC-T-GGA--GG-A-AT-A-CC-GG--T--G--GC-GAA-G--G-C---G----G--C-C-CCCTG------G-AC-GA--------------------------------------------------------------AG-A-C-T--GA--CG-----CT-CA-GG--T-G-CGA--AA-G-C--------------G-TGGG-GAG-C-A-AACA--GG-ATTA-G-ATA-C-----CC-T-G-GTA-G-T----C-CA--C-G-CCG-T-AAA--C-GATG-TC--GA-CT---------T-GG--A--G-G-TT-G-TG-C--C--------------------------------------------------------------------------------------CTT-GA--------------------------------------------------------------------------------------------------------------------------------------------------G-G-C-GT--G-G-C-T-TC-C------GG--A----GC-TAA--CG-C-G-T--T--AA-GT--C----G-ACC-GCC-T-G-GG-GAG-TA---CGG-----C-C--G-C-A-A-GGT-T--AAA-ACTC-AAA---------TGAA-TTG-ACGGG-G-G-CCCG----C-A--C-A-A-GCG-GT-G--G--AG-CA-T--GT-GGT-TT-AATT-C-G-ATG-CAAC-G-CG-A-AG-A-A-CC-TT-A-CC-TGGTC-TT-G-AC-A-T-C--------------CAC-G-G-------------A-AG-T-T-T--TC--A-GA-G-A-T--G-A-G--A-A-T-G--T-G-----CC-------------------------------------T--TC-G------------------------------------------GG----A----A---CC-GTG---A--GA---------------------------------------------------C-A-G-G-T-GCTG-CA-TGG-CT--GTC-GTC-A-GC-TC---G-TG-TT-G--TGA-AA-TGT-T-GG-G-TT-AA-GT-CCCGC-AA--------C-GAG-CGC-A-ACC-C-T-TA--TC--C-TTTG--T-T-G-C-C---AG-C-G-----G-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------TCC------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------GG---C----C-G------------G----G---A-A--CT---------------C-A-A-A-G-GA-G--AC-T-G-CCA--G-T------------------------------------G-A---TAA----------------------------------A-C-T-G--G-A-GG-A--AGG-T--GGGG-A-TGAC-GTC--AAGT-C---ATC-A-T-G-G-C-C-CTT----AC-G--AC-C-A-GG-GC-TA-CAC-ACGTG-C--TA--CAATG---G-CGCA-T-A--C-AAA-GA-GA--------------------------------------------------------------------------------------------------A-G-C-G-A--C-CTCG-C--G---------------------------------------A-GA-G-C-----------A--A-G-CG---G----------A--CCT-C------A-T-AAAGT-GC-G-T-C-G-TAG-TCC--------GGA-T-TGGAG-TC--T-GCAA-CT-C-------------------------------------------------------------------------------------------------G-ACTCC-A-T-G-AA-G-TC-GGAAT-CG-C-TA--G-TA-AT-C-G-T----GGA-TC-A-G--A------AT--GCC-AC-G-GT-G-AAT-ACGT-T-CCCGGGCCT-TGTA----CACACCG-CCC-GTC-----A---CA--CCA-TG-GG-A--G---TGG-G-TT-GC-AAA--A-GAA------G--T-AGG-TA-G-C-T-T-AA-C-C--------------------------------------------------------------TT----C-------------------------------------------------------------------------------------------------G--GG-A--GG-G--C---GC-TTA--CC--ACT-T----T-GTG-AT-TCA------------------------TG--ACT-GGGG-TG-AAG-TCGTAACAA-GGTAA-CCGT-AGGGGAA-CCTG-CGGT-TGGATCACCTCCTTA................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................";
 137                 }else if (filename == "") {  //do nothing
 138                         seqMask = "";
 139                 }else{
 140                         ifstream infile;
 141                         openInputFile(filename, infile);
 142
 143                         while (!infile.eof()) {
 144                                 Sequence temp(infile);
 145                                 seqMask = temp.getAligned();
 146
 147                                 gobble(infile);
 148                         }
 149
 150                         infile.close();
 151                 }
 152         }
 153         catch(exception& e) {
 154                 m->errorOut(e, "Chimera", "setMask");
 155                 exit(1);
 156         }
 157 }
 158 //***************************************************************************************************************
 159
 160 vector< vector<float> > Chimera::readQuantiles() {
 161         try {
 162
 163                 ifstream in;
 164                 openInputFile(quanfile, in);
 165
 166                 vector< vector<float> > quan;
 167                 vector <float> temp; temp.resize(6, 0);
 168
 169                 //to fill 0
 170                 quan.push_back(temp);
 171
 172                 int num; float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine;
 173
 174                 while(!in.eof()){
 175
 176                         in >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine;
 177
 178                         temp.clear();
 179
 180                         temp.push_back(ten);
 181                         temp.push_back(twentyfive);
 182                         temp.push_back(fifty);
 183                         temp.push_back(seventyfive);
 184                         temp.push_back(ninetyfive);
 185                         temp.push_back(ninetynine);
 186
 187                         quan.push_back(temp);
 188
 189                         gobble(in);
 190                 }
 191
 192                 in.close();
 193                 return quan;
 194
 195         }
 196         catch(exception& e) {
 197                 m->errorOut(e, "Chimera", "readQuantiles");
 198                 exit(1);
 199         }
 200 }
 201 //***************************************************************************************************************
 202 Sequence* Chimera::getSequence(string name) {
 203         try{
 204                 Sequence* temp;
 205
 206                 //look through templateSeqs til you find it
 207                 int spot = -1;
 208                 for (int i = 0; i < templateSeqs.size(); i++) {
 209                         if (name == templateSeqs[i]->getName()) {
 210                                 spot = i;
 211                                 break;
 212                         }
 213                 }
 214
 215                 if(spot == -1) { m->mothurOut("Error: Could not find sequence."); m->mothurOutEndLine(); return NULL; }
 216
 217                 temp = new Sequence(templateSeqs[spot]->getName(), templateSeqs[spot]->getAligned());
 218
 219                 return temp;
 220         }
 221         catch(exception& e) {
 222                 m->errorOut(e, "Chimera", "getSequence");
 223                 exit(1);
 224         }
 225 }
 226 //***************************************************************************************************************
 227
 228
 229
 230