Performs local alignment

@param alignmentScoring scoring system
@param seq1             first sequence
@param seq2             second sequence
@return result of alignment with information about alignment positions in both sequences and array of mutations
519147
public static <S extends Sequence<S>> Alignment<S> alignLocal(AlignmentScoring<S> alignmentScoring,
                                                                  S seq1, S seq2) {
        if (alignmentScoring instanceof AffineGapAlignmentScoring)
            return alignLocalAffine((AffineGapAlignmentScoring<S>) alignmentScoring, seq1, seq2);
        if (alignmentScoring instanceof LinearGapAlignmentScoring)
            return alignLocalLinear((LinearGapAlignmentScoring<S>) alignmentScoring, seq1, seq2);
        throw new RuntimeException("Unknown scoring type.");
    }
852623
def build_alignment(self,score,pieces):
		""""""
	 	# build text
		self.open_seqs()
		text1 = text2 = ""
		end1  = end2  = None
		for (start1,start2,length,pctId) in pieces:
			if (end1 != None):
				if (start1 == end1): # insertion in sequence 2
					text1 += self.seq1_gap * (start2-end2)
					text2 += self.seq2_file.get(end2,start2-end2)
				else: # insertion in sequence 1
					text1 += self.seq1_file.get(end1,start1-end1)
					text2 += self.seq2_gap * (start1-end1)

			text1 += self.seq1_file.get(start1,length)
			text2 += self.seq2_file.get(start2,length)
			end1 = start1 + length
			end2 = start2 + length
		# create alignment
		start1 = pieces[0][0]
		start2 = pieces[0][1]
		end1   = pieces[-1][0] + pieces[-1][2]
		end2   = pieces[-1][1] + pieces[-1][2]
		size1  = end1 - start1
		size2  = end2 - start2
		a = Alignment(score=score,species_to_lengths=self.species_to_lengths)
		#if (self.seq1_strand == "-"): start1 = self.seq1_file.length - end1
		a.add_component(Component(self.seq1_src,start1,size1,self.seq1_strand,text=text1))
		#if (self.seq2_strand == "-"): start2 = self.seq2_file.length - end2
		a.add_component(Component(self.seq2_src,start2,size2,self.seq2_strand,text=text2))
		return a
1058525
public static <S extends Sequence<S>> Alignment<S> alignGlobal(AlignmentScoring<S> alignmentScoring,
                                                                   S seq1, S seq2) {
        if (alignmentScoring instanceof AffineGapAlignmentScoring)
            return alignGlobalAffine((AffineGapAlignmentScoring<S>) alignmentScoring, seq1, seq2);
        if (alignmentScoring instanceof LinearGapAlignmentScoring)
            return alignGlobalLinear((LinearGapAlignmentScoring<S>) alignmentScoring, seq1, seq2);
        throw new RuntimeException("Unknown scoring type.");
    }
369231
private final float[] align(final int sResidue, final int dIndex) {
	int dResidue, r;
	float maxScore = -1000000;
	float rho1 = 0;
	int maxIdx = 0;
	float rho0 = 0;
	short[] dbAARow = model.dbAA[dIndex];
	int numOfIterations = model.Length[dIndex] - ORonnModel.AA_ALPHABET;
	for (dResidue = 0; dResidue <= numOfIterations; dResidue++) {
		// go though the database sequence for maximised alignment
		rho1 = 0.0f;
		for (r = 0; r < ORonnModel.AA_ALPHABET; r++) {
		// go through the query sequence for one alignment
		rho1 += RonnConstraint.Blosum62[seqAA[sResidue + r]][dbAARow[dResidue
			+ r]];
		}
		if (rho1 > maxScore) {
		maxScore = rho1;
		maxIdx = dResidue;
		}
	}
	for (r = 0; r < ORonnModel.AA_ALPHABET; r++) {
		rho0 += RonnConstraint.Blosum62[dbAARow[maxIdx + r]][dbAARow[maxIdx
			+ r]];
	}
	return new float[] { rho0, maxScore };
	}
