Aligns the N-CA and CA-CB vector of the target monomer.

    Parameters
    ----------
    tar: ampal.Residue
        The residue that will be aligned to the reference.
    ref: ampal.Residue
        The reference residue for the alignment.
964934
def align_nab(tar, ref):
    
    rot_trans_1 = find_transformations(
        tar['N'].array, tar['CA'].array, ref['N'].array, ref['CA'].array)
    apply_trans_rot(tar, *rot_trans_1)
    rot_ang_ca_cb = dihedral(tar['CB'], ref['CA'], ref['N'], ref['CB'])
    tar.rotate(rot_ang_ca_cb, ref['N'].array - ref['CA'].array, ref['N'].array)
    return
647909
def pistacking(rings_bs, rings_lig):
    """"""
    data = namedtuple(
        'pistack', 'proteinring ligandring distance angle offset type restype resnr reschain restype_l resnr_l reschain_l')
    pairings = []
    for r, l in itertools.product(rings_bs, rings_lig):
        # DISTANCE AND RING ANGLE CALCULATION
        d = euclidean3d(r.center, l.center)
        b = vecangle(r.normal, l.normal)
        a = min(b, 180 - b if not 180 - b < 0 else b)  # Smallest of two angles, depending on direction of normal

        # RING CENTER OFFSET CALCULATION (project each ring center into the other ring)
        proj1 = projection(l.normal, l.center, r.center)
        proj2 = projection(r.normal, r.center, l.center)
        offset = min(euclidean3d(proj1, l.center), euclidean3d(proj2, r.center))

        # RECEPTOR DATA
        resnr, restype, reschain = whichresnumber(r.atoms[0]), whichrestype(r.atoms[0]), whichchain(r.atoms[0])
        resnr_l, restype_l, reschain_l = whichresnumber(l.orig_atoms[0]), whichrestype(
            l.orig_atoms[0]), whichchain(l.orig_atoms[0])

        # SELECTION BY DISTANCE, ANGLE AND OFFSET
        passed = False
        if not config.MIN_DIST < d < config.PISTACK_DIST_MAX:
            continue
        if 0 < a < config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
            ptype = 'P'
            passed = True
        if 90 - config.PISTACK_ANG_DEV < a < 90 + config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
            ptype = 'T'
            passed = True
        if passed:
            contact = data(proteinring=r, ligandring=l, distance=d, angle=a, offset=offset,
                           type=ptype, resnr=resnr, restype=restype, reschain=reschain,
                           resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l)
            pairings.append(contact)
    return filter_contacts(pairings)
