nwa.cpp

ncbi源码

    if ( !x_ReadFastaFile(args["seq1"].AsString(), &seqname1, &v1)) {
        NCBI_THROW(CAppNWAException,
                   eCannotReadFile,
                   "Cannot read file " + args["seq1"].AsString());
    }
    if ( !x_ReadFastaFile(args["seq2"].AsString(), &seqname2, &v2)) {
        NCBI_THROW(CAppNWAException,
                   eCannotReadFile,
                   "Cannot read file" + args["seq2"].AsString());
    }

    // determine sequence/score matrix type
    const SNCBIPackedScoreMatrix* psm = 
      (args["matrix"].AsString() == "blosum62")?
      &NCBISM_Blosum62:
      0;

    auto_ptr<CNWAligner> aligner (
        bMrna2Dna? 
        new SPLALIGNER (&v1[0], v1.size(), &v2[0], v2.size()):
        (bMM?
         new CMMAligner (&v1[0], v1.size(), &v2[0], v2.size(), psm):
         new CNWAligner (&v1[0], v1.size(), &v2[0], v2.size(), psm))
        );

    aligner->SetWm  (args["Wm"]. AsInteger());
    aligner->SetWms (args["Wms"].AsInteger());
    aligner->SetWg  (args["Wg"]. AsInteger());
    aligner->SetWs  (args["Ws"]. AsInteger());

    if( bMrna2Dna ) {
        SPLALIGNER *aligner_mrna2dna = 
            static_cast<SPLALIGNER*> (aligner.get());

        aligner_mrna2dna->SetWi (0, args["Wi0"]. AsInteger());
        aligner_mrna2dna->SetWi (1, args["Wi1"]. AsInteger());
        aligner_mrna2dna->SetWi (2, args["Wi2"]. AsInteger());
        //        aligner_mrna2dna->SetWi (3, args["Wi3"]. AsInteger());

        aligner_mrna2dna->SetIntronMinSize(args["IntronMinSize"]. AsInteger());

        if(bGuides) {
            aligner_mrna2dna->MakePattern(args["pattern"].AsInteger());
        }
    }

    if(bMT && bMM) {
        CMMAligner* pmma = static_cast<CMMAligner*> (aligner.get());
        pmma -> EnableMultipleThreads();
    }
    
    auto_ptr<ofstream> pofs1 (0);
    auto_ptr<ofstream> pofs2 (0);
    auto_ptr<ofstream> pofsAsn (0);
    auto_ptr<ofstream> pofsFastA (0);
    auto_ptr<ofstream> pofsExons (0);

    if(output_type1) {
        pofs1.reset(open_ofstream (args["o1"].AsString()).release());
    }

    if(output_type2) {
        pofs2.reset(open_ofstream (args["o2"].AsString()).release());
    }

    if(output_asn) {
        pofsAsn.reset(open_ofstream (args["oasn"].AsString()).release());
    }

    if(output_fasta) {
        pofsFastA.reset(open_ofstream (args["ofasta"].AsString()).release());
    }

    if(output_exons) {
        pofsExons.reset(open_ofstream (args["oexons"].AsString()).release());
    }

    {{  // setup end penalties
        string ends = args["esf"].AsString();
        bool L1 = ends[0] == 'z';
        bool R1 = ends[1] == 'z';
        bool L2 = ends[2] == 'z';
        bool R2 = ends[3] == 'z';
        aligner->SetEndSpaceFree(L1, R1, L2, R2);
    }}

    int score = aligner->Run();
    cerr << "Score = " << score << endl;

    CNWFormatter formatter (*aligner);
    formatter.SetSeqIds(seqname1, seqname2);

    const size_t line_width = 50;
    string s;
    if(pofs1.get()) {
        formatter.AsText(&s, CNWFormatter::eFormatType1, line_width);
        *pofs1 << s;
    }

    if(pofs2.get()) {
        formatter.AsText(&s, CNWFormatter::eFormatType2, line_width);
        *pofs2 << s;
    }

    if(pofsAsn.get()) {
        formatter.AsText(&s, CNWFormatter::eFormatAsn, line_width);
        *pofsAsn << s;
    }

    if(pofsFastA.get()) {
        formatter.AsText(&s, CNWFormatter::eFormatFastA, line_width);
        *pofsFastA << s;
    }

    if(pofsExons.get()) {
        formatter.AsText(&s, CNWFormatter::eFormatExonTableEx, line_width);
        *pofsExons << s;
    }
    
    if(!output_type1 && !output_type2
       && !output_asn && !output_fasta
       && !output_exons) {
        formatter.AsText(&s, CNWFormatter::eFormatType2, line_width);
        cout << s;
    }
}


