Return the alignment stored in self.sequence_alignments given a seqprop + another seqprop
979960
def _get_seqprop_to_seqprop_alignment(self, seqprop1, seqprop2):
        """"""
        if isinstance(seqprop1, str):
            seqprop1_id = seqprop1
        else:
            seqprop1_id = seqprop1.id
        if isinstance(seqprop2, str):
            seqprop2_id = seqprop2
        else:
            seqprop2_id = seqprop2.id

        aln_id = '{}_{}'.format(seqprop1_id, seqprop2_id)

        if self.sequence_alignments.has_id(aln_id):
            alignment = self.sequence_alignments.get_by_id(aln_id)
            return alignment
        else:
            raise ValueError('{}: sequence alignment not found, please run the alignment first'.format(aln_id))
959650
public static <S extends Sequence<S>> Alignment<S> alignGlobalLinear(LinearGapAlignmentScoring scoring,
                                                                         S seq1, S seq2) {
        if (seq1.getAlphabet() != seq2.getAlphabet() ||
                seq1.getAlphabet() != scoring.getAlphabet())
            throw new IllegalArgumentException("Different alphabets.");

        int size1 = seq1.size() + 1,
                size2 = seq2.size() + 1;
        int[] matrix = new int[size1 * (seq2.size() + 1)];

        for (int i = 0; i < size2; ++i)
            matrix[i] = scoring.getGapPenalty() * i;

        for (int j = 1; j < size1; ++j)
            matrix[size2 * j] = scoring.getGapPenalty() * j;

        int i1, i2,
                match, delete, insert;

        for (i1 = 0; i1 < seq1.size(); ++i1)
            for (i2 = 0; i2 < seq2.size(); ++i2) {
                match = matrix[i1 * size2 + i2] +
                        scoring.getScore(seq1.codeAt(i1), seq2.codeAt(i2));
                delete = matrix[i1 * size2 + i2 + 1] + scoring.getGapPenalty();
                insert = matrix[(i1 + 1) * size2 + i2] + scoring.getGapPenalty();
                matrix[(i1 + 1) * size2 + i2 + 1] = Math.max(match, Math.max(delete, insert));
            }

        MutationsBuilder<S> builder = new MutationsBuilder<>(seq1.getAlphabet(), true);

        i1 = seq1.size() - 1;
        i2 = seq2.size() - 1;
        int score = matrix[(i1 + 1) * size2 + i2 + 1];

        while (i1 >= 0 || i2 >= 0) {
            if (i1 >= 0 && i2 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] == matrix[i1 * size2 + i2] +
                            scoring.getScore(seq1.codeAt(i1), seq2.codeAt(i2))) {
                if (seq1.codeAt(i1) != seq2.codeAt(i2))
                    builder.appendSubstitution(i1, seq1.codeAt(i1), seq2.codeAt(i2));
                --i1;
                --i2;
            } else if (i1 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] ==
                            matrix[i1 * size2 + i2 + 1] + scoring.getGapPenalty()) {
                builder.appendDeletion(i1, seq1.codeAt(i1));
                i1--;
            } else if (i2 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] ==
                            matrix[(i1 + 1) * size2 + i2] + scoring.getGapPenalty()) {
                builder.appendInsertion(i1 + 1, seq2.codeAt(i2));
                i2--;
            } else
                throw new RuntimeException();
        }

        return new Alignment<>(seq1, builder.createAndDestroy(), new Range(0, seq1.size()), new Range(0, seq2.size()),
                score);
    }
371979
def inferred_var_seqs_plus_flanks(self, ref_seq, flank_length):
        ''''''
        flank_start = max(0, self.POS - flank_length)
        flank_end = min(len(ref_seq) - 1, self.ref_end_pos() + flank_length)
        seqs = [ref_seq[flank_start:self.POS] + self.REF + ref_seq[self.ref_end_pos() + 1: flank_end + 1]]

        for alt in self.ALT:
            seqs.append(ref_seq[flank_start:self.POS] + alt + ref_seq[self.ref_end_pos() + 1: flank_end + 1])

        return flank_start, seqs
678433
def aligned_read_to_read(read, revcomp=True, qual=None, ignore_quality=False):
    ''''''
    if read.qual is None or ignore_quality:
        if qual is None or ignore_quality:
            seq = pyfastaq.sequences.Fasta(read.qname, common.decode(read.seq))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)
    else:
        if qual is None:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), common.decode(read.qual))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)

    if read.is_reverse and revcomp:
        seq.revcomp()

    return seq
