C++ Alignment类代码示例

bool ExtendJobMapperWithLinks::AlTestEdgeProximity::passes(vector<Alignment*>* matches, ScoredSeq* query){
  bool foundPass = false;
  // this allows the query to overhang the window to a degree acceptable by the percent ID,
  // while also including a -1 that compensates for the zero-indexing of the final aligned pos
  long effQsizeM1 = long(float(query->size()) * _fractId) - 1;
  long aLastPos = query->size() - 1;

  vector<Alignment*>::iterator alIt = matches->begin();
  vector<Alignment*>::iterator alEnd = matches->end();
  while (alIt != alEnd){
    Alignment* al = *alIt;
    long bLastPos;
    if ( al->isLinked(aLastPos, query) ){ bLastPos = al->getLinkage(aLastPos, query); }
    else { bLastPos = al->gapPairedAfter(aLastPos, query); }

    // test if the full sequence of the query would fit in the window given where the
    // alignment ends (it is possible that gaps push the 5p end of the query beyond the
    // window, but i don't want that to penalize the match).
    foundPass = _edgeProximity >= al->seqB()->size() - bLastPos + effQsizeM1;

    // if the test passed, skip to the end
    if (foundPass){ alIt = alEnd; }
    else { ++alIt; }
  }
  return foundPass;
}

Alignment Sampler::mutate(const Alignment& aln,
                          double base_error,
                          double indel_error) {

    if (base_error == 0 && indel_error == 0) return aln;

    string bases = "ATGC";
    uniform_real_distribution<double> rprob(0, 1);
    uniform_int_distribution<int> rbase(0, 3);

    Alignment mutaln;
    for (size_t i = 0; i < aln.path().mapping_size(); ++i) {
        auto& orig_mapping = aln.path().mapping(i);
        Mapping* new_mapping = mutaln.mutable_path()->add_mapping();
        *new_mapping->mutable_position() = orig_mapping.position();
        // for each edit in the mapping
        for (size_t j = 0; j < orig_mapping.edit_size(); ++j) {
            auto& orig_edit = orig_mapping.edit(j);
            auto new_edits = mutate_edit(orig_edit, make_pos_t(orig_mapping.position()),
                                         base_error, indel_error,
                                         bases, rprob, rbase);
            for (auto& edit : new_edits) {
                *new_mapping->add_edit() = edit;
            }
        }
    }
    // re-derive the alignment's sequence
    mutaln = simplify(mutaln);
    mutaln.set_sequence(alignment_seq(mutaln));
    return mutaln;
}

    void test2()
    {
        std::string string_a = "CAGCCCTAC";
        std::string string_b = "CCTGTACCC";

        std::vector<std::vector<int>> similarity = {
                {+2, +0, +0, +0, -1},
                {+0, +2, +0, +0, -1},
                {+0, +0, +2, +0, -1},
                {+0, +0, +0, +2, -1},
                {-1, -1, -1, -1, +2}
        };

        Alignment alignment = Alignment(string_a, string_b, Alignment::STR_WEIGHTS(1, -1, -1), similarity,
                                        [](int j) { return -(1); });

        std::pair<std::string, std::string> result = alignment.global_alignment();


        std::string expected_a = "CAGCCCTAC--";
        std::string expected_b = "C--CTGTACCC";

        if (result.first != expected_a || result.second != expected_b)
        {
            alignment.debugS();
            this->testFailed("Test2");
        }
    }

void irgen::applyLayoutAttributes(IRGenModule &IGM,
                                  CanType ASTTy,
                                  bool IsFixedLayout,
                                  Alignment &MinimumAlign) {
  assert(ASTTy && "shouldn't call applyLayoutAttributes without a type");
  
  auto &Diags = IGM.Context.Diags;
  auto decl = ASTTy->getAnyNominal();
  if (!decl)
    return;
  
  if (auto alignment = decl->getAttrs().getAttribute<AlignmentAttr>()) {
    auto value = alignment->getValue();
    assert(value != 0 && ((value - 1) & value) == 0
           && "alignment not a power of two!");
    
    if (!IsFixedLayout)
      Diags.diagnose(alignment->getLocation(),
                     diag::alignment_dynamic_type_layout_unsupported);
    else if (value < MinimumAlign.getValue())
      Diags.diagnose(alignment->getLocation(),
                   diag::alignment_less_than_natural, MinimumAlign.getValue());
    else
      MinimumAlign = Alignment(value);
  }
}

    void test5()
    {
        std::string string_a = "GTATT";
        std::string string_b = "TTT";

        std::vector<std::vector<int>> similarity = {
                //A   G   C   T   -
                {+1, -1, -1, -1, -1},
                {-1, +1, -1, -1, -1},
                {-1, -1, +1, -1, -1},
                {-1, -1, -1, +1, -1},
                {-1, -1, -1, -1, +1}
        };

        Alignment alignment = Alignment(string_a, string_b, Alignment::STR_WEIGHTS(1, -1, -1), similarity,
                                        [](int j) { return -(j+2); });

        int result = alignment.get_similarity();

        int expected_result = -3;
        if (result != expected_result)
        {
            alignment.debugS();
            this->testFailed("Test5");
        }
    }