702028
public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
			throws StructureException {

		reset();
		this.params = params;

		long timeStart = System.currentTimeMillis();

		// step 1 get all Diagonals of length X that are similar between both
		// structures
		logger.debug(" length atoms1:" + ca1.length);
		logger.debug(" length atoms2:" + ca2.length);

		logger.debug("step 1 - get fragments with similar intramolecular distances ");

		int k = params.getDiagonalDistance();
		int k2 = params.getDiagonalDistance2();
		int fragmentLength = params.getFragmentLength();

		if (ca1.length < (fragmentLength + 1)) {
			throw new StructureException("structure 1 too short (" + ca1.length
					+ "), can not align");
		}
		if (ca2.length < (fragmentLength + 1)) {
			throw new StructureException("structure 2 too short (" + ca2.length
					+ "), can not align");
		}
		int rows = ca1.length - fragmentLength + 1;
		int cols = ca2.length - fragmentLength + 1;
		distanceMatrix = new Matrix(rows, cols, 0.0);

		double[] dist1 = AlignUtils.getDiagonalAtK(ca1, k);

		double[] dist2 = AlignUtils.getDiagonalAtK(ca2, k);
		double[] dist3 = new double[0];
		double[] dist4 = new double[0];
		if (k2 > 0) {
			dist3 = AlignUtils.getDiagonalAtK(ca1, k2);
			dist4 = AlignUtils.getDiagonalAtK(ca2, k2);
		}

		double[][] utmp = new double[][] { { 0, 0, 1 } };
		Atom unitvector = new AtomImpl();
		unitvector.setCoords(utmp[0]);

		List<FragmentPair> fragments = new ArrayList<FragmentPair>();

		for (int i = 0; i < rows; i++) {

			Atom[] catmp1 = AlignUtils.getFragment(ca1, i, fragmentLength);
			Atom center1 = AlignUtils.getCenter(ca1, i, fragmentLength);

			for (int j = 0; j < cols; j++) {

				double rdd1 = AlignUtils.rms_dk_diag(dist1, dist2, i, j,
						fragmentLength, k);
				double rdd2 = 0;
				if (k2 > 0)
					rdd2 = AlignUtils.rms_dk_diag(dist3, dist4, i, j,
							fragmentLength, k2);
				double rdd = rdd1 + rdd2;
				distanceMatrix.set(i, j, rdd);

				if (rdd < params.getFragmentMiniDistance()) {
					FragmentPair f = new FragmentPair(fragmentLength, i, j);
					Atom[] catmp2 = AlignUtils.getFragment(ca2, j,
							fragmentLength);
					Atom center2 = AlignUtils.getCenter(ca2, j,
							fragmentLength);

					f.setCenter1(center1);
					f.setCenter2(center2);

					Matrix4d t = SuperPositions.superpose(
							Calc.atomsToPoints(catmp1),
							Calc.atomsToPoints(catmp2));

					Matrix rotmat = Matrices.getRotationJAMA(t);
					f.setRot(rotmat);

					Atom aunitv = (Atom) unitvector.clone();
					Calc.rotate(aunitv, rotmat);
					f.setUnitv(aunitv);

					boolean doNotAdd = false;
					if (params.reduceInitialFragments()) {
						doNotAdd = FragmentJoiner.reduceFragments(
								fragments, f, distanceMatrix);

					}
					if (doNotAdd)
						continue;

					fragments.add(f);
				}
			}
		}

		notifyFragmentListeners(fragments);

		FragmentPair[] fp = fragments
				.toArray(new FragmentPair[fragments.size()]);
		setFragmentPairs(fp);

		logger.debug(" got # fragment pairs: {}", fp.length);

		logger.debug("step 2 - join fragments");

		// step 2 combine them to possible models
		FragmentJoiner joiner = new FragmentJoiner();

		JointFragments[] frags;

		if (params.isJoinFast()) {
			// apply the quick alignment procedure.
			// less quality in alignments, better for DB searches...
			frags = joiner.approach_ap3(ca1, ca2, fp, params);

			joiner.extendFragments(ca1, ca2, frags, params);

		} else if (params.isJoinPlo()) {
			// this approach by StrComPy (peter lackner):
			frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
					params.getFragCompat(), params.getMaxrefine());

		} else {

			// my first implementation
			frags = joiner.approach_ap3(ca1, ca2, fp, params);
		}

		notifyJointFragments(frags);

		logger.debug(" number joint fragments: ", frags.length);

		logger.debug("step 3 - refine alignments");

		List<AlternativeAlignment> aas = new ArrayList<AlternativeAlignment>();
		for (int i = 0; i < frags.length; i++) {
			JointFragments f = frags[i];
			AlternativeAlignment a = new AlternativeAlignment();
			a.apairs_from_idxlst(f);
			a.setAltAligNumber(i + 1);
			a.setDistanceMatrix(distanceMatrix);

			try {
				if (params.getMaxIter() > 0) {

					a.refine(params, ca1, ca2);
				} else {

					a.finish(params, ca1, ca2);

				}
			} catch (StructureException e) {
				logger.error("Refinement of fragment {} failed", i, e);
			}
			a.calcScores(ca1, ca2);
			aas.add(a);
		}

		// sort the alternative alignments
		Comparator<AlternativeAlignment> comp = new AltAligComparator();
		Collections.sort(aas, comp);
		Collections.reverse(aas);

		alts = aas.toArray(new AlternativeAlignment[aas.size()]);
		// do final numbering of alternative solutions
		int aanbr = 0;
		for (AlternativeAlignment a : alts) {
			aanbr++;
			a.setAltAligNumber(aanbr);
		}

		logger.debug("total calculation time: {} ms.",
				(System.currentTimeMillis() - timeStart));
	}