141796
def align_two_alignments(aln1, aln2, moltype, params=None):
    
    #create SequenceCollection object from seqs
    aln1 = Alignment(aln1,MolType=moltype)
    #Create mapping between abbreviated IDs and full IDs
    aln1_int_map, aln1_int_keys = aln1.getIntMap()
    #Create SequenceCollection from int_map.
    aln1_int_map = Alignment(aln1_int_map,MolType=moltype)

    #create Alignment object from aln
    aln2 = Alignment(aln2,MolType=moltype)
    #Create mapping between abbreviated IDs and full IDs
    aln2_int_map, aln2_int_keys = aln2.getIntMap(prefix='seqn_')
    #Create SequenceCollection from int_map.
    aln2_int_map = Alignment(aln2_int_map,MolType=moltype)

    #Update aln1_int_keys with aln2_int_keys
    aln1_int_keys.update(aln2_int_keys)

    #Create Mafft app.
    app = Clustalw(InputHandler='_input_as_multiline_string',\
        params=params,
        SuppressStderr=True)
    app.Parameters['-align'].off()
    app.Parameters['-infile'].off()
    app.Parameters['-profile'].on()

    #Add aln_int_map as profile1
    app.Parameters['-profile1'].on(\
        app._tempfile_as_multiline_string(aln1_int_map.toFasta()))

    #Add seq_int_map as profile2
    app.Parameters['-profile2'].on(\
        app._tempfile_as_multiline_string(aln2_int_map.toFasta()))
    #Get results using int_map as input to app
    res = app()

    #Get alignment as dict out of results
    alignment = dict(ClustalParser(res['Align'].readlines()))

    #Make new dict mapping original IDs
    new_alignment = {}
    for k,v in alignment.items():
        new_alignment[aln1_int_keys[k]]=v
    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment,MolType=moltype)
    #Clean up
    res.cleanUp()
    remove(app.Parameters['-profile1'].Value)
    remove(app.Parameters['-profile2'].Value)
    del(aln1,aln1_int_map,aln1_int_keys,\
        aln2,aln2_int_map,aln2_int_keys,app,res,alignment)

    return new_alignment
