Add phenotype information to a VCF header.

    Encode tumor/normal relationships in VCF header.
    Could also eventually handle more complicated pedigree information if useful.
349050
def _add_vcf_header_sample_cl(in_file, items, base_file):
    
    paired = vcfutils.get_paired(items)
    if paired:
        toadd = ["##SAMPLE=<ID=%s,Genomes=Tumor>" % paired.tumor_name]
        if paired.normal_name:
            toadd.append("##SAMPLE=<ID=%s,Genomes=Germline>" % paired.normal_name)
            toadd.append("##PEDIGREE=<Derived=%s,Original=%s>" % (paired.tumor_name, paired.normal_name))
        new_header = _update_header(in_file, base_file, toadd, _fix_generic_tn_names(paired))
        if vcfutils.vcf_has_variants(in_file):
            cmd = "bcftools reheader -h {new_header} | bcftools view "
            return cmd.format(**locals())
147705
def prep_vrn_file(in_file, vcaller, work_dir, somatic_info, writer_class, seg_file=None, params=None):
    
    data = somatic_info.tumor_data
    if not params:
        params = PARAMS
    out_file = os.path.join(work_dir, "%s-%s-prep.csv" % (utils.splitext_plus(os.path.basename(in_file))[0],
                                                          vcaller))
    if not utils.file_uptodate(out_file, in_file):
        # ready_bed = _identify_heterogeneity_blocks_seg(in_file, seg_file, params, work_dir, somatic_info)
        ready_bed = None
        if ready_bed and utils.file_exists(ready_bed):
            sub_file = _create_subset_file(in_file, ready_bed, work_dir, data)
        else:
            sub_file = in_file
        max_depth = max_normal_germline_depth(sub_file, params, somatic_info)
        with file_transaction(data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                writer = writer_class(out_handle)
                writer.write_header()
                bcf_in = pysam.VariantFile(sub_file)
                for rec in bcf_in:
                    stats = _is_possible_loh(rec, bcf_in, params, somatic_info, max_normal_depth=max_depth)
                    if chromhacks.is_autosomal(rec.chrom) and stats is not None:
                        writer.write_row(rec, stats)
    return out_file
253332
def genotypesPhenotypesGenerator(self, request):
        
        # TODO make paging work using SPARQL?
        compoundId = datamodel.PhenotypeAssociationSetCompoundId.parse(
            request.phenotype_association_set_id)
        dataset = self.getDataRepository().getDataset(compoundId.dataset_id)
        phenotypeAssociationSet = dataset.getPhenotypeAssociationSet(
            compoundId.phenotypeAssociationSetId)
        featureSets = dataset.getFeatureSets()
        annotationList = phenotypeAssociationSet.getAssociations(
            request, featureSets)
        return self._protocolListGenerator(request, annotationList)
160592
def variants_to_vcf(variants, variant_to_metadata, out=sys.stdout):
    

    # TODO: The variant metadata dictionary (in the `VariantCollection`)
    # contains different data depending on the original input file format (VCF,
    # MAF, CSV).  It's easy to output variants originally in VCF format, but we
    # might want to consider how fields from MAF (for example) map to those in
    # VCF.

    # TODO: The VCF file we output doesn't contain any VCF metadata headers, as
    # the original headers were thrown away when the VCF file was parsed. We
    # may want to keep some of that information and/or infer some of the
    # headers based on the variants themselves. The former is difficult because
    # merge conflicts will inevitably occur; the latter is difficult because
    # the variants themselves don't contain all the information required for
    # these metadata headers (e.g., descriptions).
    #
    # As a side note, adding headers for certain fields will make them parse
    # correctly -- into an integer instead of a string (the default), for
    # example.

    # TODO: If we maintain headers (see above TODO), what should happen if
    # variant sources use different reference genomes?
    #
    # If we don't maintain headers, what should the default reference genome
    # be? This code chose one fairly arbitrarily.

    # TODO: If we end up needing more functions to "build" VCF record fields
    # (see functions below), we may want to abstract away the individual
    # functions and create a mapping from field to format function.

    def get_metadata_field(key, variant, default='.'):
        """Retrieve field from variant metadata dictionary."""
        val = variant_to_metadata[variant].get(key)
        if val is None:
            return default
        return val

    def build_filter_field(variant):
        """Build the filter field from the given variant.

        The `filter` field, as stored in the variants metadata dictionary,
        comes in 3 flavors:

            - empty list: 1+ filters were run and none failed
            - non-empty list: 1+ filters were run and 1+ failed
            - `None`: no filters were run

        This function maps each of these internal representations to their
        corresponding VCF representations.
        """
        filter_metadata = get_metadata_field('filter', variant)
        if type(filter_metadata) == list:
            return 'PASS' if filter_metadata == [] else ';'.join(filter_metadata)
        else:
            # TODO: Can the filter field ever be something other than the 3
            # possibilities described in this function's doc comment?
            bad_value_msg = (
                'Filter metadata field took on unexpected value `{}`. Update '
                'code to account for this value.').format(str(filter_metadata))
            assert filter_metadata == '.', bad_value_msg
            return filter_metadata

    def build_info_field(variant):
        """Build the info field from the given variant.

        Format is `<key>=<val>,...;<key>=<val>,...;<key>=<val>,...`.
        """

        def build_info_pair(key, val):
            """Build key/val pair for info field."""
            # Note: Different from `val == True`, which returns True when `val == 1`.
            if val is True:
                return key

            if type(val) == list:
                val = ','.join(map(str, val))
            else:
                val = str(val)
            return '{}={}'.format(key, val)

        info_dict = get_metadata_field('info', variant, default={})
        if not info_dict:
            return '.'

        return ';'.join(build_info_pair(k, v) for (k, v) in info_dict.items())

    def build_format_field(variant):
        """Build the sample format string from the given variant.

        Each sample column follows this format for the specified variant.
        """
        sample_info = get_metadata_field('sample_info', variant, default={})
        return ':'.join(list(sample_info.values())[0].keys()) if sample_info else '.'

    def build_sample_fields(variant):
        """Build the sample fields for the given variant."""

        def build_sample_field(sample):
            """Build a specific sample's field (i.e., one sample column)."""
            sample_vals = sample_info[sample].values()
            return ':'.join(build_sample_val(val) for val in sample_vals)

        def build_sample_val(sample_val):
            """Build a specific value for a sample (i.e., one value for one column)."""
            if type(sample_val) is list:
                return ','.join(map(str, sample_val))
            elif sample_val is not None:
                return str(sample_val)
            else:
                return '.'

        sample_info = get_metadata_field('sample_info', variant, default={})
        return list(build_sample_field(sample) for sample in sample_info)

    def build_vcf_record(variant, add_sample_info):
        """Return a list of all the variant's VCF record fields."""
        record = [
            str(variant.original_contig),
            str(variant.original_start),
            get_metadata_field('id', variant),
            variant.original_ref,
            variant.original_alt,
            str(get_metadata_field('qual', variant)),
            build_filter_field(variant),
            build_info_field(variant),
        ]
        if add_sample_info:
            record.append(build_format_field(variant))
            record.extend(build_sample_fields(variant))
        return record

    def merge_duplicate_variants():
        """Merge duplicate variants (according to ID) and return *list* of merged, sorted variants.

        Multiple `Variant`s can have the same VCF id (e.g., those variants which
        originally had > 1 alternate base), but we can't output a VCF file with multiple
        records having the same id. This function merges those duplicate variants.
        """

        def construct_id2variants():
            """Construct dictionary which maps variant IDs to `Variant`s."""
            id2variants = defaultdict(list)
            for variant in variants:
                # Note: For variants without IDs (e.g., those that came from
                # MAF files), we assign the missing value.
                variant_id = get_metadata_field('id', variant)
                id2variants[variant_id].append(variant)
            return id2variants

        def merge_variant_list(duplicates):
            """Merge duplicate variant list into one."""

            # TODO: Currently assumes that only alternate bases differ, but we may
            # want or need to merge variants that have the same ID but different
            # contigs, positions, or reference bases. If merging isn't possible (which
            # I think is the case for many situations), it may be easiest to assign
            # the "missing value" to the IDs and move on.

            assert len(duplicates) > 0
            assert all(v.original_contig == duplicates[0].original_contig and
                       v.original_start == duplicates[0].original_start and
                       v.original_ref == duplicates[0].original_ref for v in duplicates)
            import copy
            merged = copy.copy(duplicates[0])
            merged.original_alt = ','.join(duplicate.original_alt for duplicate in duplicates)
            return merged

        id2variants = construct_id2variants()
        variants_no_id = id2variants.pop('.', [])  # don't want to merge variants w/ no id
        merged_variants = list(map(merge_variant_list, id2variants.values())) + variants_no_id

        # groups variants by contig; relative ordering of contigs doesn't matter
        return sorted(
            merged_variants,
            key=lambda v: (str(v.original_contig), str(v.original_start)))

    def get_sample_names():
        """Return the sample names for all variants."""

        # TODO: For now, ensures that every variant has the same samples. If they didn't,
        # we'd have to use the missing value for VCF records with no data on a given
        # sample.  In and of itself, that isn't a problem until the format field contains
        # GT (genotype), which is required for all samples if present in the format field.
        #
        # A couple of ways to handle this:
        #
        #   1) Ignore this requirement of VCF files. I'd assume this would still be
        #   compatible with varcode, but may not be with other VCF parsers.
        #   2) Remove the GT field from those records where this rule would be violated.
        #   This is sad because we lose information.
        #
        # It's also important to note that if/when we combine all samples, we'll have to
        # add a `.` (missing value) for those samples in which a variant has no data.
        # Check the VCF spec to make sure this is valid; if not, we may have to write
        # `.:.:.:.` -- one `.` for each field in the format string.
        #
        # Moreover, we'll want to note whether we care about maintaining the relative
        # ordering of the sample names in the original VCF files. This probably isn't
        # important, but the code below does not do this (it effectively alphabetizes the
        # sample names because of the use of `set`).

        sample_names = set()
        for variant in variants:
            sample_info = get_metadata_field('sample_info', variant, default={})
            addl_sample_names = set(sample_info.keys())

            # Ensures all variants have the same samples.
            if sample_names and sample_names != addl_sample_names:
                return []

            sample_names.update(addl_sample_names)

        return list(sample_names)

    headers = ['#CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']
    sample_names = get_sample_names()
    if sample_names:
        headers += ['FORMAT'] + sample_names

    unique_variants_list = merge_duplicate_variants()

    # usually we have just one reference genome for the whole variant collection
    # but if variants from multiple sources have been merged then we might
    # not be able to write out a VCF since the individual variants may be using
    # different coordinate systems
    genome_names = list({v.ensembl.reference_name for v in unique_variants_list})
    if len(genome_names) > 1:
        raise ValueError(
            "Cannot create VCF for variants with multiple reference genomes: %s" % (
                genome_names,))
    genome_name = genome_names[0]
    print('##fileformat=VCFv4.2', file=out)
    print('##reference=%s' % genome_name, file=out)
    print('\t'.join(headers), file=out)
    for variant in unique_variants_list:
        record = build_vcf_record(variant, add_sample_info=bool(sample_names))
        print('\t'.join(record), file=out)
325134
def run(_, data, out_dir):
    
    out = []
    vcinfo = get_active_vcinfo(data)
    if vcinfo:
        if dd.get_phenotype(data) == "normal" and "germline" in vcinfo:
            out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        elif dd.get_phenotype(data) != "germline":
            out.append(_bcftools_stats(data, out_dir))
            if "germline" in vcinfo:
                out.append(_bcftools_stats(data, out_dir, "germline", germline=True))
        else:
            out.append(_bcftools_stats(data, out_dir, germline=True))
    out.append(_snpeff_stats(data, out_dir))

    out = [item for item in out if item]
    if out:
        return {"base": out[0], "secondary": out[1:]}
1087769
def _vcf_info(start, end, mate_id, info=None):
    
    out = "SVTYPE=BND;MATEID={mate};IMPRECISE;CIPOS=0,{size}".format(
        mate=mate_id, size=end-start)
    if info is not None:
        extra_info = ";".join("{0}={1}".format(k, v) for k, v in info.iteritems())
        out = "{0};{1}".format(out, extra_info)
    return out
155014
static public void addGmodGFF3GeneFeatures(LinkedHashMap<String, ChromosomeSequence> chromosomeSequenceList, FeatureList listGenes) throws Exception {


		// key off mRNA as being a known feature that may or may not have a parent gene


		FeatureList mRNAFeatures = listGenes.selectByType("mRNA");
		LinkedHashMap<String,FeatureList> featureIDHashMap = FeatureHelper.buildFeatureAtrributeIndex("ID", listGenes);
		LinkedHashMap<String,FeatureList> featureParentHashMap = FeatureHelper.buildFeatureAtrributeIndex("Parent", listGenes);

		for (FeatureI f : mRNAFeatures) {
			String geneID;
			String geneNote = null;
			String geneSource = null;
			String sequenceName = null;
			ChromosomeSequence seq = null;
			GeneSequence geneSequence = null;

			Feature mRNAFeature = (Feature) f;
			String mRNAID = mRNAFeature.getAttribute("ID");
			String mRNAsource = mRNAFeature.source();
			String mRNANote = mRNAFeature.getAttribute("Note");
			String mRNAParent = mRNAFeature.getAttribute("Parent");
			if (mRNAParent != null && mRNAParent.length() > 0) {
			   // FeatureList geneFeatureList = listGenes.selectByAttribute("ID", mRNAParent);
				FeatureList geneFeatureList = featureIDHashMap.get(mRNAParent);
				Feature geneFeature = (Feature) geneFeatureList.get(0);
				geneID = geneFeature.getAttribute("ID");
				geneNote = geneFeature.getAttribute("Note");
				geneSource = geneFeature.source();
				sequenceName = geneFeature.seqname();

				//
			} else {
				//deal with cases where no parent gene is given
				geneID = mRNAID;
				geneSource = mRNAsource;
				sequenceName = mRNAFeature.seqname();
			}

			seq = chromosomeSequenceList.get(sequenceName);

			AccessionID geneAccessionID = new AccessionID(geneID);

		  //  FeatureList mRNAChildren = listGenes.selectByAttribute("Parent", mRNAID);
			FeatureList mRNAChildren = featureParentHashMap.get(mRNAID);
			FeatureList cdsFeatures = mRNAChildren.selectByType("CDS");
			FeatureI feature = cdsFeatures.get(0);
			Strand strand = Strand.POSITIVE;

			if (feature.location().isNegative()) {
				strand = Strand.NEGATIVE;
			}
			cdsFeatures = cdsFeatures.sortByStart();







			//String seqName = feature.seqname();
			FeatureI startCodon = null;
			FeatureI stopCodon = null;
			Integer startCodonBegin = null;
			Integer stopCodonEnd = null;
			String startCodonName = "";
			String stopCodonName = "";
			FeatureList startCodonList = mRNAChildren.selectByType("five_prime_UTR");
			if (startCodonList != null && startCodonList.size() > 0) {
				startCodon = startCodonList.get(0);
				if (strand == Strand.NEGATIVE) {
					startCodonBegin = startCodon.location().bioEnd();
				} else {
					startCodonBegin = startCodon.location().bioStart();
				}
				startCodonName = startCodon.getAttribute("ID");
			}

			FeatureList stopCodonList = mRNAChildren.selectByType("three_prime_UTR");

			if (stopCodonList != null && stopCodonList.size() > 0) {
				stopCodon = stopCodonList.get(0);
				if (strand == Strand.NEGATIVE) {
					stopCodonEnd = stopCodon.location().bioStart();
				} else {
					stopCodonEnd = stopCodon.location().bioEnd();
				}
				stopCodonName = stopCodon.getAttribute("ID");

			}




			if (startCodonBegin == null) {
				if (strand == Strand.NEGATIVE) {
					FeatureI firstFeature = cdsFeatures.get(0);

					startCodonBegin = firstFeature.location().bioEnd();
				} else {
					FeatureI firstFeature = cdsFeatures.get(0);

					startCodonBegin = firstFeature.location().bioStart();
				}
			}

			if (stopCodonEnd == null) {
				if (strand == Strand.NEGATIVE) {
					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					stopCodonEnd = lastFeature.location().bioStart();
				} else {
					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					stopCodonEnd = lastFeature.location().bioEnd();
				}
			}
			//for gtf ordering can be strand based so first is last and last is first
			if (startCodonBegin > stopCodonEnd) {
				int temp = startCodonBegin;
				startCodonBegin = stopCodonEnd;
				stopCodonEnd = temp;
			}



			AccessionID transcriptAccessionID = new AccessionID(mRNAID);
			geneSequence = seq.getGene(geneID);
			if (geneSequence == null) {
				geneSequence = seq.addGene(geneAccessionID, startCodonBegin, stopCodonEnd, strand);
				geneSequence.setSource(geneSource);
				if (geneNote != null && geneNote.length() > 0) {
					geneSequence.addNote(geneNote);
				}
			} else {

				if (startCodonBegin < geneSequence.getBioBegin()) {
					geneSequence.setBioBegin(startCodonBegin);
				}
				if (stopCodonEnd > geneSequence.getBioBegin()) {
					geneSequence.setBioEnd(stopCodonEnd);
				}

			}
			TranscriptSequence transcriptSequence = geneSequence.addTranscript(transcriptAccessionID, startCodonBegin, stopCodonEnd);
			transcriptSequence.setSource(mRNAsource);
			if (mRNANote != null && mRNANote.length() > 0) {
				transcriptSequence.addNote(mRNANote);

			}
			if (startCodon != null) {
				if (startCodonName == null || startCodonName.length() == 0) {
					startCodonName = mRNAID + "-start_codon-" + startCodon.location().bioStart() + "-" + startCodon.location().bioEnd();
				}
				transcriptSequence.addStartCodonSequence(new AccessionID(startCodonName), startCodon.location().bioStart(), startCodon.location().bioEnd());
			}
			if (stopCodon != null) {
				if (stopCodonName == null || stopCodonName.length() == 0) {
					stopCodonName = mRNAID + "-stop_codon-" + stopCodon.location().bioStart() + "-" + stopCodon.location().bioEnd();
				}
				transcriptSequence.addStopCodonSequence(new AccessionID(stopCodonName), stopCodon.location().bioStart(), stopCodon.location().bioEnd());
			}

			for (FeatureI cdsFeature : cdsFeatures) {
				Feature cds = (Feature) cdsFeature;
				String cdsNote = cdsFeature.getAttribute("Note");
				String cdsSource = cds.source();
				String cdsName = cds.getAttribute("ID");
				if (cdsName == null || cdsName.length() == 0) {
					cdsName = mRNAID + "-cds-" + cds.location().bioStart() + "-" + cds.location().bioEnd();
				}
				AccessionID cdsAccessionID = new AccessionID(cdsName);
				ExonSequence exonSequence = geneSequence.addExon(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd());
				exonSequence.setSource(cdsSource);
				if (cdsNote != null && cdsNote.length() > 0) {
					exonSequence.addNote(cdsNote);
				}
				transcriptSequence.addCDS(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd(), cds.frame());
			}
			geneSequence.addIntronsUsingExons();

		}

	}
1162741
def addPhenotypeSearchOptions(parser):
    
    parser.add_argument(
        "--phenotype_association_set_id", "-s", default=None,
        help="Only return phenotypes from this phenotype_association_set.")
    parser.add_argument(
        "--phenotype_id", "-p", default=None,
        help="Only return this phenotype.")
    parser.add_argument(
        "--description", "-d", default=None,
        help="Only return phenotypes matching this description.")
    parser.add_argument(
        "--age_of_onset", "-a", default=None,
        help="Only return phenotypes with this age_of_onset.")
    parser.add_argument(
        "--type", "-T", default=None,
        help="Only return phenotypes with this type.")
774878
def finalize_vcf(in_file, variantcaller, items):
    
    out_file = "%s-annotated%s" % utils.splitext_plus(in_file)
    if not utils.file_uptodate(out_file, in_file):
        header_cl = _add_vcf_header_sample_cl(in_file, items, out_file)
        contig_cl = _add_contig_cl(in_file, items, out_file)
        cls = [x for x in (contig_cl, header_cl) if x]
        if cls:
            post_cl = " | ".join(cls) + " | "
        else:
            post_cl = None
        dbsnp_file = tz.get_in(("genome_resources", "variation", "dbsnp"), items[0])
        if dbsnp_file:
            out_file = _add_dbsnp(in_file, dbsnp_file, items[0], out_file, post_cl)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, items[0]["config"])
    else:
        return in_file