200715
private void rotateShiftAtoms(Atom[] ca){

		for (int i  = 0 ; i < ca.length; i++){
			Atom c = ca[i];

			Calc.rotate(c,currentRotMatrix);
			Calc.shift(c,currentTranMatrix);
			//System.out.println("after " + c);
			ca[i] = c;
		}
		//System.out.println("after " + ca[0]);
	}
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]]
680495
def makeCys(segID, N, CA, C, O, geo):
    ''''''
    ##R-Group
    CA_CB_length=geo.CA_CB_length
    C_CA_CB_angle=geo.C_CA_CB_angle
    N_C_CA_CB_diangle=geo.N_C_CA_CB_diangle
    
    CB_SG_length= geo.CB_SG_length
    CA_CB_SG_angle= geo.CA_CB_SG_angle
    N_CA_CB_SG_diangle= geo.N_CA_CB_SG_diangle

    carbon_b= calculateCoordinates(N, C, CA, CA_CB_length, C_CA_CB_angle, N_C_CA_CB_diangle)
    CB= Atom("CB", carbon_b, 0.0 , 1.0, " "," CB", 0,"C")
    sulfur_g= calculateCoordinates(N, CA, CB, CB_SG_length, CA_CB_SG_angle, N_CA_CB_SG_diangle)
    SG= Atom("SG", sulfur_g, 0.0, 1.0, " ", " SG", 0, "S")

    ##Create Residue Data Structure
    res= Residue((' ', segID, ' '), "CYS", '    ')
    res.add(N)
    res.add(CA)
    res.add(C)
    res.add(O)
    res.add(CB)
    res.add(SG)
    return res
454444
def bytealign(self):
        
        skipped = (8 - (self._pos % 8)) % 8
        self.pos += self._offset + skipped
        assert self._assertsanity()
        return skipped
15689
def check_for_non_canonical(residue):
        """"""
        res_label = list(residue[0])[0][2]
        atom_labels = {x[2] for x in itertools.chain(
            *residue[1].values())}  # Used to find unnatural aas
        if (all(x in atom_labels for x in ['N', 'CA', 'C', 'O'])) and (
                len(res_label) == 3):
            return Residue, True
        return None
1167499
def CA_TRACE(residue):
    

    cg_residue = deepcopy(residue)

    for atom in cg_residue.child_dict.keys():
        if cg_residue[atom].name != 'CA':
            cg_residue.detach_child(cg_residue[atom].name)
    
    return cg_residue
