Run variant quality score recalibration.
777740
def _run_vqsr(in_file, ref_file, vrn_files, sensitivity_cutoff, filter_type, data):
    
    cutoffs = ["100.0", "99.99", "99.98", "99.97", "99.96", "99.95", "99.94", "99.93", "99.92", "99.91",
               "99.9", "99.8", "99.7", "99.6", "99.5", "99.0", "98.0", "90.0"]
    if sensitivity_cutoff not in cutoffs:
        cutoffs.append(sensitivity_cutoff)
        cutoffs.sort()
    broad_runner = broad.runner_from_config(data["config"])
    gatk_type = broad_runner.gatk_type()
    base = utils.splitext_plus(in_file)[0]
    recal_file = ("%s-vqsrrecal.vcf.gz" % base) if gatk_type == "gatk4" else ("%s.recal" % base)
    tranches_file = "%s.tranches" % base
    plot_file = "%s-plots.R" % base
    if not utils.file_exists(recal_file):
        with file_transaction(data, recal_file, tranches_file, plot_file) as (tx_recal, tx_tranches, tx_plot_file):
            params = ["-T", "VariantRecalibrator",
                      "-R", ref_file,
                      "--mode", filter_type]
            if gatk_type == "gatk4":
                params += ["--variant", in_file, "--output", tx_recal,
                           "--tranches-file", tx_tranches, "--rscript-file", tx_plot_file]
            else:
                params += ["--input", in_file, "--recal_file", tx_recal,
                           "--tranches_file", tx_tranches, "--rscript_file", tx_plot_file]
            params += _get_vqsr_training(filter_type, vrn_files, gatk_type)
            resources = config_utils.get_resources("gatk_variant_recalibrator", data["config"])
            opts = resources.get("options", [])
            if not opts:
                for cutoff in cutoffs:
                    opts += ["-tranche", str(cutoff)]
                for a in _get_vqsr_annotations(filter_type, data):
                    opts += ["-an", a]
            params += opts
            cores = dd.get_cores(data)
            memscale = {"magnitude": 0.9 * cores, "direction": "increase"} if cores > 1 else None
            try:
                broad_runner.new_resources("gatk-vqsr")
                broad_runner.run_gatk(params, log_error=False, memscale=memscale, parallel_gc=True)
            except:  # Can fail to run if not enough values are present to train.
                return None, None
    if gatk_type == "gatk4":
        vcfutils.bgzip_and_index(recal_file, data["config"])
    return recal_file, tranches_file
550630
def score_pairwise(aseq, bseq):
    """"""
    assert len(aseq) == len(bseq)

    # Affine gap penalties -- default values from EMBOSS needle/water
    GAP_OPEN = -10.0
    GAP_EXTEND = -0.5
    GAP_CHARS = frozenset('-.')

    score = 0.0
    in_gap = True   # Don't apply the opening penalty to the N-terminal gap

    for ares, bres in zip(aseq.upper(), bseq.upper()):
        if ares in GAP_CHARS and bres in GAP_CHARS:
            # Both are gaps -- this happens in multiple sequence alignments
            continue
        match = blosum62.get((ares, bres), None)
        if match is None:
            assert GAP_CHARS.intersection((ares, bres)), \
                    "Expected one gap in: " + str((ares, bres))
            # Gap
            if not in_gap:
                score += GAP_OPEN
                in_gap = True
            score += GAP_EXTEND
        else:
            in_gap = False
            score += match

    if in_gap:
        # Correct for a penalty on the C-terminal gap
        score -= GAP_OPEN

    return score
