list: List of genes with no mapping to a representative sequence.
330631
def missing_representative_sequence(self):
        """"""
        return [x.id for x in self.genes if not self.genes_with_a_representative_sequence.has_id(x.id)]
1270789
def group_unique_sequences(
        allele_reads,
        max_prefix_size=None,
        max_suffix_size=None):
    
    groups = defaultdict(set)
    for r in allele_reads:
        prefix = r.prefix
        allele = r.allele
        suffix = r.suffix
        if max_prefix_size and len(prefix) > max_prefix_size:
            prefix = prefix[-max_prefix_size:]
        if max_suffix_size and len(suffix) > max_suffix_size:
            suffix = suffix[:max_suffix_size]
        key = (prefix, allele, suffix)
        groups[key].add(r)
    return groups
1120517
def map_to_precursor_biopython(seqs, names, loci, args):
    """"""
    precursor = precursor_sequence(loci, args.ref).upper()
    dat = dict()
    for s, n in itertools.izip(seqs, names):
        res = _align(str(s), precursor)
        if res:
            dat[n] = res
    logger.debug("mapped in %s: %s out of %s" % (loci, len(dat), len(seqs)))
    return dat
57053
def _align_orthologous_gene_pairwise(self, g_id, gapopen=10, gapextend=0.5, engine='needle', parse=True, force_rerun=False):
        """"""
        protein_seqs_aln_pickle_path = op.join(self.sequences_by_gene_dir, '{}_protein_withseqs_dis_aln.pckl'.format(g_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_aln_pickle_path):
            protein_seqs_pickle_path = self.gene_protein_pickles[g_id]
            protein_pickle = ssbio.io.load_pickle(protein_seqs_pickle_path)

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

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

            alignment_dir = op.join(self.sequences_by_gene_dir, g_id)
            ssbio.utils.make_dir(alignment_dir)
            protein_pickle.pairwise_align_sequences_to_representative(gapopen=gapopen, gapextend=gapextend,
                                                                      engine=engine, outdir=alignment_dir,
                                                                      parse=parse, force_rerun=force_rerun)
            protein_pickle.save_pickle(outfile=protein_seqs_aln_pickle_path)

        return g_id, protein_seqs_aln_pickle_path
1242986
def _add_genes_to_bed(in_file, gene_file, fai_file, out_file, data, max_distance=10000):
    
    try:
        input_rec = next(iter(pybedtools.BedTool(in_file)))
    except StopIteration:  # empty file
        utils.copy_plus(in_file, out_file)
        return
    # keep everything after standard chrom/start/end, 1-based
    extra_fields = list(range(4, len(input_rec.fields) + 1))
    # keep the new gene annotation
    gene_index = len(input_rec.fields) + 4
    extra_fields.append(gene_index)
    columns = ",".join([str(x) for x in extra_fields])
    max_column = max(extra_fields) + 1
    ops = ",".join(["distinct"] * len(extra_fields))
    # swap over gene name to '.' if beyond maximum distance
    # cut removes the last distance column which can cause issues
    # with bedtools merge: 'ERROR: illegal character '.' found in integer conversion of string'
    distance_filter = (r"""awk -F$'\t' -v OFS='\t' '{if ($NF > %s || $NF < -%s) $%s = "."} {print}'""" %
                       (max_distance, max_distance, gene_index))
    sort_cmd = bedutils.get_sort_cmd(os.path.dirname(out_file))
    cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
    # Ensure gene transcripts match reference genome
    ready_gene_file = os.path.join(os.path.dirname(out_file), "%s-genomeonly.bed" %
                                   (utils.splitext_plus(os.path.basename(gene_file))[0]))
    ready_gene_file = bedutils.subset_to_genome(gene_file, ready_gene_file, data)
    exports = "export TMPDIR=%s && %s" % (os.path.dirname(out_file), utils.local_path_export())
    bcbio_py = sys.executable
    gsort = config_utils.get_program("gsort", data)
    cmd = ("{exports}{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^# | "
           "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
           "{gsort} - {fai_file} | "
            "bedtools closest -g {fai_file} "
            "-D ref -t first -a - -b <({gsort} {ready_gene_file} {fai_file}) | "
            "{distance_filter} | cut -f 1-{max_column} | "
            "bedtools merge -i - -c {columns} -o {ops} -delim ',' -d -10 > {out_file}")
    do.run(cmd.format(**locals()), "Annotate BED file with gene info")
194455
def df_representative_structures(self):
        """"""
        rep_struct_pre_df = []
        df_cols = ['gene', 'id', 'is_experimental', 'file_type', 'structure_file']

        for g in self.genes_with_a_representative_structure:
            repdict = g.protein.representative_structure.get_dict(df_format=True, only_attributes=df_cols)
            repdict['gene'] = g.id
            rep_struct_pre_df.append(repdict)

        df = pd.DataFrame.from_records(rep_struct_pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)
402292
def gene_one_per_line(in_file, data):
    
    if in_file:
        # Report all duplicated annotations one-per-line
        one_per_line_file = "%s-opl.bed" % utils.splitext_plus(in_file)[0]
        if not utils.file_uptodate(one_per_line_file, in_file):
            with file_transaction(data, one_per_line_file) as tx_out_file:
                with open(tx_out_file, 'w') as out:
                    for r in pybedtools.BedTool(in_file):
                        for g in r.name.split(','):
                            out.write('\t'.join(map(str, [r.chrom, r.start, r.end, g])) + '\n')
        return one_per_line_file
782965
def genes(self):
        """"""
        return [ExpGene.from_series(g)
                for i, g in self.reset_index().iterrows()]
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
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
1231139
def add_genes(in_file, data, max_distance=10000, work_dir=None):
    
    gene_file = regions.get_sv_bed(data, "exons", out_dir=os.path.dirname(in_file))
    if gene_file and utils.file_exists(in_file):
        out_file = "%s-annotated.bed" % utils.splitext_plus(in_file)[0]
        if work_dir:
            out_file = os.path.join(work_dir, os.path.basename(out_file))
        if not utils.file_uptodate(out_file, in_file):
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            with file_transaction(data, out_file) as tx_out_file:
                _add_genes_to_bed(in_file, gene_file, fai_file, tx_out_file, data, max_distance)
        return out_file
    else:
        return in_file
389135
def _write_strain_functional_genes(self, strain_id, ref_functional_genes, orth_matrix, force_rerun=False):
        """"""
        func_genes_path = op.join(self.model_dir, '{}_funcgenes.json'.format(strain_id))

        if ssbio.utils.force_rerun(flag=force_rerun, outfile=func_genes_path):
            gene_to_func = {k:True for k in ref_functional_genes}
            # Get a list of genes which do not have orthology in the strain
            genes_to_remove = orth_matrix[pd.isnull(orth_matrix[strain_id])][strain_id].index.tolist()

            # Mark genes non-functional
            genes_to_remove = list(set(genes_to_remove).intersection(set(ref_functional_genes)))

            if len(genes_to_remove) > 0:
                for g in genes_to_remove:
                    gene_to_func[g] = False

            with open(func_genes_path, 'w') as f:
                json.dump(gene_to_func, f)
        else:
            with open(func_genes_path, 'r') as f:
                gene_to_func = json.load(f)

        return strain_id, gene_to_func
645005
def map_to_precursors(seqs, names, loci, out_file, args):
    """"""
    with make_temp_directory() as temp:
        pre_fasta = os.path.join(temp, "pre.fa")
        seqs_fasta = os.path.join(temp, "seqs.fa")
        out_sam = os.path.join(temp, "out.sam")
        pre_fasta = get_loci_fasta(loci, pre_fasta, args.ref)
        out_precursor_file = out_file.replace("tsv", "fa")
        seqs_fasta = get_seqs_fasta(seqs, names, seqs_fasta)

        # print(open(pre_fasta).read().split("\n")[1])
        if find_cmd("razers3"):
            cmd = "razers3 -dr 2 -i 80 -rr 90 -f -o {out_sam} {temp}/pre.fa  {seqs_fasta}"
            run(cmd.format(**locals()))
            out_file = read_alignment(out_sam, loci, seqs, out_file)
            shutil.copy(pre_fasta, out_precursor_file)
    return out_file
51461
def strip_unknown_features(self):
        
        if not self.db:
            return self
        ind = []
        for i, gene_id in enumerate(self.data.index):
            try:
                self.db[gene_id]
                ind.append(i)
            except gffutils.FeatureNotFoundError:
                pass
        ind = np.array(ind)
        return self.__class__(self.data.ix[ind], **self._kwargs)
693187
def find_seq_homologues(self, return_raw=False):
        
        
        
        # Get sequence from structure/chain
        # We could use Bio.PDB.Polypeptide?
        
        from Bio.SCOP.Raf import to_one_letter_code
        
        s = self
        seq_iter = s.get_residues()
        seq_str = ''

        for aa in seq_iter:
            if aa.resname in to_one_letter_code:
                seq_str += to_one_letter_code[aa.resname]
        
        # Use BLAST to find homologous sequences with associated
        # structures in the PDB database.
        # Perhaps include local BLAST?
        
        from Bio.Blast.NCBIWWW import qblast

        # Adapt for short query sequences if needed
        # From http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#shortp
        
        if len(seq_str) < 15:
            word_size = 2
            expect = 20000
            matrix_name = 'PAM30'
            filter = None
        else:
            word_size = 3
            expect = 10.0
            matrix_name = 'BLOSUM62'
            filter = 'SEG'
             
        query_result = qblast(  "blastp", "pdb", seq_str, 
                                word_size, expect, matrix_name, filter)
        
        if return_raw:
           return query_result 
            
        # Parse BLAST result to yield results
        # PDBID : (E-Value, Identity, Positives, Gaps, Alignment)
        
        from Bio.Blast import NCBIXML
        
        blast_records = NCBIXML.read(query_result)
        
        results = []
        
        for alignment in blast_records.alignments:
            for hsps in alignment.hsps:
                id_perc = "%s/%s" %(hsps.identities, alignment.length)
                pos_perc = "%s/%s" %(hsps.positives, alignment.length)
                gaps_perc = "%s/%s" %(hsps.gaps, alignment.length)

                pdb_id = alignment.title.split('|')[3]
                e_value = hsps.expect
                
                results.append( ( pdb_id,
                                "%2.5e" %e_value,
                                id_perc, 
                                pos_perc, 
                                gaps_perc,
                                '\n'.join([hsps.query, hsps.match, hsps.sbjct]) ) )

        return results
1057932
def print_table(seqs, id2name, name):
    
    itable = open('%s.itable' % (name.rsplit('.', 1)[0]), 'w')
    print('\t'.join(['#sequence', 'gene', 'model', 'insertion', 'gene position', 'model position', 'length', 'orf?', 'intron?', 'orf?intron?', 'insertion', 'orf', 'intron']), file=itable)
    for seq, info in list(seqs.items()):
        gene, model, insertions = info
        name = id2name[seq]
        for i, ins in enumerate(insertions, 1):
            gene_pos, model_pos, length, iseq, \
                    orfs, introns, orfs_b, introns_b, orf_annotations = ins
            # append annotation to orf header
            for orf in orfs:
                parts = orf[0].split()
                annotation = orf_annotations[parts[0].split('>')[1]]
                orf[0] = '%s %s %s' % (parts[0], annotation, ' '.join(parts[1:]))
            # get orf position
            gene_pos = '-'.join([str(j) for j in gene_pos])
            # check if orf, intron is present
            if orfs_b is True or introns_b is True:
                orfs_introns_b = True
            else:
                orfs_introns_b = False
            out = [name, gene, model, i, gene_pos, model_pos, length, orfs_b, introns_b, orfs_introns_b]
            out.append('|'.join(iseq))
            out.append('|'.join(['|'.join(orf) for orf in orfs]))
            out.append('|'.join(['|'.join(intron) for intron in introns]))
            print('\t'.join([str(i) for i in out]), file=itable)
    itable.close()