Use aligned prot_seq to align nucl_seq
12299
def map_alignment(self, prot_seq, nucl_seq):
        
        if prot_seq.id != nucl_seq.id:
            logging.warning(
                'ID mismatch: %s != %s. Are the sequences in the same order?',
                prot_seq.id, nucl_seq.id)

        # Ungap nucleotides
        codons = batch(str(nucl_seq.seq.ungap('-')), 3)
        codons = [''.join(i) for i in codons]
        codon_iter = iter(codons)

        ungapped_prot = str(prot_seq.seq).replace('-', '')

        if len(ungapped_prot) != len(codons):
            table = self.translation_table.forward_table
            prot_str = ' '.join(' ' + p + ' ' for p in ungapped_prot)
            codon_str = ' '.join(codons)
            trans_str = ' '.join(' ' + table.get(codon, 'X') + ' '
                                 for codon in codons)
            raise ValueError("""Length of codon sequence ({0}) does not match \
length of protein sequence ({1}) for {2}
Protein:       {3}
Codons:        {4}
Trans. Codons: {5}""".format(len(codons), len(ungapped_prot), nucl_seq.id, prot_str,
                            codon_str, trans_str))

        try:
            nucl_align = ['---' if p == '-' else next(codon_iter)
                          for p in str(prot_seq.seq)]
        except StopIteration:
            assert False  # Should be checked above

        result = SeqRecord(Seq(''.join(nucl_align)), id=nucl_seq.id,
                           description=nucl_seq.description)

        # Validate
        self._validate_translation(prot_seq.seq.upper(), result.seq.upper())

        return result
832792
def align_pipe(fastq_file, pair_file, ref_file, names, align_dir, data):
    
    pair_file = pair_file if pair_file else ""
    # back compatible -- older files were named with lane information, use sample name now
    out_file = os.path.join(align_dir, "{0}-sort.bam".format(names["lane"]))
    if not utils.file_exists(out_file):
        out_file = os.path.join(align_dir, "{0}-sort.bam".format(dd.get_sample_name(data)))
    if data.get("align_split") or fastq_file.endswith(".sdf"):
        final_file = out_file
        out_file, data = alignprep.setup_combine(final_file, data)
        fastq_file, pair_file = alignprep.split_namedpipe_cls(fastq_file, pair_file, data)
    else:
        final_file = None
    samtools = config_utils.get_program("samtools", data["config"])
    novoalign = config_utils.get_program("novoalign", data["config"])
    resources = config_utils.get_resources("novoalign", data["config"])
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    max_mem = resources.get("memory", "1G")
    extra_novo_args = " ".join(_novoalign_args_from_config(data["config"]))
    rg_info = get_rg_info(names)
    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with tx_tmpdir(data) as work_dir:
            with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
                tx_out_prefix = os.path.splitext(tx_out_file)[0]
                cmd = ("unset JAVA_HOME && "
                       "{novoalign} -o SAM '{rg_info}' -d {ref_file} -f {fastq_file} {pair_file} "
                       "  -c {num_cores} {extra_novo_args} | ")
                cmd = (cmd + tobam_cl).format(**locals())
                do.run(cmd, "Novoalign: %s" % names["sample"], None,
                       [do.file_nonempty(tx_out_file), do.file_reasonable_size(tx_out_file, fastq_file)])
    data["work_bam"] = out_file
    return data
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
1208101
def _nucmer_hits_to_assembly_coords(nucmer_hits):
        ''''''
        coords = {}
        for l in nucmer_hits.values():
            for hit in l:
                if hit.qry_name not in coords:
                    coords[hit.qry_name] = []
                coords[hit.qry_name].append(hit.qry_coords())

        for scaff in coords:
            pyfastaq.intervals.merge_overlapping_in_list(coords[scaff])

        return coords