432595
def validate_mutations(vcffile, bams, reffa, chr_reffa, samples, output_format, outfile,
                       mpileup_parameters=mpileup_parser.MPILEUP_DEFAULT_PARAMS):
    """"""
    output_header = "sample chrom pos ref cov A C G T * - + " \
        "val_ref val_alt val_al_type val_al_count val_maf "\
        "most_common_indel most_common_indel_count most_common_indel_maf most_common_indel_type most_common_al " \
        "most_common_al_count most_common_al_maf most_common_count most_common_maf".split()

    # for backwards compatibility
    # if bam or samples is a string, convert to list instead
    if isinstance(samples, six.string_types):
        samples = [samples]
    if isinstance(bams, six.string_types):
        bams = [bams]

    if output_format == 'vcf':
        vcf_reader = vcf.Reader(open(vcffile))
        vcf_reader.samples = samples
        vcf_reader.formats['GT'] = vcf.parser._Format(id='GT', num=1, type='String', desc="Genotype")
        vcf_reader.formats['AD'] = vcf.parser._Format(id='AD', num='R', type='Integer', desc="Allelic depth")
        vcf_reader.formats['DP'] = vcf.parser._Format(id='DP', num=1, type='Integer', desc="Depth")
        vcf_writer = vcf.Writer(outfile, vcf_reader)
    else:
        vcf_reader = open(vcffile)


    if output_format == "sufam":
        outfile.write("\t".join(output_header))
        outfile.write("\n")
    for record in vcf_reader:
        if output_format != 'vcf':
            line = record
            if line.startswith("#CHROM"):
                header = line[1:].rstrip('\n').split("\t")
                # create spoof pyvcf record if vcf_reader is not used
                _Record = namedtuple('Record', header)
            if line.startswith("#"):
                continue
            if len(header) == 0:
                raise(Exception("No header found in vcf file #CHROM not found"))
            # zip all column values, except alt (needs to be list in pyvcf)
            record_args = dict(zip(header, line.rstrip('\n').split("\t")))
            record_args['ALT'] = [record_args['ALT']]
            record = _Record(**record_args)

        # determine type of mutation
        record_type = "snv"
        if len(record.ALT) > 1:
            warnings.warn("Multiple ALT in one record is not implemented - using first")
        if len(record.REF) > len(record.ALT[0]):
            record_type = "deletion"
        elif len(record.ALT[0]) > len(record.REF):
            record_type = "insertion"

        # no coverage results
        no_cov = pd.Series({
            "chrom": str(record.CHROM), "pos": str(record.POS),
            "ref": str(record.REF),
            "cov": 0, "A": 0, "C": 0, "G": 0, "T": 0,
            "val_ref": str(record.REF), "val_alt": str(record.ALT[0]),
            "val_al_type": record_type, "val_al_count": 0, "val_maf": 0})

        # collect mpileup baseparser results per bam
        bps = []
        for i, bam in enumerate(bams):
            sample = samples[i]
            no_cov['sample'] = sample
            bp_lines = mpileup_parser.run_and_parse(bam, str(record.CHROM), str(record.POS), str(record.POS), reffa,
                                                    chr_reffa, mpileup_parameters)
            bpdf = get_baseparser_extended_df(sample, bp_lines, str(record.REF), str(record.ALT[0]))
            if bpdf is None:
                bp = no_cov
            else:
                bp = bpdf.ix[0, :]
            bps += [bp]

        # output call
        if output_format == "vcf":
            _write_bp_vcf(outfile, bps, vcf_writer, record)
        else:
            # only one bam file supported for outputs other than vcf
            _write_bp(outfile, bps[0], output_header, output_format)
259234
def run(vrn_info, calls_by_name, somatic_info, do_plots=True, handle_failures=True):
    
    if "seq2c" in calls_by_name:
        cnv_info = calls_by_name["seq2c"]
    elif "cnvkit" in calls_by_name:
        cnv_info = calls_by_name["cnvkit"]
    else:
        raise ValueError("BubbleTree only currently support CNVkit and Seq2c: %s" % ", ".join(calls_by_name.keys()))
    work_dir = _cur_workdir(somatic_info.tumor_data)
    class OutWriter:
        def __init__(self, out_handle):
            self.writer = csv.writer(out_handle)
        def write_header(self):
            self.writer.writerow(["chrom", "start", "end", "freq"])
        def write_row(self, rec, stats):
            self.writer.writerow([_to_ucsc_style(rec.chrom), rec.start, rec.stop, stats["tumor"]["freq"]])
    vcf_csv = prep_vrn_file(vrn_info["vrn_file"], vrn_info["variantcaller"],
                            work_dir, somatic_info, OutWriter, cnv_info["cns"])
    cnv_csv = _prep_cnv_file(cnv_info["cns"], cnv_info["variantcaller"], work_dir,
                             somatic_info.tumor_data)
    wide_lrr = cnv_info["variantcaller"] == "cnvkit" and somatic_info.normal_bam is None
    return _run_bubbletree(vcf_csv, cnv_csv, somatic_info.tumor_data, wide_lrr, do_plots,
                           handle_failures)