void CAppNWA::Exit()
{
    return;
}


bool CAppNWA::x_ReadFastaFile (const string& filename, string* seqname,
                               vector<char>* sequence) const
{
    vector<char>& vOut = *sequence;
    vOut.clear();

    ifstream ifs(filename.c_str());

    // read sequence's name
    string str;
    getline(ifs, str);
    if(!ifs)
        return false;
    istrstream iss (str.c_str());

    char c;
    iss >> c >> *seqname;
    if(!iss)
        return false;

    // read the sequence
    while ( ifs ) {
        string s;
        ifs >> s;
        NStr::ToUpper(s);
        copy(s.begin(), s.end(), back_inserter(vOut));
    }

    return true;
}

END_NCBI_SCOPE

/*
 * ===========================================================================
 * $Log: nwa.cpp,v $
 * Revision 1000.2  2004/06/01 18:05:08  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.27
 *
 * Revision 1.27  2004/05/21 21:41:02  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.26  2004/04/30 13:01:33  kuznets
 * throw -> THROWS
 *
 * Revision 1.25  2003/11/07 18:30:17  kapustin
 * mRna2Dna --> Spliced
 *
 * Revision 1.24  2003/09/30 19:50:28  kapustin
 * Adjust for standard score matrix interface
 *
 * Revision 1.23  2003/09/10 19:11:50  kapustin
 * Add eNotSupported exception for multithreading availability checking
 *
 * Revision 1.22  2003/09/02 22:38:52  kapustin
 * Adjust for the library's changes
 *
 * Revision 1.21  2003/06/17 17:20:44  kapustin
 * CNWAlignerException -> CAlgoAlignException
 *
 * Revision 1.20  2003/06/17 14:51:04  dicuccio
 * Fixed after algo/ rearragnement
 *
 * Revision 1.19  2003/05/23 18:28:27  kapustin
 * Introduce new (generic) splice type
 *
 * Revision 1.18  2003/05/06 20:27:30  kapustin
 * Specify guide size in command line argument
 *
 * Revision 1.17  2003/04/14 19:00:55  kapustin
 * Add guide creation facility.  x_Run() -> x_Align()
 *
 * Revision 1.16  2003/04/02 20:53:31  kapustin
 * Add exon table output format
 *
 * Revision 1.15  2003/03/25 22:06:02  kapustin
 * Support non-canonical splice signals
 *
 * Revision 1.14  2003/03/18 15:14:54  kapustin
 * Allow separate free end gap specification
 *
 * Revision 1.13  2003/03/17 15:32:28  kapustin
 * Enabled end-space free alignments for all currently supported methods
 *
 * Revision 1.12  2003/03/07 13:52:57  kapustin
 * Add a temporary check that -mm is not used with -esf
 *
 * Revision 1.11  2003/03/05 20:13:53  kapustin
 * Simplify FormatAsText(). Fix FormatAsSeqAlign(). Convert sequence alphabets
 * to capitals
 *
 * Revision 1.10  2003/02/11 16:06:55  kapustin
 * Add end-space free alignment support
 *
 * Revision 1.9  2003/01/28 12:46:27  kapustin
 * Format() --> FormatAsText(). Fix the flag spelling forcing ASN output.
 *
 * Revision 1.8  2003/01/24 19:43:03  ucko
 * Change auto_ptr assignment to use release and reset rather than =,
 * which not all compilers support.
 *
 * Revision 1.7  2003/01/24 16:49:59  kapustin
 * Support different output formats
 *
 * Revision 1.6  2003/01/21 16:34:22  kapustin
 * 
 *
 * Revision 1.5  2003/01/21 12:42:02  kapustin
 * Add mm parameter
 *
 * Revision 1.4  2003/01/08 15:58:32  kapustin
 * Read offset parameter from fasta reading routine
 *
 * Revision 1.3  2002/12/17 21:50:24  kapustin
 * Remove unnecesary seq type parameter from the mrna2dna constructor
 *
 * Revision 1.2  2002/12/12 17:59:30  kapustin
 * Enable spliced alignments
 *
 * Revision 1.1  2002/12/06 17:44:25  ivanov
 * Initial revision
 *
 * ===========================================================================
 */