1035696
public static <S extends Sequence<S>> Alignment<S> alignGlobalAffine(AffineGapAlignmentScoring<S> scoring,
                                                                         S seq1, S seq2) {
        if (seq1.getAlphabet() != seq2.getAlphabet() || seq1.getAlphabet() != scoring.getAlphabet())
            throw new IllegalArgumentException("Different alphabets.");

        int
                length1 = seq1.size(),
                length2 = seq2.size(),
                size1 = length1 + 1,
                size2 = length2 + 1;
        Matrix
                gapIn1 = new Matrix(size1, size2),
                gapIn2 = new Matrix(size1, size2),
                matrix = new Matrix(size1, size2);

        int
                gapExtensionPenalty = scoring.getGapExtensionPenalty(),
                gapOpenPenalty = scoring.getGapOpenPenalty() - gapExtensionPenalty;

        int i, j;
        for (i = 1; i < size1; ++i) {
            matrix.set(i, 0, gapOpenPenalty + i * gapExtensionPenalty);
            gapIn2.set(i, 0, gapOpenPenalty + i * gapExtensionPenalty);
            gapIn1.set(i, 0, MIN_VALUE);

        }
        for (j = 1; j < size2; ++j) {
            matrix.set(0, j, gapOpenPenalty + j * gapExtensionPenalty);
            gapIn1.set(0, j, gapOpenPenalty + j * gapExtensionPenalty);
            gapIn2.set(0, j, MIN_VALUE);
        }

        matrix.set(0, 0, 0);
        gapIn1.set(0, 0, MIN_VALUE);
        gapIn2.set(0, 0, MIN_VALUE);

        for (i = 1; i <= length1; ++i) {
            for (j = 1; j <= length2; ++j) {
                gapIn1.set(i, j, max(
                        matrix.get(i, j - 1) + gapOpenPenalty + gapExtensionPenalty,
                        gapIn1.get(i, j - 1) + gapExtensionPenalty));

                gapIn2.set(i, j, max(
                        matrix.get(i - 1, j) + gapOpenPenalty + gapExtensionPenalty,
                        gapIn2.get(i - 1, j) + gapExtensionPenalty));

                matrix.set(i, j, max(
                        matrix.get(i - 1, j - 1) +
                                scoring.getScore(seq1.codeAt(i - 1), seq2.codeAt(j - 1)),
                        gapIn2.get(i, j),
                        gapIn1.get(i, j)));
            }
        }

        i = length1;
        j = length2;

        MutationsBuilder<S> mutations = new MutationsBuilder<>(seq1.getAlphabet(), true);
        int score = matrix.get(length1, length2);
        boolean inGap1 = false, inGap2 = false;
        while (i > 0 || j > 0) {
            assert !inGap1 || !inGap2;
            if (!inGap2 && (inGap1 || (j > 0 && score == gapIn1.get(i, j)))) {
                inGap1 = false;
                if (score == gapIn1.get(i, j - 1) + gapExtensionPenalty) {
                    inGap1 = true;
                    score = gapIn1.get(i, j - 1);
                } else
                    score = matrix.get(i, j - 1);


                mutations.appendInsertion(i, seq2.codeAt(j - 1));
                j--;
            } else if (inGap2 || (i > 0 && score == gapIn2.get(i, j))) {
                inGap2 = false;
                if (score == gapIn2.get(i - 1, j) + gapExtensionPenalty) {
                    inGap2 = true;
                    score = gapIn2.get(i - 1, j);
                } else
                    score = matrix.get(i - 1, j);


                mutations.appendDeletion(i - 1, seq1.codeAt(i - 1));
                i--;
            } else if (i > 0 && j > 0
                    && score == matrix.get(i - 1, j - 1)
                    + scoring.getScore(seq1.codeAt(i - 1), seq2.codeAt(j - 1))) {
                score = matrix.get(i - 1, j - 1);
                if (seq1.codeAt(i - 1) != seq2.codeAt(j - 1))
                    mutations.appendSubstitution(i - 1, seq1.codeAt(i - 1), seq2.codeAt(j - 1));

                --i;
                --j;
            } else
                throw new RuntimeException();
        }

        return new Alignment<>(seq1, mutations.createAndDestroy(),
                new Range(0, seq1.size()), new Range(0, seq2.size()), matrix.get(length1, length2));
    }
931328
public static boolean isNonAmbNucleotideSequence(String sequence) {
	sequence = SequenceUtil.cleanSequence(sequence);
	if (SequenceUtil.DIGIT.matcher(sequence).find()) {
		return false;
	}
	if (SequenceUtil.NON_NUCLEOTIDE.matcher(sequence).find()) {
		return false;
		/*
		 * System.out.format("I found the text starting at " +
		 * "index %d and ending at index %d.%n", nonDNAmatcher .start(),
		 * nonDNAmatcher.end());
		 */
	}
	final Matcher DNAmatcher = SequenceUtil.NUCLEOTIDE.matcher(sequence);
	return DNAmatcher.find();
	}
486473
def tiered_alignment(in_bam, tier_num, multi_mappers, extra_args,
                     genome_build, pair_stats,
                     work_dir, dirs, config):
    
    nomap_fq1, nomap_fq2 = select_unaligned_read_pairs(in_bam, "tier{}".format(tier_num),
                                                       work_dir, config)
    if nomap_fq1 is not None:
        base_name = "{}-tier{}out".format(os.path.splitext(os.path.basename(in_bam))[0],
                                          tier_num)
        config = copy.deepcopy(config)
        dirs = copy.deepcopy(dirs)
        config["algorithm"]["bam_sort"] = "queryname"
        config["algorithm"]["multiple_mappers"] = multi_mappers
        config["algorithm"]["extra_align_args"] = ["-i", int(pair_stats["mean"]),
                                               int(pair_stats["std"])] + extra_args
        out_bam, ref_file = align_to_sort_bam(nomap_fq1, nomap_fq2,
                                              lane.rg_names(base_name, base_name, config),
                                              genome_build, "novoalign",
                                              dirs, config,
                                              dir_ext=os.path.join("hydra", os.path.split(nomap_fq1)[0]))
        return out_bam
    else:
        return None