780376
def _gaFeatureForFeatureDbRecord(self, feature):
        
        gaFeature = protocol.Feature()
        gaFeature.id = self.getCompoundIdForFeatureId(feature['id'])
        if feature.get('parent_id'):
            gaFeature.parent_id = self.getCompoundIdForFeatureId(
                    feature['parent_id'])
        else:
            gaFeature.parent_id = ""
        gaFeature.feature_set_id = self.getId()
        gaFeature.reference_name = pb.string(feature.get('reference_name'))
        gaFeature.start = pb.int(feature.get('start'))
        gaFeature.end = pb.int(feature.get('end'))
        gaFeature.name = pb.string(feature.get('name'))
        if feature.get('strand', '') == '-':
            gaFeature.strand = protocol.NEG_STRAND
        else:
            # default to positive strand
            gaFeature.strand = protocol.POS_STRAND
        gaFeature.child_ids.extend(map(
                self.getCompoundIdForFeatureId,
                json.loads(feature['child_ids'])))
        gaFeature.feature_type.CopyFrom(
            self._ontology.getGaTermByName(feature['type']))
        attributes = json.loads(feature['attributes'])
        # TODO: Identify which values are ExternalIdentifiers and OntologyTerms
        for key in attributes:
            for v in attributes[key]:
                gaFeature.attributes.attr[key].values.add().string_value = v
        if 'gene_name' in attributes and len(attributes['gene_name']) > 0:
            gaFeature.gene_symbol = pb.string(attributes['gene_name'][0])
        return gaFeature
