Wrapper for _load_sequences_to_reference_gene
176247
def load_sequences_to_reference(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _load_sequences_to_reference_gene_sc(g_id, outdir=self.sequences_by_gene_dir,
                                                 g_to_pickle=self.gene_protein_pickles,
                                                 strain_infodict=self.strain_infodict,
                                                 orth_matrix=self.df_orthology_matrix, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            from Bio import SeqIO
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs.pckl'.format(g_id))

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

                for strain, info in strain_infodict.items():
                    strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                    strain_gene_functional = info['functional_genes'][g_id]
                    if strain_gene_functional:
                        # Pull the gene ID of the strain from the orthology matrix
                        strain_gene_key = orth_matrix.at[g_id, strain]
                        new_id = '{}_{}'.format(g_id, strain)
                        if protein_pickle.sequences.has_id(new_id):
                            continue
                        protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                            ident=new_id,
                                                            set_as_representative=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_load_sequences_to_reference_gene_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        updated = []
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
            updated.append(g_id)
        not_updated = set(list(self.gene_protein_pickles.keys())).difference(updated)
        log.info('No change to {} genes, removing from gene_protein_pickles'.format(len(not_updated)))
        log.debug(not_updated)
        for rem in not_updated:
            del self.gene_protein_pickles[rem]
161104
def read_10xgenomics(cls, tarball_fpath: str, prefix: str,
                         use_ensembl_ids: bool = False):
        """"""

        _LOGGER.info('Reading file: %s', tarball_fpath)

        with tarfile.open(tarball_fpath, mode='r:gz') as tf:
            ti = tf.getmember('%smatrix.mtx' % prefix)
            with tf.extractfile(ti) as fh:
                mtx = scipy.io.mmread(fh)

            ti = tf.getmember('%sgenes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                i = 1
                if use_ensembl_ids:
                    i = 0
                gene_names = \
                        [row[i] for row in csv.reader(wrapper, delimiter='\t')]

            ti = tf.getmember('%sbarcodes.tsv' % prefix)
            with tf.extractfile(ti) as fh:
                wrapper = io.TextIOWrapper(fh, encoding='ascii')
                barcodes = \
                        [row[0] for row in csv.reader(wrapper, delimiter='\t')]


            assert mtx.shape[0] == len(gene_names)
            assert mtx.shape[1] == len(barcodes)
        
        _LOGGER.info('Matrix dimensions: %s', str(mtx.shape))
        X = mtx.todense()
        matrix = cls(X=X, genes=gene_names, cells=barcodes)
        
        return matrix
1214775
def get_transcript(self, gene_pk, refseq_id):
        ""
        if not refseq_id:
            return
        transcript_pk = self.transcripts.get(refseq_id)
        if transcript_pk:
            return transcript_pk
        gene = Gene(pk=gene_pk)
        transcript = Transcript(refseq_id=refseq_id, gene=gene)
        try:
            transcript.save()
        except IntegrityError:
            transcript = Transcript.objects.get(refseq_id=refseq_id, gene=gene)
        self.transcripts[refseq_id] = transcript.pk
        return transcript.pk
859388
def consume_seqreg(self, seqreg):
        """"""
        if not isinstance(seqreg, tag.directive.Directive) or \
                seqreg.type != 'sequence-region':
            raise ValueError('expected ##sequence-region directive')
        if seqreg.seqid in self.declared_regions:
            msg = 'duplicate sequence region "{}"'.format(seqreg.seqid)
            raise ValueError(msg)
        self.declared_regions[seqreg.seqid] = seqreg.range.copy()
589388
def load_file(self, sequence):
        """"""
        try:
            if sequence.filetype_ext == 'npy':
                sequence.series = sequence.adjust_series(
                    *self._load_npy(sequence))
            elif sequence.filetype_ext == 'asc':
                sequence.series = sequence.adjust_series(
                    *self._load_asc(sequence))
            elif sequence.filetype_ext == 'nc':
                self._load_nc(sequence)
        except BaseException:
            objecttools.augment_excmessage(
                'While trying to load the external data of sequence %s'
                % objecttools.devicephrase(sequence))
266347
def parse(handle, seqtype=None, robust=False):
    ''''''
    if seqtype is None:
        seqtype = _get_seqtype_from_ext(handle)

    if seqtype.startswith('gz-'):
        handle = _unzip_handle(handle)
        seqtype = seqtype[3:]

    # False positive from pylint, both handles are fileobj-like
    # pylint: disable=redefined-variable-type
    if robust:
        if seqtype == "embl":
            handle = sanity_check_embl(handle)
        elif seqtype == "genbank":
            handle = sanity_check_genbank(handle)
        elif seqtype == "fasta":
            handle = sanity_check_fasta(handle)
    # pylint: enable=redefined-variable-type

    return SeqIO.parse(handle, seqtype)
31071
def from_seqfeature(s, **kwargs):
    
    source = s.qualifiers.get('source', '.')[0]
    score = s.qualifiers.get('score', '.')[0]
    seqid = s.qualifiers.get('seqid', '.')[0]
    frame = s.qualifiers.get('frame', '.')[0]
    strand = _feature_strand[s.strand]

    # BioPython parses 1-based GenBank positions into 0-based for use within
    # Python.  We need to convert back to 1-based GFF format here.
    start = s.location.start.position + 1
    stop = s.location.end.position
    featuretype = s.type
    id = s.id
    attributes = dict(s.qualifiers)
    attributes.pop('source', '.')
    attributes.pop('score', '.')
    attributes.pop('seqid', '.')
    attributes.pop('frame', '.')
    return Feature(seqid, source, featuretype, start, stop, score, strand,
                   frame, attributes, id=id, **kwargs)
257105
def get_sequence(self,ref):
    """"""
    strand = '+'
    if not self._options.direction:
      sys.stderr.write("WARNING: no strand information for the transcript\n")
    if self._options.direction: strand = self._options.direction
    seq = ''
    for e in [x.range for x in self.exons]:
      seq += str(ref[e.chr][e.start-1:e.end])
    if strand == '-':  seq = rc(seq)
    return Sequence(seq.upper())
811157
def parse_codons(ref, start, end, strand):
    
    codon = []
    c = cycle([1, 2, 3])
    ref = ref[start - 1:end]
    if strand == -1:
        ref = rc_stats(ref)
    for pos in ref:
        n = next(c)
        codon.append(pos)
        if n == 3:
            yield codon
            codon = []
1078224
def _load_data_index(self, chain1, chain2):
        
        parser = PDB.PDBParser(QUIET=True)

        structure1 = parser.get_structure(chain1, self.pdb1)
        structure2 = parser.get_structure(chain2, self.pdb2)

        residues1 = list(structure1.get_residues())
        residues2 = list(structure2.get_residues())

        indexes1 = set(r.id[1] for r in residues1)
        indexes2 = set(r.id[1] for r in residues2)

        indexes = indexes1.intersection(indexes2)
        self.indexes = indexes.copy()
        self.N = len(indexes)

        coord1 = []
        indexes1 = indexes.copy()
        for r in residues1:
            if r.id[1] in indexes1 and 'CA' in r:
                coord1.append(np.concatenate((r['CA'].get_coord(), (1,)))[:, None])
                # Remove from index to avoid repeated residues
                indexes1.remove(r.id[1])
        coord1 = np.hstack(coord1).astype(DTYPE, copy=False)

        coord2 = []
        for r in residues2:
            if r.id[1] in indexes and 'CA' in r:
                coord2.append(np.concatenate((r['CA'].get_coord(), (1,)))[:, None])
                indexes.remove(r.id[1])
        coord2 = np.hstack(coord2).astype(DTYPE, copy=False)

        self.coord1 = coord1
        self.coord2 = coord2
85137
def load_just_in_time_genome_alignment(path, ref_spec, extensions=None,
                                       index_exts=None, fail_no_index=True,
                                       verbose=False):
    """"""
    if index_exts is None and fail_no_index:
      raise ValueError("Failure on no index specified for loading genome " +
                       "alignment, but no index extensions specified")

    extensions = __trim_extensions_dot(extensions)
    index_exts = __trim_extensions_dot(index_exts)

    partial_chrom_files = {}
    whole_chrom_files = {}
    for fn in os.listdir(path):
      pth = os.path.join(path, fn)
      if os.path.isfile(pth):
        base, ext = os.path.splitext(pth)
        ext = __trim_extension_dot(ext)
        if extensions is None or ext in extensions:
          idx_path = __find_index(pth, index_exts)
          if idx_path is None and fail_no_index:
            raise PyokitIOError("No index file for " + fn)
          chrom, start, end = __split_genomic_interval_filename(fn)
          assert((start is None and end is None) or
                 (start is not None and end is not None))
          if start is None:
            if chrom in whole_chrom_files:
              raise PyokitIOError("multiple files for chrom " + chrom)
            whole_chrom_files[chrom] = (pth, idx_path)
          else:
            k = (chrom, start, end)
            if k in partial_chrom_files:
              other_pth = partial_chrom_files[k][0]
              raise PyokitIOError("multiple files for " + str(k) + " --> " +
                                  pth + " and " + other_pth)
            partial_chrom_files[k] = (pth, idx_path)

    def factory(k):
      pth, idx = k
      return build_genome_alignment_from_file(pth, ref_spec, idx, verbose)

    return JustInTimeGenomeAlignment(whole_chrom_files, partial_chrom_files,
                                     factory)
1159976
def liftover_to_genome(pass_pos, gtf):
    """"""

    fixed_pos = []
    for pos in pass_pos:
        if pos["chrom"] not in gtf:
            continue
        db_pos = gtf[pos["chrom"]][0]
        mut = _parse_mut(pos["sv"])
        print([db_pos, pos])
        if db_pos[3] == "+":
            pos['pre_pos'] = db_pos[1] + pos["pre_pos"] + 1
        else:
            pos['pre_pos'] = db_pos[2] - (pos["pre_pos"] - 1)
        pos['chrom'] = db_pos[0]
        pos['nt'] = list(mut[0])
        fixed_pos.append(pos)

    _print_header(fixed_pos)
    for pos in fixed_pos:
        print_vcf(pos)
43344
def get_resources(genome, ref_file, data):
    
    base_dir = os.path.normpath(os.path.dirname(ref_file))
    resource_file = os.path.join(base_dir, "%s-resources.yaml" % genome.replace("-test", ""))
    if not os.path.exists(resource_file):
        raise IOError("Did not find resource file for %s: %s\n"
                      "To update bcbio_nextgen.py with genome resources for standard builds, run:\n"
                      "bcbio_nextgen.py upgrade -u skip"
                      % (genome, resource_file))
    with open(resource_file) as in_handle:
        resources = yaml.safe_load(in_handle)

    def resource_file_path(x):
        if isinstance(x, six.string_types) and os.path.exists(os.path.join(base_dir, x)):
            return os.path.normpath(os.path.join(base_dir, x))
        return x
    cleaned = utils.dictapply(resources, resource_file_path)
    return ensure_annotations(cleaned, data)
821106
def reprofiler(self, header, genome, sample):
        
        # Iterate through mlstseqtype - it contains genomes with partial matches to current reference profiles
        # Reset :newprofile
        newprofile = ""
        # Find the last profile entry in the dictionary of profiles
        # Opens uses the command line tool 'tail' to look at the last line of the file (-1). This last line
        # is split on tabs, and only the first entry (the sequence type number) is captured
        if sample[self.analysistype].supplementalprofile != 'NA':
            if os.path.isfile(sample[self.analysistype].supplementalprofile):
                try:
                    lastentry = int(
                        subprocess.check_output(['tail', '-1', sample[self.analysistype].supplementalprofile])
                        .split("\t")[0]) + 1
                except ValueError:
                    lastentry = 1000000
            else:
                open(sample[self.analysistype].supplementalprofile, 'w').close()
                lastentry = 1000000
            # As there can be multiple profiles in MLSTSeqType, this loop only needs to be performed once.
            seqcount = 0
            # Go through the sequence types
            try:
                sequencetype = list(self.mlstseqtype[genome].keys())[0]
            except IndexError:
                sequencetype = ''
                seqcount = 1
            # Only do this once
            if seqcount == 0:
                # Set the :newprofile string to start with the new profile name (e.g. 1000000_CFIA)
                newprofile = str(lastentry)
                # The number of matches to the reference profile
                nummatches = list(self.mlstseqtype[genome][sequencetype].keys())[0]
                # The genes in geneList - should be in the correct order
                for gene in sorted(sample[self.analysistype].allelenames):
                    # The allele for each gene in the query genome
                    allele = list(self.mlstseqtype[genome][sequencetype][nummatches][gene].keys())[0]
                    # Append the allele to newprofile
                    newprofile += '\t{}'.format(allele)
                    # Add the MLST results for the query genome as well as the new profile data
                    # to resultProfile
                    self.resultprofile[genome]['{}(new)'.format(str(lastentry))][header][gene][allele] = \
                        list(self.mlstseqtype[genome][sequencetype][nummatches][gene][allele].values())[0]
                seqcount += 1
                sample[self.analysistype].mismatchestosequencetype = 'NA'
                sample[self.analysistype].matchestosequencetype = header
            # Only perform the next loop if :newprofile exists
            if newprofile:
                # Open the profile file to append
                with open(sample[self.analysistype].supplementalprofile, 'a') as appendfile:
                    # Append the new profile to the end of the profile file
                    appendfile.write('{}\n'.format(newprofile))
                # Re-run profiler with the updated files
                self.profiler()
        else:
            sample[self.analysistype].mismatchestosequencetype = 'NA'
            sample[self.analysistype].matchestosequencetype = 'NA'
809486
def _get_genes(self, variant):
        
        transcripts = variant['transcripts']
        ensembl_ids = [transcript['ensembl_id'] for transcript in
                       transcripts if transcript['ensembl_id']]
        hgnc_symbols = [transcript['hgnc_symbol'] for transcript in
                        transcripts if transcript['hgnc_symbol']]
        genes = get_gene_info(ensembl_ids, hgnc_symbols)
        return genes
113076
def get_protein_coding_genes(
        path_or_buffer, 
        include_polymorphic_pseudogenes=True,
        remove_duplicates=True,
        **kwargs):
    

    """
    valid_biotypes = set(['protein_coding'])
    if include_polymorphic_pseudogenes:
        valid_biotypes.add('polymorphic_pseudogene')
    
    df = get_genes(path_or_buffer, valid_biotypes,
                   remove_duplicates=remove_duplicates, **kwargs)
    return df
604216
def insert_sequences_into_tree(aln, moltype, params={}):
    
    # convert aln to phy since seq_names need fixed to run through parsinsert
    new_aln=get_align_for_phylip(StringIO(aln))

    # convert aln to fasta in case it is not already a fasta file
    aln2 = Alignment(new_aln)
    seqs = aln2.toFasta()

    parsinsert_app = ParsInsert(params=params)
    result = parsinsert_app(seqs)

    # parse tree
    tree = DndParser(result['Tree'].read(), constructor=PhyloNode)

    # cleanup files
    result.cleanUp()

    return tree