406528
public AlignmentHelper getAlignmentHelper() {
        List<Boolean> matches = new ArrayList<>();

        IntArrayList pos1 = new IntArrayList(sequence1.size() + mutations.size()),
                pos2 = new IntArrayList(sequence1.size() + mutations.size());

        StringBuilder sb1 = new StringBuilder(),
                sb2 = new StringBuilder();

        Alphabet<S> alphabet = mutations.getAlphabet();

        AlignmentIteratorForward<S> iterator = forwardIterator();
        while (iterator.advance()) {
            final int mut = iterator.getCurrentMutation();
            switch (getRawTypeCode(mut)) {
                case RAW_MUTATION_TYPE_SUBSTITUTION:
                    pos1.add(iterator.getSeq1Position());
                    pos2.add(iterator.getSeq2Position());
                    sb1.append(sequence1.symbolAt(iterator.getSeq1Position()));
                    sb2.append(Mutation.getToSymbol(mut, alphabet));
                    matches.add(false);
                    break;

                case RAW_MUTATION_TYPE_DELETION:
                    pos1.add(iterator.getSeq1Position());
                    pos2.add(-1 - iterator.getSeq2Position());
                    sb1.append(sequence1.symbolAt(iterator.getSeq1Position()));
                    sb2.append("-");
                    matches.add(false);
                    break;

                case RAW_MUTATION_TYPE_INSERTION:
                    pos1.add(-1 - iterator.getSeq1Position());
                    pos2.add(iterator.getSeq2Position());
                    sb1.append("-");
                    sb2.append(Mutation.getToSymbol(mut, alphabet));
                    matches.add(false);
                    break;

                default:
                    pos1.add(iterator.getSeq1Position());
                    pos2.add(iterator.getSeq2Position());
                    char c = sequence1.symbolAt(iterator.getSeq1Position());
                    sb1.append(c);
                    sb2.append(c);
                    matches.add(true);
                    break;
            }
        }

        return new AlignmentHelper(sb1.toString(), sb2.toString(),
                pos1.toArray(), pos2.toArray(),
                new BitArray(matches));
    }
1000013
def seqannotation(self, seqrecord, allele, loc):
        
        #seqrecord = self.seqrecord(allele, loc)
        complete_annotation = get_features(seqrecord)
        annotation = Annotation(annotation=complete_annotation,
                                method='match',
                                complete_annotation=True)

        if self.alignments:
            alignment = {f: self.annoated_alignments[loc][allele][f]['Seq']
                         for f in self.annoated_alignments[loc][allele].keys()}
            annotation.aligned = alignment

        return annotation
63176
public KAlignerParameters2 setScoring(AffineGapAlignmentScoring<NucleotideSequence> scoring) {
        if (scoring != null && !scoring.uniformBasicMatchScore())
            throw new IllegalArgumentException("Use scoring with common match score.");
        this.scoring = scoring;
        return this;
    }
301176
def compare_genomes(fastas, id2desc, algorithm):
    
    hits = {}
    for pair in itertools.permutations(fastas, 2):
        if pair[0] not in hits:
            hits[pair[0]] = {}
        hits[pair[0]][pair[1]] = {}
    for pair in itertools.combinations(fastas, 2):
        A, B = pair
        blastF = open(search(A, B, method = algorithm, alignment = 'local'))
        blastR = open(search(B, A, method = algorithm, alignment = 'local'))
        for compare in [A, B, blastF], [B, A, blastR]:
            query, ref, blast = compare
            for hit in best_blast(blast, 1):
                a, b, pident, qstart, qend, e, bit = \
                        hit[0].split()[0], hit[1].split()[0], float(hit[2]), float(hit[6]), float(hit[7]), float(hit[-2]), float(hit[-1])
                alignment_length = abs(qstart - qend + 1)
                norm_bit = normalize_bit(a, b, bit, id2desc)
                length_fraction = float((alignment_length) / id2desc[a][2])
                hits[query][ref][a] = ['fbh', a, b, pident, length_fraction, e, bit, norm_bit]
    return hits