329389
def _parse_patient_variants(self, file):
        
        patient_var_map = self._convert_variant_file_to_dict(file)
        gene_coordinate_map = self._parse_gene_coordinates(
            self.map_files['gene_coord_map'])
        rs_map = self._parse_rs_map_file(self.map_files['dbsnp_map'])

        genotype = Genotype(self.graph)
        model = Model(self.graph)

        self._add_variant_gene_relationship(patient_var_map, gene_coordinate_map)

        for patient in patient_var_map:
            patient_curie = ':{0}'.format(patient)
            # make intrinsic genotype for each patient
            intrinsic_geno_bnode = self.make_id(
                "{0}-intrinsic-genotype".format(patient), "_")
            genotype_label = "{0} genotype".format(patient)
            genotype.addGenotype(
                intrinsic_geno_bnode, genotype_label,
                model.globaltt['intrinsic_genotype'])

            self.graph.addTriple(
                patient_curie, model.globaltt['has_genotype'], intrinsic_geno_bnode)
            for variant_id, variant in patient_var_map[patient].items():
                build = variant['build']
                chromosome = variant['chromosome']
                position = variant['position']
                reference_allele = variant['reference_allele']
                variant_allele = variant['variant_allele']
                genes_of_interest = variant['genes_of_interest']
                rs_id = variant['rs_id']

                variant_label = ''
                variant_bnode = self.make_id("{0}".format(variant_id), "_")

                # maybe should have these look like the elif statements below
                if position and reference_allele and variant_allele:
                    variant_label = self._build_variant_label(
                        build, chromosome, position, reference_allele,
                        variant_allele, genes_of_interest)
                elif not position and reference_allele and variant_allele \
                        and len(genes_of_interest) == 1:

                    variant_label = self._build_variant_label(
                        build, chromosome, position, reference_allele, variant_allele,
                        genes_of_interest)
                elif position and (not reference_allele or not variant_allele) \
                        and len(genes_of_interest) == 1:

                    variant_label = "{0}{1}({2}):g.{3}".format(
                        build, chromosome, genes_of_interest[0], position)
                elif len(genes_of_interest) == 1:
                    variant_label = 'variant of interest in {0} gene of patient' \
                        ' {1}'.format(genes_of_interest[0], patient)
                else:
                    variant_label = 'variant of interest in patient {0}'.format(patient)

                genotype.addSequenceAlteration(variant_bnode, None)
                # check if it we have built the label
                # in _add_variant_gene_relationship()
                labels = self.graph.objects(
                    BNode(re.sub(r'^_:', '', variant_bnode, 1)), RDFS['label'])

                label_list = list(labels)

                if len(label_list) == 0:
                    model.addLabel(variant_bnode, variant_label)

                self.graph.addTriple(
                    variant_bnode, self.globaltt['in taxon'],
                    self.globaltt['Homo sapiens'])
                self.graph.addTriple(
                    intrinsic_geno_bnode, self.globaltt['has_variant_part'],
                    variant_bnode)
                if rs_id:
                    dbsnp_curie = 'dbSNP:{0}'.format(rs_id)
                    model.addSameIndividual(variant_bnode, dbsnp_curie)

        self._add_variant_sameas_relationships(patient_var_map, rs_map)
        return