795201
def _align_mem(fastq_file, pair_file, ref_file, out_file, names, rg_info, data):
    
    with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
        cmd = ("unset JAVA_HOME && "
               "%s | %s" % (_get_bwa_mem_cmd(data, out_file, ref_file, fastq_file, pair_file), tobam_cl))
        do.run(cmd, "bwa mem alignment from fastq: %s" % names["sample"], None,
                [do.file_nonempty(tx_out_file), do.file_reasonable_size(tx_out_file, fastq_file)])
    return out_file
395665
def align(self, alignment_tool = 'clustalw', gap_opening_penalty = 0.2, ignore_bad_chains = False):
        ''''''
        if len(self.pdb_chains) > 1:
            sa = SequenceAligner(alignment_tool = alignment_tool, gap_opening_penalty = gap_opening_penalty)
            for pdb_chain in self.pdb_chains:
                sa.add_sequence('%s_%s' % (pdb_chain['pdb_id'], pdb_chain['chain_id']), pdb_chain['sequence'], ignore_bad_chains = ignore_bad_chains)
            best_matches = sa.align()
            return sa.alignment_output, best_matches
        else:
            raise Exception('Cannot align sequences - less than two chains were specified.')
71401
def gatk_indel_realignment_cl(runner, align_bam, ref_file, intervals,
                              tmp_dir, region=None, deep_coverage=False,
                              known_vrns=None):
    
    if not known_vrns:
        known_vrns = {}
    params = ["-T", "IndelRealigner",
              "-I", align_bam,
              "-R", ref_file,
              "-targetIntervals", intervals,
              ]
    if region:
        params += ["-L", region]
    if known_vrns.get("train_indels"):
        params += ["--knownAlleles", known_vrns["train_indels"]]
    if deep_coverage:
        params += ["--maxReadsInMemory", "300000",
                   "--maxReadsForRealignment", str(int(5e5)),
                   "--maxReadsForConsensuses", "500",
                   "--maxConsensuses", "100"]
    return runner.cl_gatk(params, tmp_dir)
125231
def alignUnalignedSeqsFromFile(filename,WorkingDir=None,SuppressStderr=None,\
    SuppressStdout=None):
    
    app = Clustalw(WorkingDir=WorkingDir,SuppressStderr=SuppressStderr,\
        SuppressStdout=SuppressStdout)
    return app(filename)
603006
def _balanced_strand_gen(base_aligns, limit):
    ''''''
    predicate=lambda align: align.is_reverse
    gen1, gen2 = itertools.tee((predicate(align), align) for align in base_aligns)
    for_strand_gen = (align for pred, align in gen1 if not pred)
    rev_strand_gen = (align for pred, align in gen2 if pred)
    combined_gen = iter_zip(for_strand_gen, rev_strand_gen)
    interleaved_gen = (align for for_rev_aligns in combined_gen for align in for_rev_aligns)
    return itertools.islice(interleaved_gen, limit)
348522
def postprocess_alignment(data):
    
    data = cwlutils.normalize_missing(utils.to_single_data(data))
    data = cwlutils.unpack_tarballs(data, data)
    bam_file = data.get("align_bam") or data.get("work_bam")
    ref_file = dd.get_ref_file(data)
    if vmulti.bam_needs_processing(data) and bam_file and bam_file.endswith(".bam"):
        out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                  dd.get_sample_name(data)))
        bam_file_ready = os.path.join(out_dir, os.path.basename(bam_file))
        if not utils.file_exists(bam_file_ready):
            utils.symlink_plus(bam_file, bam_file_ready)
        bam.index(bam_file_ready, data["config"])
        covinfo = callable.sample_callable_bed(bam_file_ready, ref_file, data)
        callable_region_bed, nblock_bed = \
            callable.block_regions(covinfo.raw_callable, bam_file_ready, ref_file, data)
        data["regions"] = {"nblock": nblock_bed,
                           "callable": covinfo.raw_callable,
                           "sample_callable": covinfo.callable,
                           "mapped_stats": readstats.get_cache_file(data)}
        data["depth"] = covinfo.depth_files
        data = coverage.assign_interval(data)
        data = samtools.run_and_save(data)
        data = recalibrate.prep_recal(data)
        data = recalibrate.apply_recal(data)
    elif dd.get_variant_regions(data):
        callable_region_bed, nblock_bed = \
            callable.block_regions(dd.get_variant_regions(data), bam_file, ref_file, data)
        data["regions"] = {"nblock": nblock_bed, "callable": dd.get_variant_regions(data),
                           "sample_callable": dd.get_variant_regions(data)}
    return [[data]]
