C++ HMM::AddContextSpecificPseudocounts方法代码示例

// Read input file (HMM, HHM, or alignment format), and add pseudocounts etc.
void ReadAndPrepare(char* infile, HMM& q, Alignment* qali=NULL)
{
    char path[NAMELEN];

    // Open query file and determine file type
    char line[LINELEN]=""; // input line
    FILE* inf=NULL;
    if (strcmp(infile,"stdin"))
    {
        inf = fopen(infile, "r");
        if (!inf) OpenFileError(infile);
        Pathname(path,infile);
    }
    else
    {
        inf = stdin;
        if (v>=2) printf("Reading HMM / multiple alignment from standard input ...\n(To get a help list instead, quit and type %s -h.)\n",program_name);
        *path='\0';
    }

    fgetline(line,LINELEN-1,inf);

    // Is it an hhm file?
    if (!strncmp(line,"NAME",4) || !strncmp(line,"HH",2))
    {
        if (v>=2) cout<<"Query file is in HHM format\n";

        // Rewind to beginning of line and read query hhm file
        rewind(inf);
        q.Read(inf,path);
        if (v>=2 && q.Neff_HMM>11.0)
            fprintf(stderr,"WARNING: HMM %s looks too diverse (Neff=%.1f>11). Better check the underlying alignment... \n",q.name,q.Neff_HMM);

        // Add transition pseudocounts to query -> q.p[i][a]
        q.AddTransitionPseudocounts();

        if (!*par.clusterfile) { //compute context-specific pseudocounts?
	  // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
	  q.PreparePseudocounts();
	  // Add amino acid pseudocounts to query:  q.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
	  q.AddAminoAcidPseudocounts(par.pcm, par.pca, par.pcb, par.pcc);;
        } else {
	  // Add context specific pseudocount to query
	  q.AddContextSpecificPseudocounts(par.pcm);
        }
        
        q.CalculateAminoAcidBackground();
    }

    // ... or is it an a2m/a3m alignment file
    else if (line[0]=='#' || line[0]=='>')
    {
        Alignment* pali;
        if (qali==NULL) pali=new(Alignment); else pali=qali;
        if (par.calibrate) {
            printf("\nError in %s: only HHM files can be calibrated.\n",program_name);
            printf("Build an HHM file from your alignment with 'hhmake -i %s' and rerun hhsearch with the hhm file\n\n",infile);
            exit(1);
        }

        if (v>=2 && strcmp(infile,"stdin")) cout<<infile<<" is in A2M, A3M or FASTA format\n";

        // Read alignment from infile into matrix X[k][l] as ASCII (and supply first line as extra argument)
        pali->Read(inf,infile,line);

        // Convert ASCII to int (0-20),throw out all insert states, record their number in I[k][i]
        // and store marked sequences in name[k] and seq[k]
        pali->Compress(infile);

        // Sort out the nseqdis most dissimilar sequences for display in the output alignments
        pali->FilterForDisplay(par.max_seqid,par.coverage,par.qid,par.qsc,par.nseqdis);

        // Remove sequences with seq. identity larger than seqid percent (remove the shorter of two)
        pali->N_filtered = pali->Filter(par.max_seqid,par.coverage,par.qid,par.qsc,par.Ndiff);

 	if (par.Neff>=0.999) 
	  pali->FilterNeff();

	// Calculate pos-specific weights, AA frequencies and transitions -> f[i][a], tr[i][a]
        pali->FrequenciesAndTransitions(q);
        if (v>=2 && q.Neff_HMM>11.0)
            fprintf(stderr,"WARNING: alignment %s looks too diverse (Neff=%.1f>11). Better check it with an alignment viewer... \n",q.name,q.Neff_HMM);

        // Add transition pseudocounts to query -> p[i][a]
        q.AddTransitionPseudocounts();

        if (!*par.clusterfile) { //compute context-specific pseudocounts?
	  // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
	  q.PreparePseudocounts();
	  // Add amino acid pseudocounts to query:  p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
	  q.AddAminoAcidPseudocounts(par.pcm, par.pca, par.pcb, par.pcc);
        } else {
	  // Add context specific pseudocount to query
	  q.AddContextSpecificPseudocounts(par.pcm);
        }

        q.CalculateAminoAcidBackground();

        if (qali==NULL) delete(pali);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (!*par.outfile && !*par.alnfile) {
        RemoveExtension(par.outfile, par.infile);
        strcat(par.outfile, ".seq");
    }

    // Prepare CS pseudocounts lib
    if (!par.nocontxt && *par.clusterfile) {
        InitializePseudocountsEngine(par, context_lib, crf, pc_hhm_context_engine,
                                     pc_hhm_context_mode, pc_prefilter_context_engine,
                                     pc_prefilter_context_mode);
    }

    // Set substitution matrix; adjust to query aa distribution if par.pcm==3
    SetSubstitutionMatrix(par.matrix, pb, P, R, S, Sim);

    // Read input file (HMM, HHM, or alignment format), and add pseudocounts etc.
    char input_format = 0;
    ReadQueryFile(par, par.infile, input_format, par.wg, q, qali, pb, S, Sim);

    // Same code as in PrepareQueryHMM(par.infile,input_format,q,qali), except that we add SS prediction
    // Add Pseudocounts, if no HMMER input
    if (input_format == 0) {
        // Transform transition freqs to lin space if not already done
        q->AddTransitionPseudocounts(par.gapd, par.gape, par.gapf, par.gapg,
                                     par.gaph, par.gapi, par.gapb, par.gapb);

        // Comput substitution matrix pseudocounts
        if (par.nocontxt) {
            // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
            q->PreparePseudocounts(R);
            // Add amino acid pseudocounts to query: p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
            q->AddAminoAcidPseudocounts(par.pc_hhm_nocontext_mode,
                                        par.pc_hhm_nocontext_a, par.pc_hhm_nocontext_b,
                                        par.pc_hhm_nocontext_c);
        }
        else {
            // Add full context specific pseudocounts to query
            q->AddContextSpecificPseudocounts(pc_hhm_context_engine,
                                              pc_hhm_context_mode);
        }
    }
    else {
        q->AddAminoAcidPseudocounts(0, par.pc_hhm_nocontext_a,
                                    par.pc_hhm_nocontext_b, par.pc_hhm_nocontext_c);
    }

    q->CalculateAminoAcidBackground(pb);

    if (par.columnscore == 5 && !q->divided_by_local_bg_freqs)
        q->DivideBySqrtOfLocalBackgroundFreqs(
            par.half_window_size_local_aa_bg_freqs, pb);

    // Write consensus sequence to sequence file
    // Consensus sequence is calculated in hhalignment.C, Alignment::FrequenciesAndTransitions()
    if (*par.outfile) {
        FILE* outf = NULL;
        if (strcmp(par.outfile, "stdout")) {
            outf = fopen(par.outfile, "a");
            if (!outf)
                OpenFileError(par.outfile, __FILE__, __LINE__, __func__);
        }
        else
            outf = stdout;
        // OLD
        //// ">name_consensus" -> ">name consensus"
        //strsubst(q->sname[q->nfirst],"_consensus"," consensus");
        //fprintf(outf,">%s\n%s\n",q->sname[q->nfirst],q->seq[q->nfirst]+1);
        // NEW (long header needed for NR30cons database)
        fprintf(outf, ">%s\n%s\n", q->longname, q->seq[q->nfirst] + 1);
        fclose(outf);
    }

    // Print A3M/A2M/FASTA output alignment
    if (*par.alnfile) {
        HalfAlignment qa;
        int n = imin(q->n_display,
                     par.nseqdis + (q->nss_dssp >= 0) + (q->nss_pred >= 0)
                     + (q->nss_conf >= 0) + (q->ncons >= 0));
        qa.Set(q->name, q->seq, q->sname, n, q->L, q->nss_dssp, q->nss_pred,
               q->nss_conf, q->nsa_dssp, q->ncons);

        if (par.outformat == 1)
            qa.BuildFASTA();
        else if (par.outformat == 2)
            qa.BuildA2M();
        else if (par.outformat == 3)
            qa.BuildA3M();
        if (qali->readCommentLine)
            qa.Print(par.alnfile, par.append, qali->longname); // print alignment to outfile
        else
            qa.Print(par.alnfile, par.append);   // print alignment to outfile
    }

    delete qali;
    delete q;

    DeletePseudocountsEngine(context_lib, crf, pc_hhm_context_engine,
                             pc_hhm_context_mode, pc_prefilter_context_engine,
                             pc_prefilter_context_mode);
}