275822
def search_genotype_phenotype(
            self, phenotype_association_set_id=None, feature_ids=None,
            phenotype_ids=None, evidence=None):
        
        request = protocol.SearchGenotypePhenotypeRequest()
        request.phenotype_association_set_id = phenotype_association_set_id
        if feature_ids:
            request.feature_ids.extend(feature_ids)
        if phenotype_ids:
            request.phenotype_ids.extend(phenotype_ids)
        if evidence:
            request.evidence.extend(evidence)
        request.page_size = pb.int(self._page_size)
        self._logger.debug("search_genotype_phenotype {}".format(request))
        return self._run_search_request(
            request, "featurephenotypeassociations",
            protocol.SearchGenotypePhenotypeResponse)
395630
def extract_germline_vcinfo(data, out_dir):
    
    supported_germline = set(["vardict", "octopus", "freebayes"])
    if dd.get_phenotype(data) in ["tumor"]:
        for v in _get_variants(data):
            if v.get("variantcaller") in supported_germline:
                if v.get("germline"):
                    return v
                else:
                    d = utils.deepish_copy(data)
                    d["vrn_file"] = v["vrn_file"]
                    gd = germline.extract(d, [d], out_dir)
                    v["germline"] = gd["vrn_file_plus"]["germline"]
                    return v