382277
def align_bam(in_bam, ref_file, names, align_dir, data):
    
    config = data["config"]
    out_file = os.path.join(align_dir, "{0}-sort.bam".format(names["lane"]))
    novoalign = config_utils.get_program("novoalign", config)
    samtools = config_utils.get_program("samtools", config)
    resources = config_utils.get_resources("novoalign", config)
    num_cores = config["algorithm"].get("num_cores", 1)
    max_mem = resources.get("memory", "4G").upper()
    extra_novo_args = " ".join(_novoalign_args_from_config(config, False))

    if not file_exists(out_file):
        with tx_tmpdir(data, base_dir=align_dir) as work_dir:
            with postalign.tobam_cl(data, out_file, bam.is_paired(in_bam)) as (tobam_cl, tx_out_file):
                rg_info = get_rg_info(names)
                tx_out_prefix = os.path.splitext(tx_out_file)[0]
                prefix1 = "%s-in1" % tx_out_prefix
                cmd = ("unset JAVA_HOME && "
                       "{samtools} sort -n -o -l 1 -@ {num_cores} -m {max_mem} {in_bam} {prefix1} "
                       "| {novoalign} -o SAM '{rg_info}' -d {ref_file} -f /dev/stdin "
                       "  -F BAMPE -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, in_bam)])
    return out_file
187976
def process_alignment(data, alt_input=None):
    
    data = cwlutils.normalize_missing(utils.to_single_data(data))
    data = cwlutils.unpack_tarballs(data, data)
    fastq1, fastq2 = dd.get_input_sequence_files(data)
    if alt_input:
        fastq1, fastq2 = alt_input
    config = data["config"]
    aligner = config["algorithm"].get("aligner", None)
    if fastq1 and objectstore.file_exists_or_remote(fastq1) and aligner:
        logger.info("Aligning lane %s with %s aligner" % (data["rgnames"]["lane"], aligner))
        data = align_to_sort_bam(fastq1, fastq2, aligner, data)
        if dd.get_umi_consensus(data):
            data["umi_bam"] = dd.get_work_bam(data)
            if fastq2:
                f1, f2, avg_cov = postalign.umi_consensus(data)
                data["config"]["algorithm"]["rawumi_avg_cov"] = avg_cov
                del data["config"]["algorithm"]["umi_type"]
                data["config"]["algorithm"]["mark_duplicates"] = False
                data = align_to_sort_bam(f1, f2, aligner, data)
            else:
                raise ValueError("Single fastq input for UMI processing; fgbio needs paired reads: %s" %
                                 dd.get_sample_name(data))
        data = _add_supplemental_bams(data)
    elif fastq1 and objectstore.file_exists_or_remote(fastq1) and fastq1.endswith(".bam"):
        sort_method = config["algorithm"].get("bam_sort")
        bamclean = config["algorithm"].get("bam_clean")
        if bamclean is True or bamclean == "picard":
            if sort_method and sort_method != "coordinate":
                raise ValueError("Cannot specify `bam_clean: picard` with `bam_sort` other than coordinate: %s"
                                 % sort_method)
            out_bam = cleanbam.picard_prep(fastq1, data["rgnames"], dd.get_ref_file(data), data["dirs"],
                                           data)
        elif bamclean == "fixrg":
            out_bam = cleanbam.fixrg(fastq1, data["rgnames"], dd.get_ref_file(data), data["dirs"], data)
        elif bamclean == "remove_extracontigs":
            out_bam = cleanbam.remove_extracontigs(fastq1, data)
        elif sort_method:
            runner = broad.runner_from_path("picard", config)
            out_file = os.path.join(data["dirs"]["work"], "{}-sort.bam".format(
                os.path.splitext(os.path.basename(fastq1))[0]))
            if not utils.file_exists(out_file):
                work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "bamclean",
                                                           dd.get_sample_name(data)))
                out_file = os.path.join(work_dir, "{}-sort.bam".format(dd.get_sample_name(data)))
            out_bam = runner.run_fn("picard_sort", fastq1, sort_method, out_file)
        else:
            out_bam = _link_bam_file(fastq1, os.path.join(dd.get_work_dir(data), "prealign",
                                                          dd.get_sample_name(data)), data)
        bam.index(out_bam, data["config"])
        bam.check_header(out_bam, data["rgnames"], dd.get_ref_file(data), data["config"])
        dedup_bam = postalign.dedup_bam(out_bam, data)
        bam.index(dedup_bam, data["config"])
        data["work_bam"] = dedup_bam
    elif fastq1 and objectstore.file_exists_or_remote(fastq1) and fastq1.endswith(".cram"):
        data["work_bam"] = fastq1
    elif fastq1 is None and not dd.get_aligner(data):
        data["config"]["algorithm"]["variantcaller"] = False
        data["work_bam"] = None
    elif not fastq1:
        raise ValueError("No 'files' specified for input sample: %s" % dd.get_sample_name(data))
    elif "kraken" in config["algorithm"]:  # kraken doesn's need bam
        pass
    else:
        raise ValueError("Could not process input file from sample configuration. \n" +
                         fastq1 +
                         "\nIs the path to the file correct or is empty?\n" +
                         "If it is a fastq file (not pre-aligned BAM or CRAM), "
                         "is an aligner specified in the input configuration?")
    if data.get("work_bam"):
        # Add stable 'align_bam' target to use for retrieving raw alignment
        data["align_bam"] = data["work_bam"]
        data = _add_hla_files(data)
    return [[data]]