1075615
def compute_quality_score(self, quality):
        """"""
        score = 0
        UNIT = 2
        if 'frequency' in quality:
            # TODO: should be related to frequency.
            if quality['update_in'] < 0:
                score += UNIT
            else:
                score -= UNIT
        if 'tags_count' in quality:
            if quality['tags_count'] > 3:
                score += UNIT
        if 'description_length' in quality:
            if quality['description_length'] > 100:
                score += UNIT
        if 'has_resources' in quality:
            if quality['has_only_closed_or_no_formats']:
                score -= UNIT
            else:
                score += UNIT
            if quality['has_unavailable_resources']:
                score -= UNIT
            else:
                score += UNIT
        if 'discussions' in quality:
            if quality['has_untreated_discussions']:
                score -= UNIT
            else:
                score += UNIT
        if score < 0:
            return 0
        return score
13048
def autoencoder_residual_discrete():
  """"""
  hparams = autoencoder_residual()
  hparams.bottleneck_bits = 1024
  hparams.bottleneck_noise = 0.05
  hparams.add_hparam("discretize_warmup_steps", 16000)
  hparams.add_hparam("bottleneck_kind", "tanh_discrete")
  hparams.add_hparam("isemhash_noise_dev", 0.5)
  hparams.add_hparam("isemhash_mix_prob", 0.5)
  hparams.add_hparam("isemhash_filter_size_multiplier", 2.0)
  hparams.add_hparam("vq_beta", 0.25)
  hparams.add_hparam("vq_decay", 0.999)
  hparams.add_hparam("vq_epsilon", 1e-5)
  return hparams
1270333
function prepareAndRunBenchmarks(benchmarkPatterns, revisions) {
  // Find all benchmark tests to be run.
  let benchmarks = findFiles(LOCAL_DIR('src'), '*/__tests__/*-benchmark.js');
  if (benchmarkPatterns.length !== 0) {
    benchmarks = benchmarks.filter(benchmark =>
      benchmarkPatterns.some(pattern =>
        path.join('src', benchmark).includes(pattern)
      )
    );
  }

  if (benchmarks.length === 0) {
    console.warn(
      'No benchmarks matching: ' +
        `\u001b[1m${benchmarkPatterns.join('\u001b[0m or \u001b[1m')}\u001b[0m`
    );
    return;
  }

  const environments = revisions.map(revision => ({
    revision,
    distPath: prepareRevision(revision),
  }));
  benchmarks.forEach(benchmark => runBenchmark(benchmark, environments));
}
1255864
def quality_comparator(video_data):
        """"""
        def parse_resolution(res: str) -> Tuple[int, ...]:
            return tuple(map(int, res.split('x')))

        raw_resolution = video_data['resolution']
        resolution = parse_resolution(raw_resolution)
        return resolution