371118
static public void addGeneIDGFF2GeneFeatures(LinkedHashMap<String, ChromosomeSequence> chromosomeSequenceList, FeatureList listGenes) throws Exception {
		Collection<String> geneIds = listGenes.attributeValues("gene_id");
		for (String geneid : geneIds) {
			FeatureList gene = listGenes.selectByAttribute("gene_id", geneid);
			FeatureI geneFeature = gene.get(0);
			ChromosomeSequence seq = chromosomeSequenceList.get(geneFeature.seqname());
			geneid = geneid.replaceAll("_", ".G");
			AccessionID geneAccessionID = new AccessionID(geneid);
			GeneSequence geneSequence = null;
			Collection<String> transcriptids = gene.attributeValues("gene_id");
			for (String transcriptid : transcriptids) {
				// get all the individual features (exons, CDS regions, etc.) of this gene
				FeatureList transcriptFeature = listGenes.selectByAttribute("gene_id", transcriptid);
				transcriptid = transcriptid.replaceAll("_", ".G");




				//      String seqName = feature.seqname();
				//FeatureI startCodon = null;
				//FeatureI stopCodon = null;
				Integer startCodonBegin = null;
				Integer stopCodonEnd = null;
				//String startCodonName = "";
				//String stopCodonName = "";


				// now select only the coding regions of this gene
				FeatureList firstFeatures = transcriptFeature.selectByType("First");
				FeatureList terminalFeatures = transcriptFeature.selectByType("Terminal");
				FeatureList internalFeatures = transcriptFeature.selectByType("Internal");
				FeatureList singleFeatures = transcriptFeature.selectByType("Single");
				FeatureList cdsFeatures = new FeatureList();
				cdsFeatures.add(firstFeatures);
				cdsFeatures.add(terminalFeatures);
				cdsFeatures.add(internalFeatures);
				cdsFeatures.add(singleFeatures);
				// sort them
				cdsFeatures = cdsFeatures.sortByStart();
				Strand strand = Strand.POSITIVE;
				FeatureI feature = cdsFeatures.get(0);
				if (feature.location().isNegative()) {
					strand = Strand.NEGATIVE;
				}
				if (startCodonBegin == null) {
					FeatureI firstFeature = cdsFeatures.get(0);
					if (strand == Strand.NEGATIVE) {
						startCodonBegin = firstFeature.location().bioEnd();
					} else {
						startCodonBegin = firstFeature.location().bioStart();
					}
				}

				if (stopCodonEnd == null) {

					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					if (strand == Strand.NEGATIVE) {
						stopCodonEnd = lastFeature.location().bioStart();
					} else {
						stopCodonEnd = lastFeature.location().bioEnd();
					}
				}
				//for gtf ordering can be strand based so first is last and last is first
				if (startCodonBegin > stopCodonEnd) {
					int temp = startCodonBegin;
					startCodonBegin = stopCodonEnd;
					stopCodonEnd = temp;
				}

				AccessionID transcriptAccessionID = new AccessionID(transcriptid);
				if (geneSequence == null) {
					geneSequence = seq.addGene(geneAccessionID, startCodonBegin, stopCodonEnd, strand);
					geneSequence.setSource(((Feature) feature).source());
				} else {
					//if multiple transcripts for one gene make sure the gene is defined as the min and max start/end

					if (startCodonBegin < geneSequence.getBioBegin()) {
						geneSequence.setBioBegin(startCodonBegin);
					}
					if (stopCodonEnd > geneSequence.getBioBegin()) {
						geneSequence.setBioEnd(stopCodonEnd);
					}

				}
				TranscriptSequence transcriptSequence = geneSequence.addTranscript(transcriptAccessionID, startCodonBegin, stopCodonEnd);
				/*
				if (startCodon != null) {
					if (startCodonName == null || startCodonName.length() == 0) {
						startCodonName = transcriptid + "-start_codon-" + startCodon.location().bioStart() + "-" + startCodon.location().bioEnd();
					}
					transcriptSequence.addStartCodonSequence(new AccessionID(startCodonName), startCodon.location().bioStart(), startCodon.location().bioEnd());
				}
				if (stopCodon != null) {
					if (stopCodonName == null || stopCodonName.length() == 0) {
						stopCodonName = transcriptid + "-stop_codon-" + stopCodon.location().bioStart() + "-" + stopCodon.location().bioEnd();
					}
					transcriptSequence.addStopCodonSequence(new AccessionID(stopCodonName), stopCodon.location().bioStart(), stopCodon.location().bioEnd());
				}
				*/

				for (FeatureI cdsFeature : cdsFeatures) {
					Feature cds = (Feature) cdsFeature;
					String cdsName = cds.getAttribute("transcript_name");
					if (cdsName == null || cdsName.length() == 0) {
						cdsName = transcriptid + "-cds-" + cds.location().bioStart() + "-" + cds.location().bioEnd();
					}
					AccessionID cdsAccessionID = new AccessionID(cdsName);
					//ExonSequence exonSequence = geneSequence.addExon(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd());
					CDSSequence cdsSequence = transcriptSequence.addCDS(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd(), cds.frame());
					cdsSequence.setSequenceScore(cds.score());
				}
			}
		}

	}
