+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-/**
- * Output routines ported from clustalx C code in 'interface.c'.
- *
- * Changes:
- *
- * 18-04-07,Nigel Brown(EMBL): clustalx code used (firsRes, length) convention
- * to define output columns, while the port uses (firsRes, lastRes). Fixed
- * problems in all 7 output routines where the conventions were mixed giving
- * over-long output blocks.
- *
- * 18-7-07: Mark (UCD), made changes to fastaOut, clustalOut, nbrfOut and gdeOut
- *
- * 9-2-08: Paul (UCD), changes to clustalOut to make output the same as 1.83
- * added const NAMESWIDTH which is calculated based on the
- * length of the longest sequence name and the upper and lower
- * limits specified in MAXNAMESTODISPLAY and MINNAMESTODISPLAY
- */
-
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include "AlignmentOutput.h"
-
-namespace clustalw
-{
-
-AlignmentOutput::AlignmentOutput()
-: clusSecStructOffset(9),
- clusSequenceOffset(15+1) // MARK made a change here for test case findRangeValues 10seqs
-{
- // NOTE in the old clustal these arrays ended with NULL as a value.
- try
- {
- strongGroup.resize(9);
- strongGroup[0] = string("STA");
- strongGroup[1] = string("NEQK");
- strongGroup[2] = string("NHQK");
- strongGroup[3] = string("NDEQ");
- strongGroup[4] = string("QHRK");
- strongGroup[5] = string("MILV");
- strongGroup[6] = string("MILF");
- strongGroup[7] = string("HY");
- strongGroup[8] = string("FYW");
-
- weakGroup.resize(11);
- weakGroup[0] = string("CSA");
- weakGroup[1] = string("ATV");
- weakGroup[2] = string("SAG");
- weakGroup[3] = string("STNK");
- weakGroup[4] = string("STPA");
- weakGroup[5] = string("SGND");
- weakGroup[6] = string("SNDEQK");
- weakGroup[7] = string("NDEQHK");
- weakGroup[8] = string("NEQHRK");
- weakGroup[9] = string("FVLIM");
- weakGroup[10] = string("HFY");
- }
- catch(const exception& e)
- {
- cerr << "An exception has occured in the contructor of AlignmentOutput.\n"
- << e.what() << "\n";
- exit(1);
- }
-}
-
-/*
- * The function createAlignmentOutput is used to call all the individual functions
- * for the different file types. It is possible to output the alignment in all file
- * types at the same time.
- */
-void AlignmentOutput::createAlignmentOutput(Alignment* alignPtr, int firstSeq, int lastSeq)
-{
- int length;
- int firstRes; // starting sequence range - Ramu
- int lastRes; // ending sequence range
- bool rangeOK;
-
- if((firstSeq <= 0) || (lastSeq < firstSeq))
- {
- utilityObject->error("Cannot produce alignment output."
- " Incorrect call to createAlignmentOutput."
- " firstSeq = %d"
- " lastSeq = %d\n", firstSeq, lastSeq);
- return;
- }
-
- length = 0;
- firstRes = 1;
- lastRes = 0;
- rangeOK = false;
-
- try
- {
- length = alignPtr->getLengthLongestSequence();
- lastRes = length;
-
- if (userParameters->getRangeFromToSet())
- {
- firstRes = userParameters->getRangeFrom();
- lastRes = userParameters->getRangeTo();
- // Check if the numbers are ok.
- if ((firstRes > lastRes) || (firstRes == -1) || (lastRes == -1))
- {
- cerr << "seqrange numbers are not set properly, using default....\n";
- firstRes = 1;
- lastRes = length;
- }
- else
- {
- rangeOK = true;
- }
- }
- if (rangeOK && (lastRes > length))
- {
- lastRes = length;
- cout << "Seqrange " << lastRes << " is more than the " << length
- << " setting it to " << length << " \n";
- }
-
- if (userParameters->getMenuFlag())
- {
- cout << "Consensus length = " << lastRes << " \n";
- }
-
- outputRegion partToOutput;
- partToOutput._firstSeq = firstSeq;
- partToOutput._lastSeq = lastSeq;
- partToOutput._firstRes = firstRes;
- partToOutput._lastRes = lastRes;
-
- if(userParameters->getOutputClustal())
- {
- if(clustalOutFile.get() != 0 && clustalOutFile->is_open())
- {
- clustalOut(alignPtr, partToOutput);
- if(clustalOutFile->is_open())
- {
- clustalOutFile->close();
- }
- utilityObject->info("CLUSTAL-Alignment file created [%s]\n",
- clustalOutName.c_str());
- }
- }
- if(userParameters->getOutputNbrf())
- {
- if(nbrfOutFile.get() != 0 && nbrfOutFile->is_open())
- {
- nbrfOut(alignPtr, partToOutput);
- if(nbrfOutFile->is_open())
- {
- nbrfOutFile->close();
- }
- utilityObject->info("NBRF/PIR-Alignment file created [%s]\n",
- nbrfOutName.c_str());
- }
- }
- if(userParameters->getOutputGCG())
- {
- if(gcgOutFile.get() != 0 && gcgOutFile->is_open())
- {
- gcgOut(alignPtr, partToOutput);
- if(gcgOutFile->is_open())
- {
- gcgOutFile->close();
- }
- utilityObject->info("GCG-Alignment file created [%s]\n",
- gcgOutName.c_str());
- }
- }
- if(userParameters->getOutputPhylip())
- {
- if(phylipOutFile.get() != 0 && phylipOutFile->is_open())
- {
- phylipOut(alignPtr, partToOutput);
- if(phylipOutFile->is_open())
- {
- phylipOutFile->close();
- }
- utilityObject->info("PHYLIP-Alignment file created [%s]\n",
- phylipOutName.c_str());
- }
- }
-
- if(userParameters->getOutputGde())
- {
- if(gdeOutFile.get() != 0 && gdeOutFile->is_open())
- {
- gdeOut(alignPtr, partToOutput);
- if(gdeOutFile->is_open())
- {
- gdeOutFile->close();
- }
- utilityObject->info("GDE-Alignment file created [%s]\n",
- gdeOutName.c_str());
- }
- }
-
- if(userParameters->getOutputNexus())
- {
- if(nexusOutFile.get() != 0 && nexusOutFile->is_open())
- {
- nexusOut(alignPtr, partToOutput);
- if(nexusOutFile->is_open())
- {
- nexusOutFile->close();
- }
- utilityObject->info("NEXUS-Alignment file created [%s]\n",
- nexusOutName.c_str());
- }
- }
-
- if(userParameters->getOutputFasta())
- {
- if(fastaOutFile.get() != 0 && fastaOutFile->is_open())
- {
- fastaOut(alignPtr, partToOutput);
- if(fastaOutFile->is_open())
- {
- fastaOutFile->close();
- }
- utilityObject->info("Fasta-Alignment file created [%s]\n",
- fastaOutName.c_str());
- }
- }
-
- // Output the alignment to the screen if it is required.
- if (userParameters->getShowAlign() && userParameters->getMenuFlag())
- {
- showAlign();
- }
- }
- catch(const exception& e)
- {
- cerr << "An exception has occured in the createAlignmentOutput function.\n"
- << e.what() << "\n";
- exit(1);
- }
-}
-
-bool AlignmentOutput::QTOpenFilesForOutput(AlignmentFileNames fileNames)
-{
- if(!userParameters->getOutputClustal() &&
- !userParameters->getOutputNbrf() && !userParameters->getOutputGCG() &&
- !userParameters->getOutputPhylip() && !userParameters->getOutputGde() &&
- !userParameters->getOutputNexus() && !userParameters->getOutputFasta())
- {
- utilityObject->error("You must select an alignment output format\n");
- return false;
- }
-
- if(fileNames.clustalFile == "" && fileNames.fastaFile == "" && fileNames.gcgFile == ""
- && fileNames.gdeFile == "" && fileNames.nexusFile == "" && fileNames.nrbfFile == ""
- && fileNames.phylipFile == "")
- {
- utilityObject->error("No names for output files. Cannot output alignment.\n");
- return false;
- }
-
- if(fileNames.clustalFile != "")
- {
- clustalOutName = fileNames.clustalFile;
- if(!openExplicitFile(clustalOutFile, clustalOutName))
- {
- return false;
- }
- }
- if(fileNames.fastaFile != "")
- {
- fastaOutName = fileNames.fastaFile;
- if(!openExplicitFile(fastaOutFile, fastaOutName))
- {
- return false;
- }
- }
- if(fileNames.gcgFile != "")
- {
- gcgOutName = fileNames.gcgFile;
- if(!openExplicitFile(gcgOutFile, gcgOutName))
- {
- return false;
- }
- }
- if(fileNames.gdeFile != "")
- {
- gdeOutName = fileNames.gdeFile;
- if(!openExplicitFile(gdeOutFile, gdeOutName))
- {
- return false;
- }
- }
- if(fileNames.nexusFile != "")
- {
- nexusOutName = fileNames.nexusFile;
- if(!openExplicitFile(nexusOutFile, nexusOutName))
- {
- return false;
- }
- }
- if(fileNames.nrbfFile != "")
- {
- nbrfOutName = fileNames.nrbfFile;
- if(!openExplicitFile(nbrfOutFile, nbrfOutName))
- {
- return false;
- }
- }
- if(fileNames.phylipFile != "")
- {
- phylipOutName = fileNames.phylipFile;
- if(!openExplicitFile(phylipOutFile, phylipOutName))
- {
- return false;
- }
- }
- return true;
-}
-/*
- * The function openAlignmentOutput opens a file for output. It returns true if it
- * has been successful and false if it has not been successful
- */
-bool AlignmentOutput::openAlignmentOutput(string path)
-{
- if(!userParameters->getOutputClustal() &&
- !userParameters->getOutputNbrf() && !userParameters->getOutputGCG() &&
- !userParameters->getOutputPhylip() && !userParameters->getOutputGde() &&
- !userParameters->getOutputNexus() && !userParameters->getOutputFasta())
- {
- utilityObject->error("You must select an alignment output format\n");
- return false;
- }
- string _fileNameToOutput = path;
- if(_fileNameToOutput == "")
- {
- /*_fileNameToOutput = userParameters->getSeqName();*/
- /* BUG 166 , file extension, FS, 2009-01-26 */
- utilityObject->getPath(userParameters->getSeqName(), &_fileNameToOutput);
- }
-
- if(userParameters->getOutputClustal())
- {
- if (userParameters->getOutfileName() != "")
- {
- clustalOutName = userParameters->getOutfileName();
- if(!openExplicitFile(clustalOutFile, clustalOutName))
- {
- return false;
- }
- }
- else
- {
- clustalOutName = openOutputFile(clustalOutFile,
- "\nEnter a name for the CLUSTAL output file ",
- _fileNameToOutput, "aln");
-
- if(clustalOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputNbrf())
- {
- if (userParameters->getOutfileName() != "")
- {
- nbrfOutName = userParameters->getOutfileName();
- if(!openExplicitFile(nbrfOutFile, nbrfOutName))
- {
- return false;
- }
- }
- else
- {
- nbrfOutName = openOutputFile(nbrfOutFile,
- "\nEnter a name for the NBRF/PIR output file", _fileNameToOutput,
- "pir");
-
- if(nbrfOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputGCG())
- {
- if (userParameters->getOutfileName() != "")
- {
- gcgOutName = userParameters->getOutfileName();
- if(!openExplicitFile(gcgOutFile, gcgOutName))
- {
- return false;
- }
- }
- else
- {
- gcgOutName = openOutputFile(gcgOutFile,
- "\nEnter a name for the GCG output file", _fileNameToOutput, "msf");
-
- if(gcgOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputPhylip())
- {
- if (userParameters->getOutfileName() != "")
- {
- phylipOutName = userParameters->getOutfileName();
- if(!openExplicitFile(phylipOutFile, phylipOutName))
- {
- return false;
- }
- }
- else
- {
- phylipOutName = openOutputFile(phylipOutFile,
- "\nEnter a name for the PHYLIP output file", _fileNameToOutput, "phy");
-
- if(phylipOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputGde())
- {
- if (userParameters->getOutfileName() != "")
- {
- gdeOutName = userParameters->getOutfileName();
- if(!openExplicitFile(gdeOutFile, gdeOutName))
- {
- return false;
- }
- }
- else
- {
- gdeOutName = openOutputFile(gdeOutFile,
- "\nEnter a name for the GDE output file ",
- _fileNameToOutput, "gde");
-
- if(gdeOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputNexus())
- {
- if (userParameters->getOutfileName() != "")
- {
- nexusOutName = userParameters->getOutfileName();
- if(!openExplicitFile(nexusOutFile, nexusOutName))
- {
- return false;
- }
- }
- else
- {
- nexusOutName = openOutputFile(nexusOutFile,
- "\nEnter a name for the NEXUS output file ",
- _fileNameToOutput, "nxs");
-
- if(nexusOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- if(userParameters->getOutputFasta())
- {
- if (userParameters->getOutfileName() != "")
- {
- fastaOutName = userParameters->getOutfileName();
- if(!openExplicitFile(fastaOutFile, fastaOutName))
- {
- return false;
- }
- }
- else
- {
- fastaOutName = openOutputFile(fastaOutFile,
- "\nEnter a name for the Fasta output file ",
- _fileNameToOutput, "fasta");
-
- if(fastaOutName == "")
- {
- return false; // We have not been successful.
- }
- }
- }
-
- return true;
-
-}
-
-
-/*
- * The function openExplicitFile is used to open a file when we have the name already.
- * It returns true if the file has been opened correctly, false otherwise.
- */
-bool AlignmentOutput::openExplicitFile(auto_ptr<ofstream>& outFile, string fileName)
-{
- if (fileName == "")
- {
- cerr << "Bad output file [" << fileName << "]\n";
- utilityObject->error("Bad output file [%s]\n", fileName.c_str());
- return false;
- }
-
- outFile.reset(new ofstream(fileName.c_str(), ofstream::trunc));
- if(!outFile->is_open())
- {
- utilityObject->error("Cannot open output file [%s]\n", fileName.c_str());
- return false;
- }
- return true;
-}
-
-/*
- * The function openOutputFile is used to open a file when we dont have the name yet.
- * This function returns the name if it has been successful. If it has not been
- * successful it returns a blank string ""
- */
-string AlignmentOutput::openOutputFile(auto_ptr<ofstream>& outFile, string prompt,
- string path, string fileExtension)
-{
- string temp;
- string _fileName; // Will return this name.
- string message;
- _fileName = path + fileExtension;
-
- if(_fileName.compare(userParameters->getSeqName()) == 0)
- {
- cout << "Output file name is the same as input file.\n";
- if (userParameters->getMenuFlag())
- {
- message = "\n\nEnter new name to avoid overwriting [" + _fileName + "]";
-
- utilityObject->getStr(message, temp);
- if(temp != "")
- {
- _fileName = temp;
- }
- }
- }
- else if (userParameters->getMenuFlag())
- {
-
- message = prompt + " [" + _fileName + "]";
- utilityObject->getStr(message, temp);
- if(temp != "")
- {
- _fileName = temp;
- }
- }
-
- outFile.reset(new ofstream(_fileName.c_str(), ofstream::trunc));
- if(!outFile->is_open())
- {
- utilityObject->error("Cannot open output file [%s]\n", _fileName.c_str());
- return "";
- }
- return _fileName;
-}
-
-/*
- * fastaOut: output the alignment in FASTA format
- */
-void AlignmentOutput::fastaOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- char residue;
- int val;
- int i, ii;
- int j, slen;
- int lineLength;
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
- rangeNum rnum;
- cout << "firstres = " << firstRes << " lastres = " << lastRes << "\n";
- try
- {
- const SeqArray* alignment = alignPtr->getSeqArray(); //NOTE june29
- vector<char> sequence;
- sequence.assign(lastRes + 1, '0');
-
- lineLength = PAGEWIDTH - alignPtr->getMaxNames();
-
- lineLength = lineLength - lineLength % 10; // round to a multiple of 10
- if (lineLength > LINELENGTH || lineLength <= 0) // Mark 18-7-07
- {
- lineLength = LINELENGTH;
- }
-
-
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- slen = 0;
- //for(j = firstRes; j < firstRes + lastRes; j++) //- nige
- for(j = firstRes; j <= lastRes; j++) //- nige
- {
- if (j <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '-';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- if (userParameters->getLowercase())
- {
- sequence[j - firstRes] = residue;
- }
- else
- {
- sequence[j - firstRes] = residue;
- }
- slen++;
- }
- (*fastaOutFile) << ">" << nameonly(alignPtr->getName(i));
- if(userParameters->getSeqRange())
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- (*fastaOutFile) << "/" << rnum.start << "-" << rnum.end;
-
- }
- (*fastaOutFile) << "\n";
- for(j = 1; j <= slen; j++)
- {
- (*fastaOutFile) << sequence[j-1];
- if((j % lineLength == 0) || (j == slen))
- {
- (*fastaOutFile) << "\n";
- }
- }
- }
- fastaOutFile->close();
- cout << "FASTA file created!\n";
- }
- catch(const bad_alloc& e)
- {
- fastaOutFile->close();
- cerr << "A bad_alloc exception has occured in the fastaOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- fastaOutFile->close();
- cerr << "An exception has occured in the fastaOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- fastaOutFile->close();
- cerr << "An exception has occured in the fastaOut function.\n";
- exit(1);
- }
-
-}
-
-/*
- * gcgOut: output the alignment in gcg file format
- */
-void AlignmentOutput::gcgOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- char residue;
- int val;
- int i, ii, chunks, block;
- int j, k, pos1, pos2;
- long grandChecksum;
-
- try
- {
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
- rangeNum rnum;
- vector<char> sequence;
- vector<int> allChecks;
- int _maxLength = alignPtr->getMaxAlnLength();
- sequence.assign(_maxLength + 1, '0');
- allChecks.assign(lastSeq + 1, 0);
- const SeqArray* alignment = alignPtr->getSeqArray();
-
- for(i = firstSeq; i <= lastSeq; i++)
- {
- //for(j = firstRes; j <= firstRes + lastRes - 1; j++) //- nige
- for(j = firstRes; j <= lastRes; j++) //- nige
- {
- if (j <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '.';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- sequence[j - firstRes + 1] = residue;
- }
- // pad any short sequences with gaps, to make all sequences the same length
- for(; j <= firstRes + lastRes - 1; j++)
- {
- sequence[j - firstRes + 1] = '.';
- }
- allChecks[i] = SeqGCGCheckSum(&sequence, lastRes);
- }
- int _index;
- grandChecksum = 0;
- for(i = 1; i <= alignPtr->getNumSeqs(); i++)
- {
- _index = alignPtr->getOutputIndex(i - 1);
- grandChecksum += allChecks[_index];
- }
- grandChecksum = grandChecksum % 10000;
-
- (*gcgOutFile) << "PileUp\n\n";
- (*gcgOutFile) << "\n\n MSF:" << setw(5) << lastRes << " Type: ";
-
- if(userParameters->getDNAFlag())
- {
- (*gcgOutFile) << "N";
- }
- else
- {
- (*gcgOutFile) << "P";
- }
-
- (*gcgOutFile) << " Check:" << setw(6) << grandChecksum << " .. \n\n";
- float _seqWeight;
-
- int length = lastRes - firstRes + 1; //- nige
-
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- _seqWeight = (alignPtr->getSeqWeight(i - 1) * 100) / (float) INT_SCALE_FACTOR;
-
- (*gcgOutFile) << " Name: " << alignPtr->getName(i) << " oo Len:"
- //<< setw(5) << lastRes << " Check:" << setw(6) << allChecks[i] //- nige
- << setw(5) << length << " Check:" << setw(6) << allChecks[i]
- << " Weight: " << fixed << setprecision(1) << _seqWeight << "\n";
- }
- (*gcgOutFile) << "\n//\n";
-
- chunks = length / GCG_LINELENGTH; //- nige
- if(length % GCG_LINELENGTH != 0) //- nige
- {
- ++chunks;
- }
-
- for(block = 1; block <= chunks; block++)
- {
- (*gcgOutFile) << "\n\n";
- pos1 = ((block - 1) * GCG_LINELENGTH) + 1;
- pos2 = (length < pos1 + GCG_LINELENGTH - 1)? length : pos1 + GCG_LINELENGTH - 1;//- nige
-
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- if (!userParameters->getSeqRange())
- {
- (*gcgOutFile) << "\n" << setw(alignPtr->getMaxNames() + 5) << left
- << alignPtr->getName(i) << " ";
- }
- else
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- std::stringstream ss;
- std::stringstream ss2;
- string rangeStart;
- string rangeEnd;
- ss << rnum.start;
- ss >> rangeStart;
- ss2 << rnum.end;
- ss2 >> rangeEnd;
- string nameAndRange = nameonly(alignPtr->getName(i)) + "/" + rangeStart;
- nameAndRange += "-" + rangeEnd;
- (*gcgOutFile) << "\n" << setw(alignPtr->getMaxNames() + 15) << left
- << nameAndRange;
- }
- for(j = pos1, k = 1; j <= pos2; j++, k++)
- {
-
- //JULIE -
- //check for int sequences - pad out with '.' characters to end
-
- if (j + firstRes - 1 <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j + firstRes - 1);
- val = (*alignment)[i][j + firstRes - 1]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- residue = '.';
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '.';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- (*gcgOutFile) << residue;
- if(j % 10 == 0)
- {
- (*gcgOutFile) << " ";
- }
- }
- }
- }
- (*gcgOutFile) << "\n\n";
- gcgOutFile->close();
- }
- catch(const bad_alloc& e)
- {
- gcgOutFile->close();
- cerr << "A bad_alloc exception has occured in the gcgOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- gcgOutFile->close();
- cerr << "An exception has occured in the gcgOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- gcgOutFile->close();
- cerr << "An exception has occured in the gcgOut function.\n";
- exit(1);
- }
-}
-
-void AlignmentOutput::phylipOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
- try
- {
- char residue;
- int val;
- int i, ii, chunks, block;
- int j, k, pos1, pos2;
- int nameLen;
- bool warn;
- vector<string> _seqNames;
- const SeqArray* alignment = alignPtr->getSeqArray();
- rangeNum rnum;
-
- // Push on a blank string. This is because the index starts at 1 for the seqs etc..
- _seqNames.push_back("");
-
- nameLen = 0;
-
- for(i = firstSeq; i <= lastSeq; i++)
- {
- _seqNames.push_back(alignPtr->getName(i));
- ii = _seqNames.size();
- if(nameLen < ii)
- {
- nameLen = ii;
- }
- }
- if(nameLen > 10)
- {
- warn = false;
- for(i = 0; i < (int)_seqNames.size() - 1; i++)
- {
- for(j = i + 1; j < (int)_seqNames.size(); j++)
- {
- if (_seqNames[i].compare(_seqNames[j]) == 0)
- {
- warn = true;
- }
- }
- }
- if(warn)
- {
- utilityObject->warning("Truncating sequence names to 10 characters for PHYLIP output.\nNames in the PHYLIP format file are NOT unambiguous.");
- }
- else
- {
- utilityObject->warning("Truncating sequence names to 10 characters for PHYLIP output.");
- }
- }
-
- int length = lastRes - firstRes + 1; //- nige
-
- chunks = length / GCG_LINELENGTH; //- nige
- if(length % GCG_LINELENGTH != 0) //- nige
- {
- ++chunks;
- }
-
- (*phylipOutFile) << setw(6) << alignPtr->getNumSeqs() << " "
- << setw(6) << length; //-nige
-
- for(block = 1; block <= chunks; block++)
- {
- pos1 = ((block - 1) * GCG_LINELENGTH) + 1;
- //pos2 = (lastRes < pos1 + GCG_LINELENGTH - 1)? lastRes : pos1 + GCG_LINELENGTH - 1; //- nige
- pos2 = (length < pos1 + GCG_LINELENGTH - 1)? length : pos1 + GCG_LINELENGTH - 1;
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- if(block == 1)
- {
- if(!userParameters->getSeqRange())
- {
- string name = alignPtr->getName(i);
- (*phylipOutFile) << "\n" << setw(10) << left
- << name.substr(0,10) << " ";
- }
- else
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- std::stringstream ss;
- std::stringstream ss2;
- string rangeStart;
- string rangeEnd;
- ss << rnum.start;
- ss >> rangeStart;
- ss2 << rnum.end;
- ss2 >> rangeEnd;
- string nameAndRange = nameonly(alignPtr->getName(i)) + "/";
- nameAndRange += rangeStart + "-" + rangeEnd;
- (*phylipOutFile) << "\n" << setw(alignPtr->getMaxNames() + 15)
- << left
- << nameAndRange;
- }
- }
- else
- {
- (*phylipOutFile) << "\n ";
- }
- for(j = pos1, k = 1; j <= pos2; j++, k++)
- {
- if (j + firstRes - 1 <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j + firstRes - 1);
- val = (*alignment)[i][j + firstRes - 1]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '-';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- (*phylipOutFile) << residue;
- if(j % 10 == 0)
- {
- (*phylipOutFile) << " ";
- }
- }
- }
- (*phylipOutFile) << "\n";
- }
- }
- catch(const bad_alloc& e)
- {
- phylipOutFile->close();
- cerr << "A bad_alloc exception has occured in the phylipOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- phylipOutFile->close();
- cerr << "An exception has occured in the phylipOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- phylipOutFile->close();
- cerr << "An exception has occured in the phylipOut function.\n";
- exit(1);
- }
-}
-
-void AlignmentOutput::nexusOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
- try
- {
- char residue;
- int val;
- int i, ii, chunks, block;
- int j, k, pos1, pos2;
- const SeqArray* alignment = alignPtr->getSeqArray();
- rangeNum rnum;
-
- int length = lastRes - firstRes + 1; //- nige
-
- chunks = length / GCG_LINELENGTH; //- nige
- if(length % GCG_LINELENGTH != 0) //- nige
- {
- ++chunks;
- }
-
- (*nexusOutFile) << "#NEXUS\n";
- (*nexusOutFile) << "BEGIN DATA;\n";
- (*nexusOutFile) << "dimensions ntax=" << alignPtr->getNumSeqs()
- << " nchar=" << length << ";\n"; //- nige
- (*nexusOutFile) << "format missing=?\n";
- (*nexusOutFile) << "symbols=\"";
-
- for(i = 0; i <= userParameters->getMaxAA(); i++)
- {
- (*nexusOutFile) << userParameters->getAminoAcidCode(i);
- }
- (*nexusOutFile) << "\"\n";
- (*nexusOutFile) << "interleave datatype=";
- bool _dnaFlag = userParameters->getDNAFlag();
- string _type = _dnaFlag ? "DNA " : "PROTEIN ";
- (*nexusOutFile) << _type;
- (*nexusOutFile) << "gap= -;\n";
- (*nexusOutFile) << "\nmatrix";
-
- for(block = 1; block <= chunks; block++)
- {
- pos1 = ((block - 1) * GCG_LINELENGTH) + 1;
- pos2 = (length < pos1 + GCG_LINELENGTH - 1)? length : pos1 + GCG_LINELENGTH - 1; //- nige
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- if (!userParameters->getSeqRange())
- {
- (*nexusOutFile) << "\n" << setw(alignPtr->getMaxNames() + 1) << left
- << alignPtr->getName(i) << " ";
- }
- else
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
-
- std::stringstream ss;
- std::stringstream ss2;
- string rangeStart;
- string rangeEnd;
- ss << rnum.start;
- ss >> rangeStart;
- ss2 << rnum.end;
- ss2 >> rangeEnd;
- string nameAndRange = nameonly(alignPtr->getName(i)) + "/";
- nameAndRange += rangeStart + "-" + rangeEnd;
- (*nexusOutFile) << "\n" << setw(alignPtr->getMaxNames() + 15)
- << left
- << nameAndRange;
- }
- for(j = pos1, k = 1; j <= pos2; j++, k++)
- {
- if (j + firstRes - 1 <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j + firstRes - 1);
- val = (*alignment)[i][j + firstRes - 1]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '-';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- (*nexusOutFile) << residue;
- }
- }
- (*nexusOutFile) << "\n";
- }
- (*nexusOutFile) << ";\nend;\n";
- }
- catch(const bad_alloc& e)
- {
- nexusOutFile->close();
- cerr << "A bad_alloc exception has occured in the nexusOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- nexusOutFile->close();
- cerr << "An exception has occured in the nexusOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- nexusOutFile->close();
- cerr << "An exception has occured in the nexusOut function.\n";
- exit(1);
- }
-}
-
-/*
- * gdeOut: print out the alignment in gde file format.
- */
-void AlignmentOutput::gdeOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
-
- int length = lastRes - firstRes + 1; //- nige
-
- try
- {
- char residue;
- int val;
- vector<char> _ssMask1, _ssMask2, seq;
- int i, ii;
- int j, slen;
- int lineLength;
- const SeqArray* alignment = alignPtr->getSeqArray();
- rangeNum rnum;
-
- seq.assign(alignPtr->getMaxAlnLength() + 1, '0');
-
- // decide the line length for this alignment - maximum is LINELENGTH
- lineLength = PAGEWIDTH - alignPtr->getMaxNames();
- lineLength = lineLength - lineLength % 10; // round to a multiple of 10
- if (lineLength > LINELENGTH || lineLength <= 0)
- {
- lineLength = LINELENGTH;
- }
-
- // If we are using the secondary structure info from profile 1, set it up
- if (userParameters->getStructPenalties1() == SECST &&
- userParameters->getUseSS1() == true)
- {
- int _lengthSeq1 = alignPtr->getSeqLength(1);
- vector<char>* _secStructMask1 = alignPtr->getSecStructMask1();
- _ssMask1.assign(_lengthSeq1 + 10, 0);
-
- for (i = 0;i < _lengthSeq1;i++)
- {
- _ssMask1[i] = _secStructMask1->at(i);
- }
- printSecStructMask(_lengthSeq1, _secStructMask1, &_ssMask1);
- }
-
- // If we are using the secondary structure info from profile 2, set it up
- if (userParameters->getStructPenalties2() == SECST &&
- userParameters->getUseSS2() == true)
- {
- int indexProfile2FirstSeq = alignPtr->getProfile1NumSeqs() + 1;
- int lengthSeqProfile2 = alignPtr->getSeqLength(indexProfile2FirstSeq);
- _ssMask2.assign(lengthSeqProfile2 + 10, 0);
- vector<char>* _secStructMask2 = alignPtr->getSecStructMask2();
-
- for (i=0;i < lengthSeqProfile2;i++)
- {
- _ssMask2[i] = _secStructMask2->at(i);
- }
- printSecStructMask(lengthSeqProfile2, _secStructMask2, &_ssMask2);
- }
-
- bool _dnaFlag = userParameters->getDNAFlag();
- string _prefix = _dnaFlag ? "#" : "%";
-
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
-
- (*gdeOutFile) << _prefix;
- if(!userParameters->getSeqRange())
- {
- (*gdeOutFile) << alignPtr->getName(i) << "\n";
- }
- else
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- (*gdeOutFile) << nameonly(alignPtr->getName(i)) << "/"
- << rnum.start << "-" << rnum.end << "\n";
- }
- slen = 0;
-
- //for(j = firstRes; j < firstRes + lastRes; j++) //- nige
- for(j = firstRes; j <= lastRes; j++) //- nige
- {
- if (j <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '-';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- if (userParameters->getLowercase())
- {
- seq[j - firstRes] = (char)tolower((int)residue);
- }
- else
- {
- seq[j - firstRes] = residue;
- }
- slen++;
- }
- for(j = 1; j <= slen; j++)
- {
- (*gdeOutFile) << seq[j-1];
- if((j % lineLength == 0) || (j == slen))
- {
- (*gdeOutFile) << "\n";
- }
- }
- }
-
- if (userParameters->getOutputStructPenalties() == 0 ||
- userParameters->getOutputStructPenalties() == 2)
- {
- if (userParameters->getStructPenalties1() == SECST &&
- userParameters->getUseSS1() == true)
- {
- (*gdeOutFile) << "\"SS_" << setw(alignPtr->getMaxNames()) << left
- << alignPtr->getSecStructName1() << "\n";
-
- //for(i = firstRes; i < firstRes + lastRes; i++) //- nige
- for(i = firstRes; i <= lastRes; i++) //- nige
- {
- val = _ssMask1[i - 1];
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- seq[i - firstRes] = '-';
- }
- else
- {
- seq[i - firstRes] = val;
- }
- }
-
- for(i = 1; i <= lastRes; i++)
- {
- (*gdeOutFile) << seq[i-1];
- if((i % lineLength == 0) || (i == lastRes))
- {
- (*gdeOutFile) << "\n";
- }
- }
- }
-
- if (userParameters->getStructPenalties2() == SECST &&
- userParameters->getUseSS2() == true)
- {
- (*gdeOutFile) << "\"SS_" << setw(alignPtr->getMaxNames()) << left
- << alignPtr->getSecStructName2() << "\n";
-
- //for(i = firstRes; i < firstRes + lastRes; i++) //- nige
- for(i = firstRes; i <= lastRes; i++) //- nige
- {
- val = _ssMask2[i - 1];
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- seq[i - firstRes] = '-';
- }
- else
- {
- seq[i - firstRes] = val;
- }
- }
-
- //for(i = 1; i <= lastRes; i++) //- nige
- for(i = 1; i <= length; i++) //- nige
- {
- (*gdeOutFile) << seq[i - 1];
- //if((i % lineLength == 0) || (i == lastRes)) //- nige
- if((i % lineLength == 0) || (i == length)) //- nige
- {
- (*gdeOutFile) << "\n";
- }
- }
- }
- }
- if (userParameters->getOutputStructPenalties() == 1 ||
- userParameters->getOutputStructPenalties() == 2)
- {
- if (userParameters->getStructPenalties1() != NONE &&
- userParameters->getUseSS1() == true)
- {
- (*gdeOutFile) << "\"GM_" << setw(alignPtr->getMaxNames()) << left
- << alignPtr->getSecStructName1() << "\n";
-
- vector<char>* _gapPenaltyMask1 = alignPtr->getGapPenaltyMask1();
- //for(i = firstRes; i < firstRes + lastRes; i++) //- nige
- for(i = firstRes; i <= lastRes; i++) //- nige
- {
- val = _gapPenaltyMask1->at(i - 1);
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- seq[i - firstRes] = '-';
- }
- else
- {
- seq[i - firstRes] = val;
- }
- }
-
- //for(i = 1; i <= lastRes; i++) //- nige
- for(i = 1; i <= length; i++) //- nige
- {
- (*gdeOutFile) << seq[i - 1];
- //if((i % lineLength == 0) || (i == lastRes)) //- nige
- if((i % lineLength == 0) || (i == length)) //- nige
- {
- (*gdeOutFile) << "\n";
- }
- }
- }
- if (userParameters->getStructPenalties2() != NONE &&
- userParameters->getUseSS2() == true)
- {
- (*gdeOutFile) << "\"GM_" << setw(alignPtr->getMaxNames()) << left
- << alignPtr->getSecStructName2() << "\n";
-
- vector<char>* _gapPenaltyMask2 = alignPtr->getGapPenaltyMask2();
- //for(i = firstRes; i < firstRes + lastRes; i++) //- nige
- for(i = firstRes; i < length; i++) //- nige
- {
- val = _gapPenaltyMask2->at(i-1);
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- seq[i - firstRes] = '-';
- }
- else
- {
- seq[i - firstRes] = val;
- }
- }
-
- //for(i = 1; i <= lastRes; i++) //- nige
- for(i = 1; i <= length; i++) //- nige
- {
- (*gdeOutFile) << seq[i - 1];
- //if((i % lineLength == 0) || (i == lastRes)) //- nige
- if((i % lineLength == 0) || (i == length)) //- nige
- {
- (*gdeOutFile) << "\n";
- }
- }
- }
- }
- gdeOutFile->close();
- }
- catch(const bad_alloc& e)
- {
- gdeOutFile->close();
- cerr << "A bad_alloc exception has occured in the gdeOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- gdeOutFile->close();
- cerr << "An exception has occured in the gdeOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- gdeOutFile->close();
- cerr << "An exception has occured in the gdeOut function.\n";
- exit(1);
- }
-}
-
-void AlignmentOutput::nbrfOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- char residue;
- int val;
- int i, ii;
- int j, slen;
- int lineLength;
- rangeNum rnum;
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
-
- try
- {
- int _maxLength = alignPtr->getMaxAlnLength();
- vector<char> sequence;
- sequence.assign(_maxLength + 1, '0');
- const SeqArray* alignment = alignPtr->getSeqArray();
- // decide the line length for this alignment - maximum is LINELENGTH
- lineLength = PAGEWIDTH - alignPtr->getMaxNames();
- lineLength = lineLength - lineLength % 10; // round to a multiple of 10
- if (lineLength > LINELENGTH || lineLength <= 0) // Mark, 18-7-07
- {
- lineLength = LINELENGTH;
- }
-
- // Get name prefix. DL if DNA, P1 if protein
- string namePrefix;
- bool _dnaFlag = userParameters->getDNAFlag();
- namePrefix = _dnaFlag ? ">DL;" : ">P1;";
-
- //int length = lastRes - firstRes + 1; //- nige
-
- for(ii = firstSeq; ii <= lastSeq; ii++)
- {
- i = alignPtr->getOutputIndex(ii - 1);
-
- (*nbrfOutFile) << namePrefix;
-
- if (!userParameters->getSeqRange())
- {
- (*nbrfOutFile) << alignPtr->getName(i) << "\n"
- << alignPtr->getTitle(i) << "\n";
- }
- else
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- (*nbrfOutFile) << nameonly(alignPtr->getName(i)) << "/" << rnum.start
- << "-"<< rnum.end << "\n"
- << alignPtr->getTitle(i) << "\n";
- }
- slen = 0;
- //for(j = firstRes; j < firstRes + lastRes; j++) //- nige
- for(j = firstRes; j <= lastRes; j++) //- nige
- {
- // NOTE I changed this here!!!!!
- if (j <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- residue = '-';
- }
- else
- {
- residue = userParameters->getAminoAcidCode(val);
- }
- sequence[j - firstRes] = residue;
- slen++;
- }
- for(j = 1; j <= slen; j++)
- {
- (*nbrfOutFile) << sequence[j - 1];
- if((j % lineLength == 0) || (j == slen))
- {
- (*nbrfOutFile) << "\n";
- }
- }
- (*nbrfOutFile) << "*\n";
- }
- nbrfOutFile->close();
- }
- catch(const bad_alloc& e)
- {
- nbrfOutFile->close();
- cerr << "A bad_alloc exception has occured in the nbrfOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- nbrfOutFile->close();
- cerr << "An exception has occured in the nbrfOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- nbrfOutFile->close();
- cerr << "An exception has occured in the nbrfOut function.\n";
- exit(1);
- }
-}
-
-/*
- * The function clustalOut is used to ouput the Alignment into a clustal format
- * file.
- */
-void AlignmentOutput::clustalOut(Alignment* alignPtr, outputRegion partToOutput)
-{
- int firstRes = partToOutput._firstRes;
- int lastRes = partToOutput._lastRes;
- int firstSeq = partToOutput._firstSeq;
- int lastSeq = partToOutput._lastSeq;
-
- try
- {
- vector<char> seq1;
- vector<int> seqNo;
- vector<int> printSeqNo;
- vector<char> ss_mask1, ss_mask2;
- const SeqArray* alignment = alignPtr->getSeqArray();
- char temp[MAXLINE];
- char c;
- int val;
- int ii, lv1, catident1[NUMRES], catident2[NUMRES], ident, chunks;
- int i, j, k, l;
- int pos, ptr;
- int lineLength;
-
- rangeNum rnum;
- string tmpStr = "";
- string sequenceLine = "";
-
-
- int _numSequences = alignPtr->getNumSeqs();
- // PMcG revert to Clustalw1.83 output format for name width
- const int NAMESWIDTH=std::max(std::min(MAXNAMESTODISPLAY,alignPtr->getMaxNames()) , MINNAMESTODISPLAY);
-
- seqNo.assign(_numSequences + 1, 0);
- printSeqNo.assign(_numSequences + 1, 0);
-
- // check that lastSeq <= _numSequences
- if(lastSeq > _numSequences)
- {
- lastSeq = _numSequences;
- }
-
- for (i = firstSeq; i <= lastSeq; i++)
- {
- printSeqNo[i] = seqNo[i] = 0;
- for(j = 1; j < firstRes; j++)
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- if((val >= 0) || (val <= userParameters->getMaxAA()))
- {
- seqNo[i]++;
- }
- }
- }
-
- seq1.assign(alignPtr->getMaxAlnLength() + 1, 0);
- // Check if we have secondary structure in file 1 and if we want to output it.
- if (userParameters->getStructPenalties1() == SECST &&
- userParameters->getUseSS1() == true)
- {
- int lengthSeq1 = alignPtr->getSeqLength(1);
- ss_mask1.assign(lengthSeq1 + 10, 0);
- vector<char>* _secStructMask1 = alignPtr->getSecStructMask1();
-
- for (i = 0; i < lengthSeq1; i++)
- {
- ss_mask1[i] = _secStructMask1->at(i);
- }
- printSecStructMask(lengthSeq1, _secStructMask1, &ss_mask1);
-
- }
-
- // Check if we have secondary structure in file 2 and if we want to output it.
- if (userParameters->getStructPenalties2() == SECST &&
- userParameters->getUseSS2() == true &&
- userParameters->getProfile2Empty() == false)
- {
- // AW added test getProfile2Empty was meant to fix bug
- // 141, but only prevents segfault, output still faulty.
- // Has to be fixed somewhere else
- int indexProfile2FirstSeq = alignPtr->getProfile1NumSeqs() + 1;
- int lengthSeqProfile2 = alignPtr->getSeqLength(indexProfile2FirstSeq);
- ss_mask2.assign(lengthSeqProfile2 + 10, 0);
- vector<char>* _secStructMask2 = alignPtr->getSecStructMask2();
- for (i = 0; i < lengthSeqProfile2; i++)
- {
- ss_mask2[i] = _secStructMask2->at(i);
- }
- printSecStructMask(lengthSeqProfile2, _secStructMask2, &ss_mask2);
- }
-
- (*clustalOutFile) << "CLUSTAL "<< userParameters->getRevisionLevel()
- << " multiple sequence alignment\n\n";
-
- // decide the line length for this alignment - maximum is LINELENGTH
-
- //PMcG 9-2-2008 make line output same as 1.83
- //lineLength = PAGEWIDTH - alignPtr->getMaxNames();
- lineLength = PAGEWIDTH - NAMESWIDTH;
- lineLength = lineLength - lineLength % 10; // round to a multiple of 10
- if (lineLength > LINELENGTH || lineLength <= 0) // Mark 18-7-07
- {
- lineLength = LINELENGTH;
- }
-
- int length = lastRes - firstRes + 1; //- nige
-
- chunks = length / lineLength; //- nige
- if(length % lineLength != 0) //- nige
- {
- ++chunks;
- }
- //printf("firstRes=%d,lastRes=%d,length=%d,chunks=%d\n",
- // firstRes, lastRes, length, chunks);
-
- // This will loop through each of the blocks.
- for(lv1 = 1; lv1 <= chunks; ++lv1)
- {
- // pos is begining of chunk, ptr is the end of the chunk to be displayed.
- pos = ((lv1 - 1) * lineLength) + 1;
- ptr = (length < pos + lineLength - 1) ? length : pos + lineLength - 1; //- nige
-
- (*clustalOutFile) << "\n";
- int _outStructPenalty = userParameters->getOutputStructPenalties();
-
- string secStructName1 = alignPtr->getSecStructName1(); // Mark 18-7-07
-
- if((int)secStructName1.size() > MAXNAMESTODISPLAY)
- {
- //secStructName1 = secStructName1.substr(0, MAXNAMESTODISPLAY);
- secStructName1 = secStructName1.substr(0, NAMESWIDTH);
- }
- if (_outStructPenalty == 0 || _outStructPenalty == 2)
- {
- if (userParameters->getStructPenalties1() == SECST &&
- userParameters->getUseSS1() == true)
- {
- for(i = pos; i <= ptr; ++i)
- {
- // Check if we can access mask position first
- if((int)ss_mask1.size() > i + firstRes - 2)
- {
- val = ss_mask1[i + firstRes - 2];
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- temp[i - pos] = '-';
- }
- else
- {
- temp[i - pos] = val;
- }
- }
- else
- {
- break;
- }
- }
- temp[i - pos] = EOS;
-
- if(userParameters->getSeqRange())
- {
- //(*clustalOutFile) << "!SS_" << setw(MAXNAMESTODISPLAY +
- (*clustalOutFile) << "!SS_" << setw(NAMESWIDTH +
- clusSecStructOffset)
- << left << secStructName1
- << " " << temp << "\n";
- }
- else
- {
- //(*clustalOutFile) << "!SS_" << setw(MAXNAMESTODISPLAY)
- (*clustalOutFile) << "!SS_" << setw(NAMESWIDTH)
- << left << secStructName1 << " "
- << temp << "\n";
- }
- }
- }
- if (_outStructPenalty == 1 || _outStructPenalty == 2)
- {
- if (userParameters->getStructPenalties1() != NONE &&
- userParameters->getUseSS1() == true)
- {
- vector<char>* _gapPenaltyMask1 = alignPtr->getGapPenaltyMask1();
- for(i = pos; i <= ptr; ++i)
- {
- if((int)_gapPenaltyMask1->size() > i + firstRes - 2)
- {
- val = _gapPenaltyMask1->at(i + firstRes - 2);
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- temp[i - pos] = '-';
- }
- else
- {
- temp[i - pos] = val;
- }
- }
- else
- {
- break;
- }
- }
- temp[i - pos] = EOS;
-
- //(*clustalOutFile) << "!GM_" << setw(MAXNAMESTODISPLAY) << left
- (*clustalOutFile) << "!GM_" << setw(NAMESWIDTH) << left
- << secStructName1 << " "
- << temp << "\n";
- }
- }
-
-
- string secStructName2 = alignPtr->getSecStructName2(); // Mark 18-7-07
- //if((int)secStructName2.size() > MAXNAMESTODISPLAY)
- if((int)secStructName2.size() > NAMESWIDTH)
- {
- //secStructName2 = secStructName2.substr(0, MAXNAMESTODISPLAY);
- secStructName2 = secStructName2.substr(0, NAMESWIDTH);
- }
-
- if (_outStructPenalty == 0 || _outStructPenalty == 2)
- {
- if (userParameters->getStructPenalties2() == SECST &&
- userParameters->getUseSS2() == true)
- {
- for(i = pos; i <= ptr; ++i)
- {
- if((int)ss_mask2.size() > i + firstRes - 2)
- {
- val = ss_mask2[i + firstRes - 2];
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- temp[i - pos] = '-';
- }
- else
- {
- temp[i - pos] = val;
- }
- }
- else
- {
- break;
- }
- }
- temp[i - pos] = EOS;
-
-
- if (userParameters->getSeqRange())
- {
- //(*clustalOutFile) << "!SS_" << setw(MAXNAMESTODISPLAY +
- (*clustalOutFile) << "!SS_" << setw(NAMESWIDTH +
- clusSecStructOffset)
- << left << secStructName2;
- (*clustalOutFile) << " " << temp << "\n";
- }
- else
- {
- //(*clustalOutFile) << "!SS_" << setw(MAXNAMESTODISPLAY)
- (*clustalOutFile) << "!SS_" << setw(NAMESWIDTH)
- << left << secStructName2 << " "
- << temp << "\n";
- }
- }
- }
-
- if (_outStructPenalty == 1 || _outStructPenalty == 2)
- {
- if (userParameters->getStructPenalties2() != NONE &&
- userParameters->getUseSS2() == true)
- {
- vector<char>* _gapPenaltyMask2 = alignPtr->getGapPenaltyMask2();
- for(i = pos; i <= ptr; ++i)
- {
- if((int)_gapPenaltyMask2->size() > i + firstRes - 2)
- {
- val = _gapPenaltyMask2->at(i + firstRes - 2);
- if (val == userParameters->getGapPos1() ||
- val == userParameters->getGapPos2())
- {
- temp[i - pos] = '-';
- }
- else
- {
- temp[i - pos] = val;
- }
- }
- else
- {
- break;
- }
- }
- temp[i - pos] = EOS;
-
- //(*clustalOutFile) << "!GM_" << setw(MAXNAMESTODISPLAY) << left
- (*clustalOutFile) << "!GM_" << setw(NAMESWIDTH) << left
- << secStructName2 << " "
- << temp << "\n";
- }
- }
-
- for(ii = firstSeq; ii <= lastSeq; ++ii)
- {
- i = alignPtr->getOutputIndex(ii - 1);
- printSeqNo[i] = 0;
- for(j = pos; j <= ptr; ++j)
- {
- if (j + firstRes - 1 <= alignPtr->getSeqLength(i))
- {
- //val = alignPtr.getResidue(i, j + firstRes - 1);
- val = (*alignment)[i][j + firstRes - 1]; // NOTE june29
- }
- else
- {
- val = -3;
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- seq1[j] = '-';
- }
- else
- {
- seq1[j] = userParameters->getAminoAcidCode(val);
- seqNo[i]++;
- printSeqNo[i] = 1;
- }
- }
- for(; j <= ptr; ++j)
- {
- seq1[j]='-';
- }
-
- sequenceLine = "";
-
- for(int index = pos; index < ptr + 1; index++)
- {
- sequenceLine += seq1[index];
- }
-
- string seqName = alignPtr->getName(i);
- //if((int)seqName.size() > MAXNAMESTODISPLAY)
- if((int)seqName.size() > NAMESWIDTH)
- {
- //seqName = seqName.substr(0, MAXNAMESTODISPLAY);
- seqName = seqName.substr(0, NAMESWIDTH);
- }
-
- if (!userParameters->getSeqRange())
- {
- // NOTE I made a change here from + 5, to + 6.
- //(*clustalOutFile) << setw(alignPtr->getMaxNames() + 6) << left
- // << alignPtr->getName(i);
-
- //(*clustalOutFile) << setw(MAXNAMESTODISPLAY + 6) << left
- // << seqName;
- // PMcG changed this to revert to behaviour of clustalw1.83
- //
- //(*clustalOutFile) << setw(std::max(std::min(MAXNAMESTODISPLAY,alignPtr->getMaxNames()) , MINNAMESTODISPLAY) + 6)
- (*clustalOutFile) << setw(NAMESWIDTH + 6) << left << seqName;
-
-
- }
- else // Put the sequence range information in!
- {
- findRangeValues(alignPtr, &rnum, firstRes, lastRes, ii);
- string nameOnly = nameonly(seqName);
- std::stringstream ss;
- std::stringstream ss2;
- string rangeStart;
- string rangeEnd;
- ss << rnum.start;
- ss >> rangeStart;
- ss2 << rnum.end;
- ss2 >> rangeEnd;
- tmpStr = nameOnly + "/" + rangeStart + "-" + rangeEnd;
- //(*clustalOutFile) << setw(MAXNAMESTODISPLAY + clusSequenceOffset)
- (*clustalOutFile) << setw(NAMESWIDTH + clusSequenceOffset)
- << left << tmpStr;
- }
-
- (*clustalOutFile) << sequenceLine;
-
- if (userParameters->getClSeqNumbers() && printSeqNo[i])
- {
- (*clustalOutFile) << " " << seqNo[i];
- }
-
- (*clustalOutFile) << "\n";
- }
-
- // Now print out the conservation information!
- for(i = pos; i <= ptr; ++i)
- {
- seq1[i] = ' ';
- ident = 0;
-
- for(j = 1; j <= (int)strongGroup.size(); j++)
- {
- catident1[j - 1] = 0;
- }
- for(j = 1; j <= (int)weakGroup.size(); j++)
- {
- catident2[j - 1] = 0;
- }
-
- for(j = firstSeq; j <= lastSeq; ++j)
- {
- if((i + firstRes - 1 <= alignPtr->getSeqLength(j)) &&
- (i + firstRes - 1 <= alignPtr->getSeqLength(firstSeq)))
- {// NOTE june29
- if(((*alignment)[firstSeq][i + firstRes - 1] >= 0) &&
- ((*alignment)[firstSeq][i + firstRes - 1] <=
- userParameters->getMaxAA()))
- {
- // Count how many are identical to the first sequence
- if((*alignment)[firstSeq][i + firstRes - 1] ==
- (*alignment)[j][i + firstRes - 1])
- {
- ++ident;
- }
- // Count how many are in the same category.
- for(k = 1; k <= (int)strongGroup.size(); k++)
- {
- for(l = 0; (c = strongGroup[k - 1][l]); l++)
- {
- int resCode = (*alignment)[j][i + firstRes - 1];
- if(resCode <= userParameters->getMaxAA() + 1)
- {
- if (userParameters->getAminoAcidCode(resCode) == c)
- {
- catident1[k - 1]++;
- break;
- }
- }
- }
- }
- for(k = 1; k <= (int)weakGroup.size(); k++)
- {
- for(l = 0; (c = weakGroup[k - 1][l]); l++)
- {//NOTE june29
- int resCode = (*alignment)[j][i + firstRes - 1];
- if(resCode <= userParameters->getMaxAA() + 1)
- {
- if (userParameters->getAminoAcidCode(resCode) == c)
- {
- catident2[k - 1]++;
- break;
- }
- }
- }
- }
- }
- }
- }
-
- // Now do the conservation part for each block.
- if(ident == lastSeq - firstSeq + 1)
- {
- seq1[i] = '*'; // All residues the same!
- }
- else if (!userParameters->getDNAFlag())
- {
- for(k = 1; k <= (int)strongGroup.size(); k++)
- {
- if (catident1[k - 1] == lastSeq - firstSeq + 1)
- {
- seq1[i] = ':'; // All residues member of the same category
- break;
- }
- }
- if(seq1[i] == ' ')
- {
- for(k = 1; k <= (int)weakGroup.size(); k++)
- {
- if (catident2[k - 1] == lastSeq - firstSeq + 1)
- {
- seq1[i] = '.'; // All residues member of the same category
- break;
- }
- }
- }
- }
- }
-
- sequenceLine = "";
- for(int index = pos; index < ptr + 1; index++)
- {
- sequenceLine += seq1[index];
- }
-
- //for(k = 0; k < MAXNAMESTODISPLAY + 6; k++)
- for(k = 0; k < NAMESWIDTH + 6; k++)
- {
- (*clustalOutFile) << " ";
- }
- if(userParameters->getSeqRange())
- {
-
- (*clustalOutFile) << " ";
- }
-
- (*clustalOutFile) << sequenceLine << "\n";
- }
- clustalOutFile->close();
- }
- catch(const bad_alloc& e)
- {
- clustalOutFile->close();
- cerr << "A bad_alloc exception has occured in the clustalOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(VectorOutOfRange e)
- {
- clustalOutFile->close();
- cerr << "An exception has occured in the clustalOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(exception e)
- {
- clustalOutFile->close();
- cerr << "An exception has occured in the clustalOut function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- clustalOutFile->close();
- cerr << "An exception has occured in the clustalOut function.\n";
- exit(1);
- }
-}
-
-/*
- * The function printSecStructMask takes in the 2 mask vectors. It makes changes to
- * mask and puts the result in structMask. This is the one that will be printed out.
- * NOTE I have had to use (*structMask)[i] syntax. This is not good. But there is no
- * way to do random access setting in vectors otherwise.
- */
-void AlignmentOutput::printSecStructMask(int prfLength, vector<char>* mask,
- vector<char>* structMask)
-{
- int i, j;
-// calculate the gap penalty mask from the secondary structures
-
- i = 0;
- try
- {
- while (i < prfLength)
- {
- if (tolower(mask->at(i)) == 'a' || mask->at(i) == '$')
- {
- for (j = 0; j < userParameters->getHelixEndMinus(); j++)
- {
- if (i + j >= prfLength || (tolower(mask->at(i+j)) != 'a'
- && mask->at(i+j) != '$'))
- {
- break;
- }
- (*structMask)[i+j] = 'a';
- }
- i += j;
- while (tolower(mask->at(i)) == 'a' || mask->at(i) == '$')
- {
- if (i >= prfLength)
- break;
- if (mask->at(i) == '$')
- {
- (*structMask)[i] = 'A';
- i++;
- break;
- }
- else
- (*structMask)[i] = (*mask)[i];
- i++;
- }
- for (j = 0; j < userParameters->getHelixEndMinus(); j++)
- {
- if ((i-j-1>=0) && (tolower(mask->at(i-j-1)) == 'a'
- || mask->at(i-j-1) == '$'))
- {
- (*structMask)[i - j - 1] = 'a';
- }
- }
- }
- else if (tolower(mask->at(i)) == 'b' || mask->at(i) == '%')
- {
- for (j = 0; j < userParameters->getStrandEndMinus(); j++)
- {
- if (i + j >= prfLength || (tolower(mask->at(i+j)) != 'b'
- && mask->at(i+j) != '%'))
- {
- break;
- }
- (*structMask)[i+j] = 'b';
- }
- i += j;
- while (tolower(mask->at(i)) == 'b' || mask->at(i) == '%')
- {
- if (i >= prfLength)
- break;
- if (mask->at(i) == '%')
- {
- (*structMask)[i] = 'B';
- i++;
- break;
- }
- else
- {
- (*structMask)[i] = mask->at(i);
- }
- i++;
- }
- for (j = 0; j < userParameters->getStrandEndMinus(); j++)
- {
- if ((i - j - 1 >= 0) && (tolower(mask->at(i - j - 1)) == 'b' ||
- mask->at(i-j-1) == '%'))
- {
- (*structMask)[i-j-1] = 'b';
- }
- }
- }
- else
- {
- i++;
- }
- }
- }
- catch(const exception& e)
- {
- // Catch all the exceptions
- cerr << "There has been an exception in the function printSecStructMask\n"
- << e.what() << "\n";
- exit(1);
- }
-}
-
-/*
- * The function findRangeValues is used to find the range to be printed out.
- *
- *
- */
-void AlignmentOutput::findRangeValues(Alignment* alignPtr, rangeNum *rnum, int firstRes,
- int len, int firstSeq)
-{
- assert(rnum);
- assert(firstRes > 0);
- assert(len > 0);
- assert(firstSeq > 0);
-
- try
- {
- int val;
- int i, ii;
- int j, slen;
-
- char tmpName[FILENAMELEN + 15];
- int iStart = 0;
- int iEnd = 0; // to print sequence start-end with names
- int found = 0;
- int ngaps = 0;
- int tmpStart = 0;
- int tmpEnd = 0;
- int ntermgaps = 0;
- int pregaps = 0;
- int tmpk = 0;
- int isRange = 0;
- int formula = 0;
- const SeqArray* alignment = alignPtr->getSeqArray();
-
- tmpName[0] = '\0';
- slen = 0;
-
- ii = firstSeq ;
- i = alignPtr->getOutputIndex(ii - 1); // NOTE Same as elsewhere!
- string name = alignPtr->getName(i);
-
- if( (sscanf(name.c_str(),"%[^/]/%d-%d",tmpName, &tmpStart, &tmpEnd) == 3))
- {
- isRange = 1;
- }
-
- for(tmpk = 1; tmpk < firstRes; tmpk++)
- {
- // do this irrespective of above sscanf
- //val = alignPtr.getResidue(i, tmpk);
- val = (*alignment)[i][tmpk]; // NOTE june29
- if ((val < 0) || (val > userParameters->getMaxAA()))
- {
- //it is gap
- pregaps++;
- }
- }
-
- for(j = firstRes; j < firstRes + len; j++)
- {
- if(j > alignPtr->getSeqLength(i))
- {
- val = -3; // Cant get it so break out.
- }
- else
- {
- //val = alignPtr.getResidue(i, j);
- val = (*alignment)[i][j]; // NOTE june29
- }
- if((val == -3) || (val == 253))
- {
- break;
- }
- else if((val < 0) || (val > userParameters->getMaxAA()))
- {
- ngaps++;
- }
- else
- {
- found = j;
- }
- if ( found && (iStart == 0) )
- {
- iStart = found;
- ntermgaps = ngaps;
- }
- slen++;
- }
- if(userParameters->getSeqRange())
- {
- cout << "Name : " << alignPtr->getName(i) << " "
- << "\n firstRes = "<< firstRes << " "
- << " len = " << len << " "
- << "\n iStart = " << iStart << " "
- << "\n tmpStart = " << tmpStart << " "
- << "\n ngaps = " << ngaps << " "
- << "\n pregaps = " << pregaps << " ";
- if (!isRange)
- {
- formula = iStart - pregaps;
- }
- else
- {
- formula = iStart - pregaps + ( tmpStart == 1 ? 0: tmpStart-1) ;
- }
-
- cout << "\n\nsuggestion iStart - pregaps + tmpStart - ntermgaps = "
- << iStart << " - " << pregaps << " + " << tmpStart << " - " << ntermgaps
- << " formula " << formula << " ";
- }
- else
- {
- cerr << "\n no range found .... strange, iStart = " << iStart;
- formula = 1;
- }
- if (pregaps == firstRes - 1) // all gaps - now the conditions........
- {
- formula = tmpStart ; // keep the previous start...
- }
- formula = (formula <= 0) ? 1: formula;
- if (pregaps == 0 && tmpStart == 0)
- {
- formula = firstRes;
- }
- iEnd = formula + len - ngaps -1;
- rnum->start = formula;
- rnum->end = iEnd;
- cout << "\n check... " << alignPtr->getName(i) << " " << rnum->start
- << " - " << rnum->end;
- cout << " Done checking.........";
- }
- catch(VectorOutOfRange e)
- {
- cerr << "An exception has occured in the findRangeValues function.\n"
- << e.what() << "\n";
- exit(1);
- }
- catch(...)
- {
- cerr << "An exception has occured in findRangeValues function\n";
- exit(1);
- }
-}
-
-string AlignmentOutput::nameonly(string s)
-{
- string tmp;
- try
- {
- int i = 0;
-
- while (i < (int)s.size())
- {
- if(s.at(i) != '/')
- {
- tmp += s.at(i);
- i++;
- }
- else
- {
- break;
- }
- }
- }
- catch(const exception& e)
- {
- cerr << "An exception has occured in the function nameonly\n"
- << e.what();
- exit(1);
- }
- return tmp;
-}
-
-int AlignmentOutput::SeqGCGCheckSum(vector<char>* seq, int length)
-{
- int i;
- long check;
- int seqResIndex = 1;
- //int seqLength = seq->size();
-
- for (i = 0, check = 0; i < length; i++, seqResIndex++)
- {
- char _toUpperCase = (*seq)[seqResIndex];
- check += ((i % 57) + 1) * toupper(_toUpperCase);
- }
-
- return (check % 10000);
-}
-
-void AlignmentOutput::showAlign()
-{
- //FILE *file;
- int numLines;
- char temp[MAXLINE + 1];
- temp[0] = '\0';
- string fileName;
- string answer;
-
- if(userParameters->getOutputClustal())
- {
- fileName = clustalOutName;
- }
- else if(userParameters->getOutputNbrf())
- {
- fileName = nbrfOutName;
- }
- else if(userParameters->getOutputGCG())
- {
- fileName = gcgOutName;
- }
- else if(userParameters->getOutputPhylip())
- {
- fileName = phylipOutName;
- }
- else if(userParameters->getOutputGde())
- {
- fileName = gdeOutName;
- }
- else if(userParameters->getOutputNexus())
- {
- fileName = nexusOutName;
- }
- else if(userParameters->getOutputFasta())
- {
- fileName = fastaOutName;
- }
- else
- {
- return; // No file output type!
- }
-
- ifstream _fileIn;
- _fileIn.open(fileName.c_str(), ios::in);
- _fileIn.seekg(0, std::ios::beg); // start at the beginning
-
- cout << "\n\n";
- numLines = 0;
-
- while(_fileIn.getline(temp, MAXLINE + 1, '\n'))
- {
- //fputs(temp,stdout);
- cout << temp << "\n";
- ++numLines;
- if(numLines >= PAGE_LEN)
- {
- cout << "\n";
- utilityObject->getStr(string("Press [RETURN] to continue or X to stop"), answer);
- if(toupper(answer[0]) == 'X')
- {
- _fileIn.close();
- return;
- }
- else
- {
- numLines = 0;
- }
- }
- }
- _fileIn.close();
- cout << "\n";
- utilityObject->getStr(string("Press [RETURN] to continue"), answer);
-}
-
-}
-
git://source.jalview.org / jabaws.git / blobdiff