1053237
def align(s1,s2,test=False):
    
    import operator
    from Bio import pairwise2
    alignments = pairwise2.align.localms(s1.upper(),s2.upper(),1,-1,-5,-5)
    if test:
        print(alignments)
    alignsymb=np.nan
    score=np.nan
    sorted_alignments = sorted(alignments, key=operator.itemgetter(2))
    for a in alignments:
        alignstr=pairwise2.format_alignment(*a)
        alignsymb=alignstr.split('\n')[1]
        score=a[2]
        if test:
            print(alignstr)
        break
    return alignsymb,score
1357336
def _align_backtrack(fastq_file, pair_file, ref_file, out_file, names, rg_info, data):
    
    bwa = config_utils.get_program("bwa", data["config"])
    config = data["config"]
    sai1_file = "%s_1.sai" % os.path.splitext(out_file)[0]
    sai2_file = "%s_2.sai" % os.path.splitext(out_file)[0] if pair_file else ""
    if not utils.file_exists(sai1_file):
        with file_transaction(data, sai1_file) as tx_sai1_file:
            _run_bwa_align(fastq_file, ref_file, tx_sai1_file, config)
    if sai2_file and not utils.file_exists(sai2_file):
        with file_transaction(data, sai2_file) as tx_sai2_file:
            _run_bwa_align(pair_file, ref_file, tx_sai2_file, config)
    with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
        align_type = "sampe" if sai2_file else "samse"
        cmd = ("unset JAVA_HOME && {bwa} {align_type} -r '{rg_info}' {ref_file} {sai1_file} {sai2_file} "
               "{fastq_file} {pair_file} | ")
        cmd = cmd.format(**locals()) + tobam_cl
        do.run(cmd, "bwa %s" % align_type, data)
    return out_file
1071843
def subalignment(alnfle, subtype, alntype="fasta"):
    
    aln = AlignIO.read(alnfle, alntype)
    alnlen = aln.get_alignment_length()
    nseq = len(aln)
    subaln = None
    subalnfile = alnfle.rsplit(".", 1)[0] + "_{0}.{1}".format(subtype, alntype)

    if subtype == "synonymous":
        for j in range( 0, alnlen, 3 ):
            aa = None
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in CODON_TRANSLATION:
                    break
                if aa and CODON_TRANSLATION[codon] != aa:
                    break
                else:
                    aa = CODON_TRANSLATION[codon]
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subtype == "fourfold":
        for j in range( 0, alnlen, 3 ):
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in FOURFOLD:
                    break
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subaln:
        AlignIO.write(subaln, subalnfile, alntype)
        return subalnfile
    else:
        print("No sites {0} selected.".format(subtype), file=sys.stderr)
        return None