32178
public static function fromASN1(Sequence $seq): self
    {
        $idx = 0;
        if ($seq->hasTagged(0)) {
            $idx++;
            $version = $seq->getTagged(0)
                ->asExplicit()
                ->asInteger()
                ->intNumber();
        } else {
            $version = self::VERSION_1;
        }
        $serial = $seq->at($idx++)
            ->asInteger()
            ->number();
        $algo = AlgorithmIdentifier::fromASN1($seq->at($idx++)->asSequence());
        if (!$algo instanceof SignatureAlgorithmIdentifier) {
            throw new \UnexpectedValueException(
                "Unsupported signature algorithm " . $algo->name() . ".");
        }
        $issuer = Name::fromASN1($seq->at($idx++)->asSequence());
        $validity = Validity::fromASN1($seq->at($idx++)->asSequence());
        $subject = Name::fromASN1($seq->at($idx++)->asSequence());
        $pki = PublicKeyInfo::fromASN1($seq->at($idx++)->asSequence());
        $tbs_cert = new self($subject, $pki, $issuer, $validity);
        $tbs_cert->_version = $version;
        $tbs_cert->_serialNumber = $serial;
        $tbs_cert->_signature = $algo;
        if ($seq->hasTagged(1)) {
            $tbs_cert->_issuerUniqueID = UniqueIdentifier::fromASN1(
                $seq->getTagged(1)
                    ->asImplicit(Element::TYPE_BIT_STRING)
                    ->asBitString());
        }
        if ($seq->hasTagged(2)) {
            $tbs_cert->_subjectUniqueID = UniqueIdentifier::fromASN1(
                $seq->getTagged(2)
                    ->asImplicit(Element::TYPE_BIT_STRING)
                    ->asBitString());
        }
        if ($seq->hasTagged(3)) {
            $tbs_cert->_extensions = Extensions::fromASN1(
                $seq->getTagged(3)
                    ->asExplicit()
                    ->asSequence());
        }
        return $tbs_cert;
    }
797876
public static function fromV2ASN1(Sequence $seq): self
    {
        $issuer = null;
        $cert_id = null;
        $digest_info = null;
        if ($seq->has(0, Element::TYPE_SEQUENCE)) {
            $issuer = GeneralNames::fromASN1($seq->at(0)->asSequence());
        }
        if ($seq->hasTagged(0)) {
            $cert_id = IssuerSerial::fromASN1(
                $seq->getTagged(0)
                    ->asImplicit(Element::TYPE_SEQUENCE)
                    ->asSequence());
        }
        if ($seq->hasTagged(1)) {
            $digest_info = ObjectDigestInfo::fromASN1(
                $seq->getTagged(1)
                    ->asImplicit(Element::TYPE_SEQUENCE)
                    ->asSequence());
        }
        $obj = new self($issuer);
        $obj->_baseCertificateID = $cert_id;
        $obj->_objectDigestInfo = $digest_info;
        return $obj;
    }
85435
public static <S extends Sequence<S>> int calculateScore(S seq1, Range seq1Range, Mutations<S> mutations,
                                                             LinearGapAlignmentScoring<S> scoring) {
        if (!mutations.isEmpty() && mutations.getPositionByIndex(0) < seq1Range.getFrom() - 1)
            throw new IllegalArgumentException();

        final AlignmentIteratorForward<S> iterator = new AlignmentIteratorForward<>(mutations, seq1Range);

        int score = 0;

        while (iterator.advance()) {
            final int mut = iterator.getCurrentMutation();
            switch (Mutation.getRawTypeCode(mut)) {
                case RAW_MUTATION_TYPE_SUBSTITUTION:
                    score += scoring.getScore(Mutation.getFrom(mut), Mutation.getTo(mut));
                    break;

                case RAW_MUTATION_TYPE_DELETION:
                case RAW_MUTATION_TYPE_INSERTION:
                    score += scoring.getGapPenalty();
                    break;

                default:
                    byte c = seq1.codeAt(iterator.getSeq1Position());
                    score += scoring.getScore(c, c);
                    break;
            }
        }

        return score;
    }