data_STEMDIRECTION Layer::GetDrawingStemDir(const ArrayOfBeamElementCoords *coords)
{
    assert(!coords->empty());

    // Adjust the x position of the first and last element for taking into account the stem width
    LayerElement *first = dynamic_cast<LayerElement *>(coords->front()->m_element);
    LayerElement *last = dynamic_cast<LayerElement *>(coords->back()->m_element);

    if (!first || !last) {
        return m_drawingStemDir;
    }

    Measure *measure = dynamic_cast<Measure *>(this->GetFirstParent(MEASURE));
    assert(measure);

    // First check if there is any <space> in the measure - if not we can return the layer stem direction
    if (!measure->FindChildByType(SPACE)) {
        return m_drawingStemDir;
    }

    Alignment *alignmentFirst = first->GetAlignment();
    assert(alignmentFirst);
    Alignment *alignmentLast = last->GetAlignment();
    assert(alignmentLast);

    // We are ignoring cross-staff situation here because this should not be called if we have one
    Staff *staff = dynamic_cast<Staff *>(first->GetFirstParent(STAFF));
    assert(staff);

    double time = alignmentFirst->GetTime();
    double duration = alignmentLast->GetTime() - time + last->GetAlignmentDuration();
    duration = durRound(duration);

    return GetDrawingStemDir(time, duration, measure, staff->GetN());
}

Alignment Sampler::alignment_with_error(size_t length,
                                        double base_error,
                                        double indel_error) {
    size_t maxiter = 100;
    Alignment aln;
    if (base_error > 0 || indel_error > 0) {
        // sample a longer-than necessary alignment, then trim
        size_t iter = 0;
        while (iter++ < maxiter) {
            aln = mutate(
                alignment(length + 2 * ((double) length * indel_error)),
                base_error, indel_error);
            if (aln.sequence().size() == length) {
                break;
            } else if (aln.sequence().size() > length) {
                aln = strip_from_end(aln, aln.sequence().size() - length);
                break;
            }
        }
        if (iter == maxiter) {
            cerr << "[vg::Sampler] Warning: could not generate alignment of sufficient length. "
                 << "Graph may be too small, or indel rate too high." << endl;
        }
    } else {
        aln = alignment(length);
    }
    return aln;
}

Alignment SymForceAligner::cIntersection(int *a, int m, int* b, int n){
    Alignment out;
    for (int j = 1; j <= m; j++)
        if (a[j] && b[a[j]] == j)
            out.insert(AlignmentPoint(a[j]-1, j-1));
    return out;
}

    void test1()
    {
        std::string string_a = "GCATGCA";
        std::string string_b = "GATTACA";

        std::vector<std::vector<int>> similarity = {
                //A   G   C   T   -
                {+2, +0, +0, +0, -1},
                {+0, +2, +0, +0, -1},
                {+0, +0, +2, +0, -1},
                {+0, +0, +0, +2, -1},
                {-1, -1, -1, -1, +2},
        };

        Alignment alignment = Alignment(string_a, string_b, Alignment::STR_WEIGHTS(1, -1, -1), similarity,
                                        [](int j) { return -(j); });

        std::pair<std::string, std::string> result = alignment.global_alignment();

        std::string expected_a = "GCATG-CA";
        std::string expected_b = "G-ATTACA";
        std::cout << alignment.get_similarity();

        if (result.first != expected_a || result.second != expected_b)
        {
            alignment.debugS();
            this->testFailed("Test1");
        }
    }

void Pileups::compute_from_alignment(VG& graph, Alignment& alignment) {
    if (alignment.is_reverse()) {
        flip_alignment(alignment);
    }
    const Path& path = alignment.path();
    int64_t read_offset = 0;
    for (int i = 0; i < path.mapping_size(); ++i) {
        const Mapping& mapping = path.mapping(i);
        if (graph.has_node(mapping.position().node_id())) {
            const Node* node = graph.get_node(mapping.position().node_id());
            NodePileup* pileup = get_create(node->id());
            int64_t node_offset = mapping.position().offset();
            for (int j = 0; j < mapping.edit_size(); ++j) {
                const Edit& edit = mapping.edit(j);
                // process all pileups in edit.
                // update the offsets as we go
                compute_from_edit(*pileup, node_offset, read_offset, *node,
                                  alignment, mapping, edit);
            }
        }
    }
    assert(alignment.sequence().empty() ||
           alignment.path().mapping_size() == 0 ||
           read_offset == alignment.sequence().length());
}