1015431
def align_pipe(fastq_file, pair_file, ref_file, names, align_dir, data):
    
    pair_file = pair_file if pair_file else ""
    # back compatible -- older files were named with lane information, use sample name now
    if names["lane"] != dd.get_sample_name(data):
        out_file = os.path.join(align_dir, "{0}-sort.bam".format(names["lane"]))
    else:
        out_file = None
    if not out_file or not utils.file_exists(out_file):
        umi_ext = "-cumi" if "umi_bam" in data else ""
        out_file = os.path.join(align_dir, "{0}-sort{1}.bam".format(dd.get_sample_name(data), umi_ext))
    qual_format = data["config"]["algorithm"].get("quality_format", "").lower()
    min_size = None
    if data.get("align_split") or fastq_file.endswith(".sdf"):
        if fastq_file.endswith(".sdf"):
            min_size = rtg.min_read_size(fastq_file)
        final_file = out_file
        out_file, data = alignprep.setup_combine(final_file, data)
        fastq_file, pair_file = alignprep.split_namedpipe_cls(fastq_file, pair_file, data)
    else:
        final_file = None
        if qual_format == "illumina":
            fastq_file = alignprep.fastq_convert_pipe_cl(fastq_file, data)
            if pair_file:
                pair_file = alignprep.fastq_convert_pipe_cl(pair_file, data)
    rg_info = novoalign.get_rg_info(names)
    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        # If we cannot do piping, use older bwa aln approach
        if ("bwa-mem" not in dd.get_tools_on(data) and
              ("bwa-mem" in dd.get_tools_off(data) or not _can_use_mem(fastq_file, data, min_size))):
            out_file = _align_backtrack(fastq_file, pair_file, ref_file, out_file,
                                        names, rg_info, data)
        else:
            out_file = _align_mem(fastq_file, pair_file, ref_file, out_file,
                                  names, rg_info, data)
    data["work_bam"] = out_file
    return data
100763
def libseq_parse(self, line):
        """"""
        count = line.count(',')
        if count > 1:
            raise ValueError("""ERROR: There are too many commas in your
            lib_seq line. There should only be one, between the two filepath
            globs.""")
        if count == 0:
            raise ValueError("""ERROR: There are no commas in your
            lib_seq line. There should be one, between the two filepath
            globs.""")
        line = [x.strip() for x in line.split() if x]
        #print(line)

        #This block tries to fix user input if there is an inappropriate
        # number of commas or spaces associated with commas
        if len(line) > 13:
            commas = [x for x in line if "," in x]
            if commas:
                #this handles the case if there was an extra space added on
                # both sides of the comma
                if (len(commas) == 1) and (commas[0] != ','):
                    if (commas[0][-1] == ',') or (commas[0][0] == ','):
                        line = [line[0]] + ["".join(line[1:3])] + line[3:]
                #This handles the case where there is an extra space on either
                #side of the comma
                elif line[2] == ',':
                    line = [line[0]] + ["".join(line[1:4])] + line[4::]

        #make sure there is an appropriate numer of elements in list
        if (len(line) ) != 13:
            raise ValueError("""ERROR: There are {0} values in the lib_seq line when
                             there should be 12. Check your input meraculous
                             config file for errors and make sure that it matches
                             the specification for the Meraculous manual.""".format(
                                 len(line)))

        for index in self.indices:
            if index == 1:
                globs = [os.path.abspath(x.strip()) for x in line[1].split(',') if x]
                #make sure there are only two filepaths to search for globs
                if len(globs) != 2:
                    raise ValueError("""ERROR: remember to split your glob
                    filepaths in lib_seq with a comma. This error occured because
                    the program was looking for two glob filepaths, but instead
                    found {}: {}""".format(len(globs), globs))
                self.data["globs"] = globs
                #Note that glob.glob returns an empty list if the file or glob
                # filepath does not exist.
                forwards = glob.glob(os.path.abspath(globs[0]))
                reverses = glob.glob(os.path.abspath(globs[1]))
                # print(os.path.abspath(globs[0]))
                # print("globs: ", globs)
                # print("forwards: ", forwards)
                # print("reverses: ", reverses)
                #At this point the
                if (len(forwards) < 1) or (len(reverses) < 1):
                    # print("forwards: ", forwards)
                    # print("reverses: ", reverses)

                    report = []
                    if len(forwards) < 1:
                        report.append(globs[0])
                    if len(reverses) < 1:
                        report.append(globs[1])
                    raise ValueError("""ERROR: no read files were found for the glob strings,
                    {}.

                    Please check the glob string and/or use absolute
                    file paths. A common source of this error is using
                    relative filepaths from the incorrect directory, or
                    referring to files that do not exist.""".format(report))
                if len(forwards) != len(reverses):
                    raise ValueError("""ERROR: the number of reverse read files does not
                    match the number of forward read files""")
                #zip the pairs together into tuples, not lists
                self.data["pairs"] = [x for x in zip(sorted(forwards), sorted(reverses))]

            self.data[self.indices.get(index)] = line[index]