297588
def _get_alignment_lines(self):
        ''' '''

        # Strip the boilerplate lines
        lines = self.alignment_output.split("\n")
        assert(lines[0].startswith('CLUSTAL'))
        lines = '\n'.join(lines[1:]).lstrip().split('\n')

        # The sequence IDs should be unique. Reassert this here
        assert(len(self.sequence_ids.values()) == len(set(self.sequence_ids.values())))

        # Create the list of sequence IDs
        id_list = [v for k, v in sorted(self.sequence_ids.iteritems())]

        # Determine the indentation level
        first_id = id_list[0]
        header_regex = re.compile("(.*?\s+)(.*)")
        alignment_regex = re.compile("^([A-Z\-]+)\s*$")
        mtchs = header_regex.match(lines[0])
        assert(mtchs.group(1).strip() == first_id)
        indentation = len(mtchs.group(1))
        sequence = mtchs.group(2)
        assert(sequence)
        assert(alignment_regex.match(sequence))

        # Create empty lists for the sequences
        sequences = {}
        for id in id_list:
            sequences[id] = []
        sequences[-1] = []

        # Get the lists of sequences
        num_ids = len(id_list)
        for x in range(0, len(lines), num_ids + 2):
            for y in range(num_ids):
                id = id_list[y]
                assert(lines[x + y][:indentation].strip() == id)
                assert(lines[x + y][indentation - 1] == ' ')
                sequence = lines[x + y][indentation:].strip()
                assert(alignment_regex.match(sequence))
                sequences[id].append(sequence)

            # Get the length of the sequence lines
            length_of_sequences = list(set(map(len, [v[-1] for k, v in sequences.iteritems() if k != -1])))
            assert(len(length_of_sequences) == 1)
            length_of_sequences = length_of_sequences[0]

            # Parse the Clustal match line
            assert(lines[x + num_ids][:indentation].strip() == '')
            match_sequence = lines[x + num_ids][indentation:indentation + length_of_sequences]
            assert(match_sequence.strip() == lines[x + num_ids].strip())
            assert(lines[x + y][indentation - 1] == ' ')
            sequences[-1].append(match_sequence)

            # Check for the empty line
            assert(lines[x + num_ids + 1].strip() == '')

        # Create the sequences, making sure that all sequences are the same length
        lengths = set()
        for k, v in sequences.iteritems():
            sequences[k] = "".join(v)
            lengths.add(len(sequences[k]))
        assert(len(lengths) == 1)

        return sequences
835536
def alignment_to_partials(alignment, missing_data=None):
    
    partials_dict = {}
    for (name, sequence) in alignment.get_sequences():
        datatype = 'dna' if alignment.is_dna() else 'protein'
        partials_dict[name] = seq_to_partials(sequence, datatype)

    if missing_data is not None:
        l = len(alignment)
        for name in missing_data:
            if name not in partials_dict:
                partials_dict[name] = seq_to_partials('-'*l, datatype)
    return partials_dict
1149759
def _iterative_bridged_contig_pair_merge(self, outprefix):
        ''''''
        if self.reads is None:
            if self.verbose:
                print('Skipping iterative contig merging because no reads given (see --reads option)')
            return self.original_fasta, None, None

        log_file = outprefix + '.iterations.log'
        log_fh = pyfastaq.utils.open_file_write(log_file)
        genome_fasta = self.original_fasta
        nucmer_coords = outprefix + '.iter.1.coords'
        reads_to_map = self.reads
        act_prefix = None
        made_a_join = True
        iteration = 1

        while made_a_join:
            this_log_prefix = '[' + self.log_prefix + ' iterative_merge ' + str(iteration) + ']'
            print(this_log_prefix, '\tUsing nucmer matches from ', nucmer_coords, sep='', file=log_fh)
            self._run_nucmer(genome_fasta, self.reassembly.contigs_fasta, nucmer_coords)
            act_prefix = outprefix + '.iter.' + str(iteration)
            print(this_log_prefix, '\tYou can view the nucmer matches with ACT using: ./', act_prefix, '.start_act.sh', sep='', file=log_fh)
            self._write_act_files(genome_fasta, self.reassembly.contigs_fasta, nucmer_coords, act_prefix)
            nucmer_hits_by_ref = self._load_nucmer_hits(nucmer_coords)
            made_a_join = self._merge_all_bridged_contigs(nucmer_hits_by_ref, self.original_contigs, self.reassembly_contigs, log_fh, this_log_prefix)
            iteration += 1

            if made_a_join:
                print(this_log_prefix, '\tMade at least one merge. Remapping reads and reassembling',sep='', file=log_fh)
                nucmer_coords = outprefix + '.iter.' + str(iteration) + '.coords'
                genome_fasta = outprefix + '.iter.' + str(iteration) + '.merged.fasta'
                self._contigs_dict_to_file(self.original_contigs, genome_fasta)
                bam = outprefix + '.iter.' + str(iteration) + '.bam'

                circlator.mapping.bwa_mem(
                  genome_fasta,
                  reads_to_map,
                  bam,
                  threads=self.threads,
                  verbose=self.verbose,
                )

                reads_prefix = outprefix + '.iter.' + str(iteration) + '.reads'
                reads_to_map =  reads_prefix + ('.fasta' if self.spades_only_assembler else '.fastq')
                bam_filter = circlator.bamfilter.BamFilter(bam, reads_prefix, fastq_out=not self.spades_only_assembler, split_all_reads=self.split_all_reads)
                bam_filter.run()
                assembler_dir = outprefix + '.iter.' + str(iteration) + '.assembly'
                a = circlator.assemble.Assembler(
                    reads_prefix + ('.fasta' if self.spades_only_assembler else '.fastq'),
                    assembler_dir,
                    threads=self.threads,
                    careful=self.spades_careful,
                    only_assembler=self.spades_only_assembler,
                    verbose=self.verbose,
                    spades_kmers=self.spades_kmers,
                    spades_use_first_success=self.spades_use_first_success,
                    assembler=self.assembler,
                    genomeSize=self.length_cutoff,
                    data_type=self.data_type
                )
                a.run()
                self.reassembly = circlator.assembly.Assembly(assembler_dir, assembler=self.assembler)
                self.reassembly_contigs = self.reassembly.get_contigs()
            elif iteration <= 2:
                print(this_log_prefix, '\tNo contig merges were made',sep='', file=log_fh)

        pyfastaq.utils.close(log_fh)
        return genome_fasta, nucmer_coords, act_prefix + '.start_act.sh'