265057
public static <S extends Sequence<S>> Alignment<S> shiftIndelsAtHomopolymers(Alignment<S> alignment) {
        return new Alignment<>(alignment.sequence1,
                MutationsUtil.shiftIndelsAtHomopolymers(alignment.sequence1, alignment.sequence1Range.getFrom(), alignment.mutations),
                alignment.sequence1Range, alignment.sequence2Range, alignment.score);
    }
958031
def alignments(self):
        
        referenceIds = self.referenceIds
        dropUnmapped = self.dropUnmapped
        dropSecondary = self.dropSecondary
        dropSupplementary = self.dropSupplementary
        dropDuplicates = self.dropDuplicates
        keepQCFailures = self.keepQCFailures
        storeQueryIds = self.storeQueryIds
        filterRead = self.filterRead
        minScore = self.minScore
        maxScore = self.maxScore
        scoreTag = self.scoreTag

        if storeQueryIds:
            self.queryIds = queryIds = set()

        lastAlignment = None
        count = 0
        with samfile(self.filename) as samAlignment:
            for count, alignment in enumerate(samAlignment.fetch(), start=1):
                if storeQueryIds:
                    queryIds.add(alignment.query_name)

                if minScore is not None or maxScore is not None:
                    try:
                        score = alignment.get_tag(scoreTag)
                    except KeyError:
                        continue
                    else:
                        if ((minScore is not None and score < minScore) or
                                (maxScore is not None and score > maxScore)):
                            continue

                # Secondary and supplementary alignments may have a '*'
                # (pysam returns this as None) SEQ field, indicating that
                # the previous sequence should be used. This is best
                # practice according to section 2.5.2 of
                # https://samtools.github.io/hts-specs/SAMv1.pdf So we use
                # the last alignment query and quality strings if we get
                # None as a query sequence.
                if alignment.query_sequence is None:
                    if lastAlignment is None:
                        raise InvalidSAM(
                            'pysam produced an alignment (number %d) with no '
                            'query sequence without previously giving an '
                            'alignment with a sequence.' % count)
                    # Use the previous query sequence and quality. I'm not
                    # making the call to _hardClip dependent on
                    # alignment.cigartuples (as in the else clause below)
                    # because I don't think it's possible for
                    # alignment.cigartuples to be None in this case. If we
                    # have a second match on a query, then it must be
                    # aligned to something (i.e., it cannot be unmapped
                    # with no CIGAR string). The assertion will tell us if
                    # this is ever not the case.
                    assert alignment.cigartuples
                    (alignment.query_sequence,
                     alignment.query_qualities, _) = _hardClip(
                         lastAlignment.query_sequence,
                         lastAlignment.query_qualities,
                         alignment.cigartuples)
                else:
                    lastAlignment = alignment
                    if alignment.cigartuples:
                        (alignment.query_sequence,
                         alignment.query_qualities, _) = _hardClip(
                             alignment.query_sequence,
                             alignment.query_qualities,
                             alignment.cigartuples)

                if ((filterRead is None or
                     filterRead(Read(alignment.query_name,
                                     alignment.query_sequence,
                                     alignment.qual))) and
                    not (
                        (referenceIds and
                         alignment.reference_name not in referenceIds) or
                        (alignment.is_unmapped and dropUnmapped) or
                        (alignment.is_secondary and dropSecondary) or
                        (alignment.is_supplementary and dropSupplementary) or
                        (alignment.is_duplicate and dropDuplicates) or
                        (alignment.is_qcfail and not keepQCFailures))):
                    yield alignment

        self.alignmentCount = count