57053
def _align_orthologous_gene_pairwise(self, g_id, gapopen=10, gapextend=0.5, engine='needle', parse=True, force_rerun=False):
        """"""
        protein_seqs_aln_pickle_path = op.join(self.sequences_by_gene_dir, '{}_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 = self.gene_protein_pickles[g_id]
            protein_pickle = ssbio.io.load_pickle(protein_seqs_pickle_path)

            if not protein_pickle.representative_sequence:
                log.error('{}: no representative sequence to align to'.format(g_id))
                return

            if len(protein_pickle.sequences) < 1:
                log.error('{}: no other sequences to align to'.format(g_id))
                return

            alignment_dir = op.join(self.sequences_by_gene_dir, 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
959604
def align(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(align.__doc__)
    p.set_firstN(firstN=0)
    p.add_option("--full", default=False, action="store_true",
                 help="Enforce end-to-end alignment [default: local]")
    p.add_option("--reorder", default=False, action="store_true",
                 help="Keep the input read order [default: %default]")
    p.add_option("--null", default=False, action="store_true",
                 help="Do not write to SAM/BAM output")
    p.add_option("--fasta", default=False, action="store_true",
                 help="Query reads are FASTA")
    p.set_cutoff(cutoff=800)
    p.set_mateorientation(mateorientation="+-")
    p.set_sam_options(bowtie=True)

    opts, args = p.parse_args(args)
    extra = opts.extra
    mo = opts.mateorientation
    if mo == '+-':
        extra += ""
    elif mo == '-+':
        extra += "--rf"
    else:
        extra += "--ff"

    PE = True
    if len(args) == 2:
        logging.debug("Single-end alignment")
        PE = False
    elif len(args) == 3:
        logging.debug("Paired-end alignment")
    else:
        sys.exit(not p.print_help())

    firstN = opts.firstN
    mapped = opts.mapped
    unmapped = opts.unmapped
    fasta = opts.fasta
    gl = "--end-to-end" if opts.full else "--local"

    dbfile, readfile = args[0:2]
    dbfile = check_index(dbfile)
    prefix = get_prefix(readfile, dbfile)
    samfile, mapped, unmapped = get_samfile(readfile, dbfile, bowtie=True,
                                            mapped=mapped, unmapped=unmapped,
                                            bam=opts.bam)
    logfile = prefix + ".log"
    if not fasta:
        offset = guessoffset([readfile])

    if not need_update(dbfile, samfile):
        logging.error("`{0}` exists. `bowtie2` already run.".format(samfile))
        return samfile, logfile

    cmd = "bowtie2 -x {0}".format(dbfile)
    if PE:
        r1, r2 = args[1:3]
        cmd += " -1 {0} -2 {1}".format(r1, r2)
        cmd += " --maxins {0}".format(opts.cutoff)
        mtag, utag = "--al-conc", "--un-conc"
    else:
        cmd += " -U {0}".format(readfile)
        mtag, utag = "--al", "--un"

    if mapped:
        cmd += " {0} {1}".format(mtag, mapped)
    if unmapped:
        cmd += " {0} {1}".format(utag, unmapped)

    if firstN:
        cmd += " --upto {0}".format(firstN)
    cmd += " -p {0}".format(opts.cpus)
    if fasta:
        cmd += " -f"
    else:
        cmd += " --phred{0}".format(offset)
    cmd += " {0}".format(gl)
    if opts.reorder:
        cmd += " --reorder"

    cmd += " {0}".format(extra)
    # Finally the log
    cmd += " 2> {0}".format(logfile)

    if opts.null:
        samfile = "/dev/null"

    cmd = output_bam(cmd, samfile)
    sh(cmd)
    print(open(logfile).read(), file=sys.stderr)

    return samfile, logfile