828166
def chi2_adaptive_binning(features_0,features_1,number_of_splits_list,systematics_fraction=0.0,title = "title", name="name", PLOT = True, DEBUG = False, transform='StandardScalar'):
	""""""

	max_number_of_splits = np.max(number_of_splits_list)
	#determine how many data points are in each sample
	no_0=features_0.shape[0]
	no_1=features_1.shape[0]

	print("features_0.shape : ", features_0.shape)
	no_dim = features_0.shape[1]
	#Give all samples in file 0 the label 0 and in file 1 the feature 1
	label_0=np.zeros((no_0,1))
	label_1=np.ones((no_1,1))

	#Create an array containing samples and features.
	data_0=np.c_[features_0,label_0]
	data_1=np.c_[features_1,label_1]

	features= np.r_[features_0,features_1]
	labels=    np.r_[label_0,   label_1]

	data=np.r_[data_0,data_1]
	data_same=np.c_[features,labels]

	#print("data : ",data)
	#print("data_same : ", data_same)
	#print("np.sum(data!=data_same) : ",np.sum(data!=data_same))
	assert np.sum(data!=data_same)==0
	assert (no_dim == data.shape[1]-1)

	if no_dim==2:
		plt.scatter(features[:,0],features[:,1], 0.1)
		plt.savefig('test.png')
		plt.clf()

	if transform=='StandardScalar':
		features = preprocessing.scale(features)
		data = np.c_[features,labels]

        if transform=='uniform':
		#data_new2 = data[:,0]
                data_new  = norm_highD_searchsorted(data[:,0])
                for D in range(1,no_dim):
                        temp = norm_highD_searchsorted(data[:,D])
                        data_new = np.c_[data_new,temp]
			#data_new2= np.c_[data_new2,data[:,D]]
		data_new = np.c_[data_new, np.r_[label_0,label_1]]
		#data_new2= np.c_[data_new2,np.r_[label_0,label_1]]

		print("data : ", data)
		data = data_new
		print("data new : ", data)
		#print("data_new2 : ", data_new2)
		#print("np.sum(data!=data_new2) : ",np.sum(data!=data_new2))

	
	np.random.shuffle(data)

	assert (no_dim == data.shape[1]-1)
	labels=data[:,-1]

	X_values= data[:,:-1]
	X_max   = np.amax(data,axis=0)[:-1] 
	X_min   = np.amin(data,axis=0)[:-1]
	X_total_width = (np.subtract(X_max,X_min))
	del data

	if transform=='fill01':
		#Scaling	
		X_values = X_values - X_min[None,:]
		X_values = X_values / X_total_width[None,:]
	
	if True:
		X_min = [0.]*no_dim
		X_total_width = [1.]*no_dim
	#b = X_values[:,0]
	#print("b[b[:]>2].shape[0] : \n", b[b[:]>2].shape[0] )	
	data = np.concatenate((X_values, labels[:,None]), axis=1)

	if no_dim==2:
		plt.scatter(data[:,0],data[:,1],0.1)
		plt.savefig('test_scaled.png')

	#print("X_values.shape : ",X_values.shape)
	starting_boundary = []
	for i in range(no_dim):
		starting_boundary.append([0.0,1.0])
	#Each key has the following stricture: # of splits and for each split if it was closer (a) or further away from (b) the origin. The original bin is "0"
	#For example "2ab" means it is the bin that was closer to the origin for the first split and further away for the second one.
	bin_boundaries_dict = {'0' : np.array(starting_boundary)}
	bin_points_dict = {'0' : data}

	for split_number in range(1,1+max_number_of_splits):
		for bin_key, bin_boundary in bin_boundaries_dict.items():
			if str(split_number-1) in bin_key:
				variances= np.var(bin_points_dict[bin_key][:,:-1], axis=0)
				#print("\nvariances : ", variances)
				dim_to_be_sliced = np.argmax(variances)
				#print("dim_to_be_sliced : ",dim_to_be_sliced)
				#print("bin_points_dict[bin_key] : ",bin_points_dict[bin_key])
				#print("bin_points_dict[bin_key][:,dim_to_be_sliced] : ",bin_points_dict[bin_key][:,dim_to_be_sliced])
				median = np.median(bin_points_dict[bin_key][:,dim_to_be_sliced])
				#print("median : ",median)

				a_bin_boundary, b_bin_boundary = bin_boundary.copy(), bin_boundary.copy()
				#print("a_bin_boundary : ",a_bin_boundary)
				a_bin_boundary[dim_to_be_sliced,1] = median
				b_bin_boundary[dim_to_be_sliced,0] = median
				bin_boundaries_dict[str(split_number)+bin_key[1:]+'a'] = a_bin_boundary
				bin_boundaries_dict[str(split_number)+bin_key[1:]+'b'] = b_bin_boundary
				
				a_points, b_points = [],[] 
				for event_number in range(bin_points_dict[bin_key].shape[0]):
					if bin_points_dict[bin_key][event_number,dim_to_be_sliced] < median: a_points.append(bin_points_dict[bin_key][event_number,:].tolist())
					else: b_points.append(bin_points_dict[bin_key][event_number,:].tolist())
				
				bin_points_dict[str(split_number)+bin_key[1:]+'a'] = np.array(a_points)
				bin_points_dict[str(split_number)+bin_key[1:]+'b'] = np.array(b_points)
				#If a bin contains no particles it should be deleted
				if len(a_points)==0: 
					del bin_points_dict[str(split_number)+bin_key[1:]+'a']
					del bin_boundaries_dict[str(split_number)+bin_key[1:]+'a']
				if len(b_points)==0: 
					del bin_points_dict[str(split_number)+bin_key[1:]+'b']
					del bin_boundaries_dict[str(split_number)+bin_key[1:]+'b']

	if PLOT:	pickle.dump( bin_boundaries_dict, open( "bin_boundaries_dict.p", "wb" ) )	
	bins_sample01_dict= {}
	signed_Scp2_dict= {}
	results_list = []

	for number_of_splits in number_of_splits_list:
		print("\nnumber_of_splits : ",number_of_splits,"\nsystematics_fraction : ",systematics_fraction)
		bins_sample0, bins_sample1 = [] , [] 
		for bin_key, bin_points in bin_points_dict.items():
			if str(number_of_splits) in bin_key:
				labels_in_bin = bin_points[:,-1]
				#print("labels_in_bin : ",labels_in_bin)
				bin_sample0 = np.count_nonzero( labels_in_bin == 0)
				bin_sample1 = np.count_nonzero( labels_in_bin == 1)
				#print("bin_sample0 : ",bin_sample0) 
				#print("bin_sample1 : ",bin_sample1) 
				#simulate uncertainties 
				if(systematics_fraction*float(bin_sample0)!=0.): bin_sample0 += int(round(np.random.normal(0.,systematics_fraction*float(bin_sample0))))
				if(systematics_fraction*float(bin_sample1)!=0.): bin_sample1 += int(round(np.random.normal(0.,systematics_fraction*float(bin_sample1))))
				bins_sample01_dict[bin_key]=[bin_sample0,bin_sample1]
				signed_Scp2_dict[bin_key] = np.square(float(bin_sample1-bin_sample0))/(float(bin_sample1)+float(bin_sample0)+np.square(float(bin_sample1)*systematics_fraction)+np.square(float(bin_sample1)*systematics_fraction))*np.sign(bin_sample1-bin_sample0)
				#print("\n\nbin_sample0 : ",bin_sample0, "\n bins_sample0 : ", bins_sample0 )
				#print("type(bin_sample0) : ",type(bin_sample0), "  type(bins_sample0) : ",type(bins_sample0))
				bins_sample0.append(bin_sample0)
				#print("  bins_sample0 : ", bins_sample0, "\n\n" )
				bins_sample1.append(bin_sample1)
		bins_sample0, bins_sample1 = np.array(bins_sample0,dtype=float), np.array(bins_sample1, dtype=float)
		print("bins_sample0 : ",bins_sample0,"\n bins_sample1 : ",bins_sample1)
		#element wise subtraction and division
		Scp2 = ((bins_sample1-bins_sample0)**2)/ (bins_sample1+bins_sample0+(systematics_fraction*bins_sample1)**2+(systematics_fraction*bins_sample0)**2 )
		#Scp2 =  np.divide(np.square(np.subtract(bins_sample1,bins_sample0)),np.add(bins_sample1,bins_sample0))
		if DEBUG:
			print(Scp2)

		#nansum ignores all the contributions that are Not A Number (NAN)
		Chi2 = np.nansum(Scp2)
		if DEBUG:
			print("Chi2")
			print(Chi2)
		dof=bins_sample0.shape[0]-1
	
		pvalue= 1 - stats.chi2.cdf(Chi2,dof)
		
		print("\nThe p value for Scp2 = ",Scp2," and Chi2 = ", Chi2, " is ",pvalue,"\n\n")
		if DEBUG:
			print(bins_sample0)
			print(bins_sample1)
			print("Chi2/dof : {0}".format(str(Chi2/dof)))

			print("pvalue : {0}".format(str(pvalue)))

		results_list.append(pvalue)

		if PLOT: 
			if no_dim==1: chi2_plots.adaptive_binning_1Dplot(bin_boundaries_dict,data,number_of_splits,title+" "+str(no_dim) + "D "+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits")
			if no_dim==2: chi2_plots.adaptive_binning_2Dplot(bin_boundaries_dict,signed_Scp2_dict,number_of_splits,X_values,title+" "+str(no_dim) + "D"+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits", X_min= X_min,X_total_width=X_total_width )
                        if no_dim>1:  chi2_plots.adaptive_binning_2D1Dplot(bin_boundaries_dict,bins_sample01_dict,number_of_splits,X_values,title+" "+str(no_dim) + "D"+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits", no_dim)



	return results_list
281637
def iter_variants(self):
        """"""
        if not self.has_index:
            raise NotImplementedError("Not implemented when IMPUTE2 file is "
                                      "not indexed (see genipe)")

        for name, row in self._impute2_index.iterrows():
            # Seeking to the right place in the file
            f = self._impute2_file
            f.seek(int(row.seek))
            chrom, name, pos, a1, a2 = f.read(1024).split(" ")[:5]
            pos = int(pos)

            yield Variant(name, CHROM_STR_ENCODE.get(chrom, chrom), pos,
                          [a1, a2])
1326005
def _select_better_fit(self, matching_venvs):
        """"""
        # keep the venvs in a separate array, to pick up the winner, and the (sorted, to compare
        # each dependency with its equivalent) in other structure to later compare
        venvs = []
        to_compare = []
        for matching, venv in matching_venvs:
            to_compare.append(sorted(matching, key=lambda req: getattr(req, 'key', '')))
            venvs.append(venv)

        # compare each n-tuple of dependencies to see which one is bigger, and add score to the
        # position of the winner
        scores = [0] * len(venvs)
        for dependencies in zip(*to_compare):
            if not isinstance(dependencies[0], Distribution):
                # only distribution URLs can be compared
                continue

            winner = dependencies.index(max(dependencies))
            scores[winner] = scores[winner] + 1

        # get the rightmost winner (in case of ties, to select the latest venv)
        winner_pos = None
        winner_score = -1
        for i, score in enumerate(scores):
            if score >= winner_score:
                winner_score = score
                winner_pos = i
        return venvs[winner_pos]
103912
def _get_variant_effect(cls, variants, ref_sequence):
        ''''''
        assert len(variants) != 0

        var_types = [x.var_type for x in variants]
        if len(set(var_types)) != 1:
            return 'MULTIPLE', '.', '.'

        var_type = var_types[0]

        assert set([x.ref_name for x in variants]) == set([ref_sequence.id])
        codon_starts = [AssemblyVariants._get_codon_start(0, x.ref_start) for x in variants]
        assert len(set(codon_starts)) == 1
        codon_start = codon_starts[0]
        aa_start = codon_start // 3
        ref_codon = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start:codon_start+3])
        ref_aa = ref_codon.translate()

        if var_type == pymummer.variant.SNP:
            new_codon = list(ref_codon.seq)
            for v in variants:
                new_codon[v.ref_start - codon_start] = v.qry_base
            new_codon = pyfastaq.sequences.Fasta('new', ''.join(new_codon))
            qry_aa = new_codon.translate()

            if ref_aa.seq == qry_aa.seq:
                return ('SYN', '.', aa_start)
            elif qry_aa.seq == '*':
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            else:
                return ('NONSYN', ref_aa.seq + str(aa_start + 1) + qry_aa.seq, aa_start)
        elif var_type in [pymummer.variant.INS, pymummer.variant.DEL]:
            if len(variants) > 1:
                return 'INDELS', '.', aa_start

            var = variants[0]

            if var_type == pymummer.variant.INS:
                new_seq = pyfastaq.sequences.Fasta('seq', var.qry_base)
            else:
                new_seq = pyfastaq.sequences.Fasta('seq', var.ref_base)

            if len(new_seq) % 3 != 0:
                return ('FSHIFT', ref_aa.seq + str(aa_start + 1) + 'fs', aa_start)

            new_seq_aa = new_seq.translate()
            if '*' in new_seq_aa.seq:
                return ('TRUNC', ref_aa.seq + str(aa_start + 1) + 'trunc', aa_start)
            elif var_type == pymummer.variant.INS:
                ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                aa_after_ins = ref_codon_after_ins.translate()
                return ('INS', ref_aa.seq + str(aa_start + 1) + '_' + aa_after_ins.seq + str(aa_start + 2) + 'ins' + new_seq_aa.seq , aa_start)
            else:
                if len(new_seq) == 3:
                    return ('DEL', ref_aa.seq + str(aa_start + 1) + 'del', aa_start)
                else:
                    assert len(new_seq) % 3 == 0
                    ref_codon_after_ins = pyfastaq.sequences.Fasta('codon', ref_sequence[codon_start+3:codon_start+6])
                    aa_after_ins = ref_codon_after_ins.translate()
                    return ('DEL', ref_aa.seq + str(aa_start + 1)+ '_' + aa_after_ins.seq + str(aa_start + 2) + 'del', aa_start)

        else:
            return ('UNKNOWN', '.', aa_start)
868830
def _run_variantcall_batch_multicore(items, regions, final_file):
    
    batch_name = _get_batch_name(items)
    variantcaller = _get_batch_variantcaller(items)
    work_bams = [dd.get_work_bam(d) or dd.get_align_bam(d) for d in items]
    def split_fn(data):
        out = []
        for region in regions:
            region = _region_to_coords(region)
            chrom, start, end = region
            region_str = "_".join(str(x) for x in region)
            out_file = os.path.join(dd.get_work_dir(items[0]), variantcaller, chrom,
                                    "%s-%s.vcf.gz" % (batch_name, region_str))
            out.append((region, work_bams, out_file))
        return final_file, out
    parallel = {"type": "local", "num_jobs": dd.get_num_cores(items[0]), "cores_per_job": 1}
    run_parallel = dmulti.runner(parallel, items[0]["config"])
    to_run = copy.deepcopy(items[0])
    to_run["sam_ref"] = dd.get_ref_file(to_run)
    to_run["group_orig"] = items
    parallel_split_combine([[to_run]], split_fn, run_parallel,
                           "variantcall_sample", "concat_variant_files",
                           "vrn_file", ["region", "sam_ref", "config"])
    return final_file
132642
def iter_variants(self):
        """"""
        for v in self.get_vcf():
            yield Variant(v.ID, v.CHROM, v.POS, {v.REF} | set(v.ALT))
190783
def _convert_to_sequence_variants(self, variant, acc):
        
        start = variant['start']
        insertion = ''.join(variant['ins'])
        deletion = ''.join(variant['del'])

        # defaults
        is_dup = False    # assume not dup
        fsext_len = None    # fs or ext length
        is_sub = False
        is_ext = False

        if start == 1:    # initial methionine is modified
            aa_start = aa_end = AAPosition(base=start, aa=deletion)
            ref = ''
            alt = ''
            self._is_init_met = True
            self._is_ambiguous = True    # side-effect

        if insertion and insertion.find("*") == 0:    # stop codon at variant position
            aa_start = aa_end = AAPosition(base=start, aa=deletion[0])
            ref = ''
            alt = '*'
            is_sub = True

        elif start == len(self._ref_seq):    # extension
            if self._alt_seq[-1] == '*':
                fsext_len = len(insertion) - len(deletion)    # don't include the former stop codon
            else:
                fsext_len = '?'
            subst_at_stop_codon = insertion[0]

            aa_start = aa_end = AAPosition(base=start, aa='*')
            ref = ''
            alt = subst_at_stop_codon
            is_ext = True

        elif self._is_frameshift:    # frameshift
            aa_start = aa_end = AAPosition(base=start, aa=deletion[0])
            ref = ''

            try:
                fsext_len = str(
                    insertion.index("*") + 1)    # start w/ 1st change; ends w/ * (inclusive)
            except ValueError:
                fsext_len = "?"

            alt = insertion[0]

        else:    # no frameshift - sub/delins/dup
            if insertion == deletion:    # silent
                aa_start = aa_end = AAPosition(base=start, aa=deletion)
                ref = alt = ''

            elif len(insertion) == len(deletion) == 1:    # substitution
                aa_start = aa_end = AAPosition(base=start, aa=deletion)
                ref = ''
                alt = insertion
                is_sub = True

            elif len(deletion) > 0:    # delins OR deletion OR stop codon at variant position
                ref = deletion
                end = start + len(deletion) - 1
                if len(insertion) > 0:    # delins
                    aa_start = AAPosition(base=start, aa=deletion[0])
                    if end > start:
                        aa_end = AAPosition(base=end, aa=deletion[-1])
                    else:
                        aa_end = aa_start
                    alt = insertion

                else:    # deletion OR stop codon at variant position
                    if len(deletion) + start == len(
                            self._ref_seq):    # stop codon at variant position
                        aa_start = AAPosition(base=start, aa=deletion[0])
                        aa_end = AAPosition(base=start, aa=deletion[0])
                        ref = ''
                        alt = '*'
                        is_sub = True
                    else:    # deletion
                        aa_start = AAPosition(base=start, aa=deletion[0])
                        if end > start:
                            aa_end = AAPosition(base=end, aa=deletion[-1])
                        else:
                            aa_end = aa_start
                        alt = None

            elif len(deletion) == 0:    # insertion OR duplication OR extension

                is_dup, dup_start = self._check_if_ins_is_dup(start, insertion)

                if is_dup:    # duplication
                    dup_end = dup_start + len(insertion) - 1
                    aa_start = AAPosition(base=dup_start, aa=insertion[0])
                    aa_end = AAPosition(base=dup_end, aa=insertion[-1])
                    ref = alt = None

                else:    # insertion
                    start -= 1
                    end = start + 1

                    aa_start = AAPosition(base=start, aa=self._ref_seq[start - 1])
                    aa_end = AAPosition(base=end, aa=self._ref_seq[end - 1])
                    ref = None
                    alt = insertion

            else:    # should never get here
                raise ValueError("unexpected variant: {}".format(variant))

        var_p = self._create_variant(
            aa_start,
            aa_end,
            ref,
            alt,
            fsext_len=fsext_len,
            is_dup=is_dup,
            acc=acc,
            is_ambiguous=self._is_ambiguous,
            is_sub=is_sub,
            is_ext=is_ext,
            is_init_met=self._is_init_met
        )

        return var_p
448871
def runExperiment():
  
  trainSeqN = [5, 10, 20, 50, 100, 200]
  rptPerCondition = 20
  correctRateAll = np.zeros((len(trainSeqN), rptPerCondition))
  missRateAll = np.zeros((len(trainSeqN), rptPerCondition))
  fpRateAll = np.zeros((len(trainSeqN), rptPerCondition))
  for i in xrange(len(trainSeqN)):
    for rpt in xrange(rptPerCondition):
      numTrainSequence = trainSeqN[i]
      correctRate, missRate, fpRate = runSingleExperiment(numTrainSequence=numTrainSequence)
      correctRateAll[i, rpt] = correctRate
      missRateAll[i, rpt] = missRate
      fpRateAll[i, rpt] = fpRate

  plt.figure()
  plt.subplot(2,2,1)
  plt.semilogx(trainSeqN, 100*np.mean(correctRateAll,1),'-*')
  plt.xlabel(' Training Sequence Number')
  plt.ylabel(' Hit Rate - Best Match (%)')
  plt.subplot(2,2,2)
  plt.semilogx(trainSeqN, 100*np.mean(missRateAll,1),'-*')
  plt.xlabel(' Training Sequence Number')
  plt.ylabel(' Miss Rate (%)')
  plt.subplot(2,2,3)
  plt.semilogx(trainSeqN, 100*np.mean(fpRateAll,1),'-*')
  plt.xlabel(' Training Sequence Number')
  plt.ylabel(' False Positive Rate (%)')
  plt.savefig('result/ReberSequence_HMMperformance.pdf')
  plt.show()
1196008
def evaluate(data_loader):
    
    translation_out = []
    all_inst_ids = []
    avg_loss_denom = 0
    avg_loss = 0.0
    for _, (src_seq, tgt_seq, src_valid_length, tgt_valid_length, inst_ids) \
            in enumerate(data_loader):
        src_seq = src_seq.as_in_context(ctx)
        tgt_seq = tgt_seq.as_in_context(ctx)
        src_valid_length = src_valid_length.as_in_context(ctx)
        tgt_valid_length = tgt_valid_length.as_in_context(ctx)
        # Calculating Loss
        out, _ = model(src_seq, tgt_seq[:, :-1], src_valid_length, tgt_valid_length - 1)
        loss = loss_function(out, tgt_seq[:, 1:], tgt_valid_length - 1).mean().asscalar()
        all_inst_ids.extend(inst_ids.asnumpy().astype(np.int32).tolist())
        avg_loss += loss * (tgt_seq.shape[1] - 1)
        avg_loss_denom += (tgt_seq.shape[1] - 1)
        # Translate
        samples, _, sample_valid_length =\
            translator.translate(src_seq=src_seq, src_valid_length=src_valid_length)
        max_score_sample = samples[:, 0, :].asnumpy()
        sample_valid_length = sample_valid_length[:, 0].asnumpy()
        for i in range(max_score_sample.shape[0]):
            translation_out.append(
                [tgt_vocab.idx_to_token[ele] for ele in
                 max_score_sample[i][1:(sample_valid_length[i] - 1)]])
    avg_loss = avg_loss / avg_loss_denom
    real_translation_out = [None for _ in range(len(all_inst_ids))]
    for ind, sentence in zip(all_inst_ids, translation_out):
        real_translation_out[ind] = sentence
    return avg_loss, real_translation_out