int main()
{
  init();

  ifstream fi("res/streptococcus_references.fasta");
  vector<string> vs;
  string l;
  while (getline(fi, l)) {
    if (l[0] == '>') {
      vs.push_back(string());
    } else {
      vs.back() += l;
    }
  }
  fi.close();

  int cutoff = 10000;
  for (auto &s : vs) {
    s = s.substr(0, cutoff);
  }

  SquareAA algo;
  //NaiveCubeAA algo;
  Scoring sc;

  int start = clock();
  Alignment sol = algo.align(vs[0], vs[1], sc);
  cout << "Run time: " << (float)(clock() - start) / CLOCKS_PER_SEC << endl;

  sol.output(cout, vs[0], vs[1]);

  system("pause");
}

    void test0()
    {
        std::string string_a = "AGAGTCAATCCATAG";
        std::string string_b = "CAGAGGTCCATCATG";

        std::vector<std::vector<int>> similarity = {
                //A   G   C   T   -
                {+2, +0, +0, +0, -1},
                {+0, +2, +0, +0, -1},
                {+0, +0, +2, +0, -1},
                {+0, +0, +0, +2, -1},
                {-1, -1, -1, -1, +2},
        };

        Alignment alignment = Alignment(string_a, string_b, Alignment::STR_WEIGHTS(1, -1, -1), similarity);

        std::pair<std::string, std::string> result = alignment.global_alignment();

        std::string expected_a = "-AGAG-TCAATCCATAG";
        std::string expected_b = "CAGAGGTCCATC-AT-G";
        std::cout << alignment.get_similarity();

        if (result.first != expected_a || result.second != expected_b)
        {
            alignment.debugS();
            this->testFailed("Test0");
        }
    }

void QualAdjAligner::align(Alignment& alignment, Graph& g, bool print_score_matrices) {

    gssw_graph* graph = create_gssw_graph(g, 0, nullptr);

    const string& sequence = alignment.sequence();
    const string& quality = alignment.quality();

    if (quality.length() != sequence.length()) {
        cerr << "error:[Aligner] sequence and quality strings different lengths, cannot perform base quality adjusted alignmenterror:[Aligner] sequence and quality strings different lengths, cannot perform base quality adjusted alignment" << endl;
    }

    gssw_graph_fill_qual_adj(graph, sequence.c_str(), quality.c_str(),
                             nt_table, adjusted_score_matrix,
                             scaled_gap_open, scaled_gap_extension, 15, 2);

    gssw_graph_mapping* gm = gssw_graph_trace_back_qual_adj (graph,
                             sequence.c_str(),
                             quality.c_str(),
                             sequence.size(),
                             nt_table,
                             adjusted_score_matrix,
                             scaled_gap_open,
                             scaled_gap_extension);

    gssw_mapping_to_alignment(graph, gm, alignment, print_score_matrices);

#ifdef debug
    gssw_print_graph_mapping(gm, stderr);
#endif

    gssw_graph_mapping_destroy(gm);
    gssw_graph_destroy(graph);
}

    void test4()
    {
        std::string string_a = "AGAGTCAATCCATAG";
        std::string string_b = "CAGAGGTCCATCATG";
        std::vector<std::vector<int>> similarity = {
                {+2, +0, +0, +0, -1},
                {+0, +2, +0, +0, -1},
                {+0, +0, +2, +0, -1},
                {+0, +0, +0, +2, -1},
                {-1, -1, -1, -1, +2}
        };

        Alignment alignment = Alignment(string_a, string_b, Alignment::STR_WEIGHTS(1, -1, -1), similarity,
                                        [](int j) { return -(j+2); });

        std::pair<std::string, std::string> result = alignment.global_alignment();

        std::string expected_a = "-AGAG-TCAATCCATAG";
        std::string expected_b = "CAGAGGTCCATC-AT-G";

        if (result.first != expected_a || result.second != expected_b)
        {
            alignment.debugS();
            this->testFailed("Test4");
        }
    }

    /**
     * Filter reads that are less than <PCTID> reference.
     * I.E. if a read matches the reference along 80% of its
     * length, and your cutoff is 90% PCTID, throw it out.
     */
    Alignment Filter::percent_identity_filter(Alignment& aln){
        double read_pctid = 0.0;
        //read pct_id = len(matching sequence / len(total sequence)

        int64_t aln_total_len = aln.sequence().size();
        int64_t aln_match_len = 0;

        std::function<double(int64_t, int64_t)> calc_pct_id = [](int64_t rp, int64_t ttlp){
            return ((double) rp / (double) ttlp);
        };



        Path path = aln.path();
        //TODO handle reversing mappings

        for (int i = 0; i < path.mapping_size(); i++){
            Mapping mapping = path.mapping(i);

            for (int j = 0; j < mapping.edit_size(); j++){
                Edit ee = mapping.edit(j);
                if (ee.from_length() == ee.to_length() && ee.sequence() == ""){
                    aln_match_len += ee.to_length();
                }

            }
        }
        if (calc_pct_id(aln_match_len, aln_total_len) < min_percent_identity){
            return inverse ? aln : Alignment();
        }

        return inverse ? Alignment() : aln;


    }