1344089
def align_orthologous_genes_pairwise(self, sc=None, joblib=False, cores=1, gapopen=10, gapextend=0.5,
                                         engine='needle', parse=True, force_rerun=False):
        """"""
        log.info('Aligning sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _align_orthologous_gene_pairwise_sc(g_id, g_to_pickle=self.gene_protein_pickles,
                                                gapopen=gapopen, gapextend=gapextend, engine=engine, parse=parse,
                                                outdir=self.sequences_by_gene_dir,
                                                force_rerun=force_rerun):
            """Align orthologous strain sequences to representative Protein sequence, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            import os.path as op

            protein_seqs_aln_pickle_path = op.join(outdir, '{}_protein_withseqs_dis_aln.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_aln_pickle_path):
                protein_seqs_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_seqs_pickle_path)

                if not protein_pickle.representative_sequence:
                    return

                if len(protein_pickle.sequences) < 1:
                    return

                alignment_dir = op.join(outdir, g_id)
                ssbio.utils.make_dir(alignment_dir)
                protein_pickle.pairwise_align_sequences_to_representative(gapopen=gapopen, gapextend=gapextend,
                                                                          engine=engine, outdir=alignment_dir,
                                                                          parse=parse, force_rerun=force_rerun)
                protein_pickle.save_pickle(outfile=protein_seqs_aln_pickle_path)

            return g_id, protein_seqs_aln_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result_raw = genes_rdd.map(_align_orthologous_gene_pairwise_sc).collect()
        else:
            result_raw = []
            for g in tqdm(g_ids):
                result_raw.append(self._align_orthologous_gene_pairwise(g, gapopen=gapopen, gapextend=gapextend,
                                                                    engine=engine, parse=parse,
                                                                    force_rerun=force_rerun))

        result = [x for x in result_raw if x is not None]
        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
174729
def ref_string_matches_dict_of_ref_sequences(self, ref_sequences):
        ''''''
        return self.CHROM in ref_sequences and self.ref_string_matches_ref_sequence(ref_sequences[self.CHROM])
706195
def seq_info(names, id2names, insertions, sequences):
    
    seqs = {} # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns]], ...]]
    for name in names:
        id = id2names[name]
        gene = name.split('fromHMM::', 1)[0].rsplit(' ', 1)[1]
        model = name.split('fromHMM::', 1)[1].split('=', 1)[1].split()[0]
        i_gene_pos = insertions[id] # coordinates of each insertion wrt gene
        i_model_pos = name.split('fromHMM::', 1)[1].split('model-pos(ins-len)=')[1].split()[0].split(';') # model overlap
        i_info = []
        for i, ins in enumerate(i_gene_pos):
            model_pos = i_model_pos[i].split('-')[1].split('(')[0]
            length = i_model_pos[i].split('(')[1].split(')')[0]
            iheader = '>%s_%s insertion::seq=%s type=insertion strand=n/a gene-pos=%s-%s model-pos=%s'\
                    % (id, (i + 1), (i + 1), ins[0], ins[1], model_pos)
            iseq = sequences[id][1][ins[0]:(ins[1] + 1)]
            iseq = [iheader, iseq]
            info = [ins, model_pos, length, iseq, [], []]
            i_info.append(info)
        seqs[id] = [gene, model, i_info]
    return seqs
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
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
1161304
def get_alignment_strings(self,min_intron_size=68):
    
    qseq = self.query_sequence
    if not qseq:
      sys.exit("ERROR: Query sequence must be accessable to get alignment strings\n")
      sys.exit()
    ref = self._options.reference
    qual = self.query_quality
    if not qual: 
      qual = 'I'*len(qseq) # for a placeholder quality
    if self.strand == '-': 
      qseq = rc(qseq)
      qual = qual[::-1]
    tarr = []
    qarr = []
    yarr = []
    tdone = ''
    qdone = ''
    ydone = '' #query quality
    for i in range(len(self.alignment_ranges)):
      [t,q] = self.alignment_ranges[i]
      textra = ''
      qextra = ''
      yextra = ''
      if i >= 1:
        dift = t.start-self.alignment_ranges[i-1][0].end-1
        difq = q.start-self.alignment_ranges[i-1][1].end-1
        if dift < min_intron_size:
          if dift > 0:
            textra = str(ref[t.chr][t.start-dift-1:t.start-1]).upper()
            qextra = '-'*dift
            yextra = '\0'*dift
          elif difq > 0:
            textra = '-'*difq
            qextra = qseq[q.start-difq-1:q.start-1].upper()
            yextra = qual[q.start-difq-1:q.start-1]
        else:
          tarr.append(tdone)
          qarr.append(qdone)
          yarr.append(ydone)
          tdone = ''
          qdone = ''
          ydone = ''
      tdone += textra+str(ref[t.chr][t.start-1:t.end]).upper()
      qdone += qextra+qseq[q.start-1:q.end].upper()
      ydone += yextra+qual[q.start-1:q.end]
    if len(tdone) > 0: 
      tarr.append(tdone)
      qarr.append(qdone)
      yarr.append(ydone)
    if self.query_quality == '*': yarr = [x.replace('I',' ') for x in yarr]
    #query, target, query_quality
    return [qarr,tarr,yarr]
1371184
def _guess_sequence_type_from_string(self, seq):
        ''''''
        # Define expected residues for each sequence type
        aa_chars = ['P','V','L','I','M','F','Y','W','H','K','R','Q','N','E','D','S','X','*']
        aas = set(itertools.chain(aa_chars, [lower(a) for a in aa_chars]))
        na_chars = ['A', 'T', 'G', 'C', 'N', 'U']
        nas = set(itertools.chain(na_chars, [lower(a) for a in na_chars]))

        num_nucleotide = 0
        num_protein = 0
        count = 0
        for residue in seq:
            if residue in nas:
                num_nucleotide += 1
            elif residue in aas:
                num_protein += 1
            else:
                raise Exception('Encountered unexpected character when attempting to guess sequence type: %s' % (residue))
            count += 1
            if count >300: break
        if float(num_protein) / (num_protein+num_nucleotide) > 0.1:
            return self.PROTEIN_SEQUENCE_TYPE
        else:
            return self.NUCLEOTIDE_SEQUENCE_TYPE
437650
def inner_join_impl(other, sequence):
    
    seq_dict = {}
    for element in sequence:
        seq_dict[element[0]] = element[1]
    seq_kv = seq_dict
    other_kv = dict(other)
    keys = seq_kv.keys() if len(seq_kv) < len(other_kv) else other_kv.keys()
    result = {}
    for k in keys:
        if k in seq_kv and k in other_kv:
            result[k] = (seq_kv[k], other_kv[k])
    return six.viewitems(result)
21460
def analyze_insertions(fa, threads = 6):
    
    safe, sequences, id2name, names, insertions = analyze_fa(fa)
    seqs = seq_info(names, id2name, insertions, sequences)
    seqs, orfs = find_orfs(safe, seqs)
    seqs = find_introns(safe, seqs, sequences, threads)
    seqs = seqs2bool(seqs)
    seqs = annotate_orfs(orfs, seqs, threads)
    return seqs, id2name
923259
def get_generic_subseq_2D(protein, cutoff, prop, condition):
        """"""
        subseq, subseq_resnums = protein.representative_sequence.get_subsequence_from_property(property_key=prop,
                                                                                               property_value=cutoff,
                                                                                               condition=condition,
                                                                                               return_resnums=True) or (
                                     None, [])

        return {'subseq_len': len(subseq_resnums), 'subseq': subseq, 'subseq_resnums': subseq_resnums}