655892
def makeArg(segID, N, CA, C, O, geo):
    ''''''
    ##R-Group
    CA_CB_length=geo.CA_CB_length
    C_CA_CB_angle=geo.C_CA_CB_angle
    N_C_CA_CB_diangle=geo.N_C_CA_CB_diangle

    CB_CG_length=geo.CB_CG_length
    CA_CB_CG_angle= geo.CA_CB_CG_angle
    N_CA_CB_CG_diangle=geo.N_CA_CB_CG_diangle
    
    CG_CD_length=geo.CG_CD_length
    CB_CG_CD_angle=geo.CB_CG_CD_angle
    CA_CB_CG_CD_diangle=geo.CA_CB_CG_CD_diangle
    
    CD_NE_length=geo.CD_NE_length
    CG_CD_NE_angle=geo.CG_CD_NE_angle
    CB_CG_CD_NE_diangle=geo.CB_CG_CD_NE_diangle

    NE_CZ_length=geo.NE_CZ_length
    CD_NE_CZ_angle=geo.CD_NE_CZ_angle
    CG_CD_NE_CZ_diangle=geo.CG_CD_NE_CZ_diangle

    CZ_NH1_length=geo.CZ_NH1_length
    NE_CZ_NH1_angle=geo.NE_CZ_NH1_angle
    CD_NE_CZ_NH1_diangle=geo.CD_NE_CZ_NH1_diangle

    CZ_NH2_length=geo.CZ_NH2_length
    NE_CZ_NH2_angle=geo.NE_CZ_NH2_angle
    CD_NE_CZ_NH2_diangle=geo.CD_NE_CZ_NH2_diangle
    
    carbon_b= calculateCoordinates(N, C, CA, CA_CB_length, C_CA_CB_angle, N_C_CA_CB_diangle)
    CB= Atom("CB", carbon_b, 0.0 , 1.0, " "," CB", 0,"C")
    carbon_g= calculateCoordinates(N, CA, CB, CB_CG_length, CA_CB_CG_angle, N_CA_CB_CG_diangle)
    CG= Atom("CG", carbon_g, 0.0, 1.0, " ", " CG", 0, "C")
    carbon_d= calculateCoordinates(CA, CB, CG, CG_CD_length, CB_CG_CD_angle, CA_CB_CG_CD_diangle)
    CD= Atom("CD", carbon_d, 0.0, 1.0, " ", " CD", 0, "C")
    nitrogen_e= calculateCoordinates(CB, CG, CD, CD_NE_length, CG_CD_NE_angle, CB_CG_CD_NE_diangle)
    NE= Atom("NE", nitrogen_e, 0.0, 1.0, " ", " NE", 0, "N")
    carbon_z= calculateCoordinates(CG, CD, NE, NE_CZ_length, CD_NE_CZ_angle, CG_CD_NE_CZ_diangle)
    CZ= Atom("CZ", carbon_z, 0.0, 1.0, " ", " CZ", 0, "C")
    nitrogen_h1= calculateCoordinates(CD, NE, CZ, CZ_NH1_length, NE_CZ_NH1_angle, CD_NE_CZ_NH1_diangle)
    NH1= Atom("NH1", nitrogen_h1, 0.0, 1.0, " ", " NH1", 0, "N")
    nitrogen_h2= calculateCoordinates(CD, NE, CZ, CZ_NH2_length, NE_CZ_NH2_angle, CD_NE_CZ_NH2_diangle)
    NH2= Atom("NH2", nitrogen_h2, 0.0, 1.0, " ", " NH2", 0, "N")

    ##Create Residue Data Structure
    res= Residue((' ', segID, ' '), "ARG", '    ')
    res.add(N)
    res.add(CA)
    res.add(C)
    res.add(O)
    res.add(CB)
    res.add(CG)
    res.add(CD)
    res.add(NE)
    res.add(CZ)
    res.add(NH1)
    res.add(NH2)
    return res
1089341
def align(fastq_file, pair_file, index_dir, names, align_dir, data):
    
    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))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    rg_info = novoalign.get_rg_info(names)
    preset = "sr"

    pair_file = pair_file if pair_file else ""
    if data.get("align_split"):
        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 not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
            index_file = None
            # Skip trying to use indices now as they provide only slight speed-ups
            # and give inconsitent outputs in BAM headers
            # If a single index present, index_dir points to that
            # if index_dir and os.path.isfile(index_dir):
            #     index_dir = os.path.dirname(index_dir)
            #     index_file = os.path.join(index_dir, "%s-%s.mmi" % (dd.get_genome_build(data), preset))
            if not index_file or not os.path.exists(index_file):
                index_file = dd.get_ref_file(data)
            cmd = ("minimap2 -a -x {preset} -R '{rg_info}' -t {num_cores} {index_file} "
                   "{fastq_file} {pair_file} | ")
            do.run(cmd.format(**locals()) + tobam_cl, "minimap2 alignment: %s" % dd.get_sample_name(data))
    data["work_bam"] = out_file
    return data
