- //read in seqs and store in vector
- while(!inFASTA.eof()){
- Sequence current(inFASTA);
-
- if (current.getAligned() == "") { current.setAligned(current.getUnaligned()); }
-
- seqs.push_back(current);
-
- gobble(inFASTA);
- }
- inFASTA.close();
-
- }
- catch(exception& e) {
- errorOut(e, "ChimeraSeqsCommand", "readSeqs");
- exit(1);
- }
-}
-
-/***************************************************************************************************************/
-int ChimeraSeqsCommand::createSparseMatrix(int startSeq, int endSeq, SparseMatrix* sparse, vector<Sequence> s){
- try {
-
- for(int i=startSeq; i<endSeq; i++){
-
- for(int j=0;j<i;j++){
-
- distCalculator->calcDist(s[i], s[j]);
- float dist = distCalculator->getDist();
-
- PCell temp(i, j, dist);
- sparse->addCell(temp);
-
- }
- }
-
-
- return 1;
- }
- catch(exception& e) {
- errorOut(e, "ChimeraSeqsCommand", "createSparseMatrix");
- exit(1);
- }
-}
-/***************************************************************************************************************/
-void ChimeraSeqsCommand::generatePreferences(vector<SeqMap> left, vector<SeqMap> right, int mid){
- try {
-
- float dme = 0.0;
- SeqMap::iterator itR;
- SeqMap::iterator itL;
-
- //initialize pref[i]
- for (int i = 0; i < pref.size(); i++) {
- pref[i].score[1] = 0.0;
- pref[i].closestLeft[1] = 100000.0;
- pref[i].closestRight[1] = 100000.0;
- pref[i].leftParent[1] = "";
- pref[i].rightParent[1] = "";
- }
-
- for (int i = 0; i < left.size(); i++) {
-
- SeqMap currentLeft = left[i]; //example i = 3; currentLeft is a map of 0 to the distance of sequence 3 to sequence 0,
- // 1 to the distance of sequence 3 to sequence 1,
- // 2 to the distance of sequence 3 to sequence 2.
- SeqMap currentRight = right[i]; // same as left but with distances on the right side.
-
- for (int j = 0; j < i; j++) {
-
- itL = currentLeft.find(j);
- itR = currentRight.find(j);
-cout << " i = " << i << " j = " << j << " distLeft = " << itL->second << endl;
-cout << " i = " << i << " j = " << j << " distright = " << itR->second << endl;
-
- //if you can find this entry update the preferences
- if ((itL != currentLeft.end()) && (itR != currentRight.end())) {
-
- if (!correction) {
- pref[i].score[1] += abs((itL->second - itR->second));
- pref[j].score[1] += abs((itL->second - itR->second));
-cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
-cout << "abs = " << abs((itL->second - itR->second)) << endl;
-cout << i << " score = " << pref[i].score[1] << endl;
-cout << j << " score = " << pref[j].score[1] << endl;
- }else {
- pref[i].score[1] += abs((sqrt(itL->second) - sqrt(itR->second)));
- pref[j].score[1] += abs((sqrt(itL->second) - sqrt(itR->second)));
-cout << "left " << i << " " << j << " = " << itL->second << " right " << i << " " << j << " = " << itR->second << endl;
-cout << "abs = " << abs((sqrt(itL->second) - sqrt(itR->second))) << endl;
-cout << i << " score = " << pref[i].score[1] << endl;
-cout << j << " score = " << pref[j].score[1] << endl;
- }
-cout << "pref[" << i << "].closestLeft[1] = " << pref[i].closestLeft[1] << " parent = " << pref[i].leftParent[1] << endl;
- //are you the closest left sequence
- if (itL->second < pref[i].closestLeft[1]) {
-
- pref[i].closestLeft[1] = itL->second;
- pref[i].leftParent[1] = seqs[j].getName();
-cout << "updating closest left to " << pref[i].leftParent[1] << endl;
- }
-cout << "pref[" << j << "].closestLeft[1] = " << pref[j].closestLeft[1] << " parent = " << pref[j].leftParent[1] << endl;
- if (itL->second < pref[j].closestLeft[1]) {
- pref[j].closestLeft[1] = itL->second;
- pref[j].leftParent[1] = seqs[i].getName();
-cout << "updating closest left to " << pref[j].leftParent[1] << endl;
- }
-
- //are you the closest right sequence
- if (itR->second < pref[i].closestRight[1]) {
- pref[i].closestRight[1] = itR->second;
- pref[i].rightParent[1] = seqs[j].getName();
- }
- if (itR->second < pref[j].closestRight[1]) {
- pref[j].closestRight[1] = itR->second;
- pref[j].rightParent[1] = seqs[i].getName();
- }
-
- }
- }
-
- }
-
-
-
- //calculate the dme
- int count0 = 0;
- for (int i = 0; i < pref.size(); i++) { dme += pref[i].score[1]; if (pref[i].score[1] == 0.0) { count0++; } }
-
- float expectedPercent = 1 / (float) (pref.size() - count0);
-cout << endl << "dme = " << dme << endl;
- //recalculate prefernences based on dme
- for (int i = 0; i < pref.size(); i++) {
-cout << "unadjusted pref " << i << " = " << pref[i].score[1] << endl;
- // gives the actual percentage of the dme this seq adds
- pref[i].score[1] = pref[i].score[1] / dme;
-
- //how much higher or lower is this than expected
- pref[i].score[1] = pref[i].score[1] / expectedPercent;
-
- //so a non chimeric sequence would be around 1, and a chimeric would be signifigantly higher.
-cout << "adjusted pref " << i << " = " << pref[i].score[1] << endl;
- }