535424
def run(vrn_info, cnvs_by_name, somatic_info):
    
    config = {"sample_size": 5000}
    work_dir = _cur_workdir(somatic_info.tumor_data)
    if "battenberg" not in cnvs_by_name:
        logger.warn("PhyloWGS requires Battenberg CNV calls, skipping %s"
                    % dd.get_sample_name(somatic_info.tumor_data))
    else:
        ssm_file, cnv_file = _prep_inputs(vrn_info, cnvs_by_name["battenberg"], somatic_info, work_dir, config)
        evolve_file = _run_evolve(ssm_file, cnv_file, work_dir, somatic_info.tumor_data)
        summary_file = _prepare_summary(evolve_file, ssm_file, cnv_file, work_dir, somatic_info)
        print(summary_file)
620417
def add_flanking_seqs(self, ref_seq, new_start, new_end):
        ''''''
        if new_start > self.POS or new_end < self.ref_end_pos():
            raise Error('new start and end positions must not try to shrink VCF record. new_start=' + str(new_start) + ', new_end=' + str(new_end) + '. VCF=' + str(self))

        new_start_nucleotides = ref_seq[new_start:self.POS]
        new_end_nucleotodes = ref_seq[self.ref_end_pos() + 1:new_end + 1]
        self.POS = new_start
        self.REF = new_start_nucleotides + self.REF + new_end_nucleotodes
        self.ALT = [new_start_nucleotides + x + new_end_nucleotodes for x in self.ALT]
1336636
def add_subject(self, ind_id, sex=None, phenotype=None):
        

        self.pedigree_data[ind_id] = len(self.phenotype_data[0])
        if phenotype != None:
            if type(self.phenotype_data) is list:
                self.phenotype_data[0].append(phenotype)
            else:
                self.phenotype_data[-1, len(self.individual_mask)] = phenotype
        self.individual_mask.append(0)

        if PhenoCovar.sex_as_covariate:
            try:
                self.covariate_data[0].append(float(sex))
            except Exception, e:
                raise MalformedInputFile("Invalid setting, %s, for sex in pedigree" % (sex))