62596
def alignprotein(self, sample, target):
        
        # Initialise dictionaries
        sample[self.analysistype].dnaseq = dict()
        sample[self.analysistype].protseq = dict()
        sample[self.analysistype].ntindex = dict()
        sample[self.analysistype].aaindex = dict()
        sample[self.analysistype].ntalign = dict()
        sample[self.analysistype].aaalign = dict()
        sample[self.analysistype].aaidentity = dict()
        # Remove any gaps incorporated into the sequence
        sample[self.analysistype].targetsequence[target] = \
            sample[self.analysistype].targetsequence[target].replace('-', '')
        # In order to properly translate the nucleotide sequence, BioPython requests that the sequence is a multiple of
        # three - not partial codons. Trim the sequence accordingly
        remainder = 0 - len(sample[self.analysistype].targetsequence[target]) % 3
        seq = sample[self.analysistype].targetsequence[target] if remainder == 0 \
            else sample[self.analysistype].targetsequence[target][:remainder]
        # Set the DNA and protein sequences of the target in the sample
        sample[self.analysistype].dnaseq[target] = Seq(seq, IUPAC.unambiguous_dna)
        # Translate the nucleotide sequence
        sample[self.analysistype].protseq[target] = str(sample[self.analysistype].dnaseq[target].translate())
        for targetfile in self.targetfiles:
            # Trim the reference sequence to multiples of three
            refremainder = 0 - len(self.records[targetfile][target].seq) % 3
            refseq = str(self.records[targetfile][target].seq) if refremainder % 3 == 0 \
                else str(self.records[targetfile][target].seq)[:refremainder]
            # Translate the nucleotide sequence of the reference sequence
            refdna = Seq(refseq, IUPAC.unambiguous_dna)
            refprot = str(refdna.translate())
            # Use pairwise2 to perform a local alignment with the following parameters:
            # x     No match parameters. Identical characters have score of 1, otherwise 0.
            # s     Same open (-1)  and extend (-.1) gap penalties for both sequences
            ntalignments = pairwise2.align.localxs(seq, refseq, -1, -.1)
            # Use format_alignment to create a formatted alignment that is subsequently split on newlines e.g.
            '''
            ACCGT
            | ||
            A-CG-
            Score=3
            '''
            ntformat = (str(format_alignment(*ntalignments[0])).split('\n'))
            # Align the nucleotide sequence of the reference (ntalignments[2]) to the sample (ntalignments[0]).
            # If the corresponding bases match, add a |, otherwise a space
            ntalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(ntformat[0], ntformat[2])))
            # Create the nucleotide alignment: the sample sequence, the (mis)matches, and the reference sequence
            sample[self.analysistype].ntalign[target] = self.interleaveblastresults(ntformat[0], ntformat[2])
            # Regex to determine location of mismatches in the sequences
            count = 0
            sample[self.analysistype].ntindex[target] = str()
            for snp in re.finditer(' ', ntalignment):
                # If there are many SNPs, then insert line breaks for every 10 SNPs
                if count <= 10:
                    sample[self.analysistype].ntindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].ntindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1
            # Perform the same steps, except for the amino acid sequence
            aaalignments = pairwise2.align.localxs(sample[self.analysistype].protseq[target], refprot, -1, -.1)
            aaformat = (str(format_alignment(*aaalignments[0])).split('\n'))
            aaalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(aaformat[0], aaformat[2])))
            sample[self.analysistype].aaidentity[target] = '{:.2f}'\
                .format(float(aaalignment.count('|')) / float(len(aaalignment)) * 100)
            sample[self.analysistype].aaalign[target] = self.interleaveblastresults(aaformat[0], aaformat[2])
            count = 0
            sample[self.analysistype].aaindex[target] = str()
            for snp in re.finditer(' ', aaalignment):
                if count <= 10:
                    sample[self.analysistype].aaindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].aaindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1
167009
def get_abi_alignment(self, target_data, context=None):
        
        llty = self._get_ll_pointer_type(target_data, context)
        return target_data.get_pointee_abi_alignment(llty)
503687
def align_to_other(self, other, mapping, self_root_pair, other_root_pair = None):
        '''
        
        '''
        if other_root_pair == None:
            other_root_pair = self_root_pair

        assert( len(self_root_pair) == len(other_root_pair) )

        unmoved_atom_names = []
        new_coords = [ None for x in xrange( len(self_root_pair) ) ]
        for atom in self.names:
            if atom in self_root_pair:
                i = self_root_pair.index(atom)
                assert( new_coords[i] == None )
                new_coords[i] = self.get_coords_for_name(atom)

            if atom in mapping:
                other_atom = mapping[atom]
                self.set_coords_for_name( atom, other.get_coords_for_name(other_atom) )
            else:
                unmoved_atom_names.append(atom)

        # Move unmoved coordinates after all other atoms have been moved (so that
        # references will have been moved already)
        if None in new_coords:
            print new_coords
            assert( None not in new_coords )
        ref_coords = [other.get_coords_for_name(x) for x in other_root_pair]

        # Calculate translation and rotation matrices
        U, new_centroid, ref_centroid = calc_rotation_translation_matrices( ref_coords, new_coords )
        for atom in unmoved_atom_names:
            original_coord = self.get_coords_for_name(atom)
            self.set_coords_for_name( atom, rotate_and_translate_coord(original_coord, U, new_centroid, ref_centroid) )
        self.chain = other.chain
585657
def partial_charge(aa, pH):
    
    difference = pH - residue_pka[aa]
    if residue_charge[aa] > 0:
        difference *= -1
    ratio = (10 ** difference) / (1 + 10 ** difference)
    return ratio