Determines which core genes from a pre-calculated database are present in each strain
10534
def cdsthreads(self):
        
        # Create and start threads
        for i in range(self.cpus):
            # Send the threads to the appropriate destination function
            threads = Thread(target=self.cds, args=())
            # Set the daemon to true - something to do with thread management
            threads.setDaemon(True)
            # Start the threading
            threads.start()
        for sample in self.metadata.samples:
            #
            sample[self.analysistype].corepresence = dict()
            self.cdsqueue.put(sample)
        self.cdsqueue.join()
139097
def total_core(self):
        
        corefile = os.path.join(self.reffilepath, self.analysistype, 'Escherichia', 'core_combined.fasta')
        for record in SeqIO.parse(corefile, 'fasta'):
            gene_name = record.id.split('-')[0]
            if gene_name not in self.coregenomes:
                self.coregenomes.append(gene_name)
1251120
def reporter(self):
        
        # Create a set of all the gene names without alleles or accessions e.g. sul1_18_AY260546 becomes sul1
        genedict = dict()
        # Load the notes file to a dictionary
        notefile = os.path.join(self.targetpath, 'notes.txt')
        with open(notefile, 'r') as notes:
            for line in notes:
                # Ignore comment lines - they will break the parsing
                if line.startswith('#'):
                    continue
                # Split the line on colons e.g. stx1Aa:  Shiga toxin 1, subunit A, variant a: has three variables after
                # the split: gene(stx1Aa), description(Shiga toxin 1, subunit A, variant a), and _(\n)
                try:
                    gene, description, _ = line.split(':')
                # There are exceptions to the parsing. Some lines only have one :, while others have three. Allow for
                # these possibilities.
                except ValueError:
                    try:
                        gene, description = line.split(':')
                    except ValueError:
                        gene, description, _, _ = line.split(':')
                # Set up the description dictionary
                genedict[gene] = description.replace(', ', '_').strip()
        # Find unique gene names with the highest percent identity
        for sample in self.runmetadata.samples:
            try:
                if sample[self.analysistype].results:
                    # Initialise a dictionary to store the unique genes, and their percent identities
                    sample[self.analysistype].uniquegenes = dict()
                    for name, identity in sample[self.analysistype].results.items():
                        # Split the name of the gene from the string e.g. stx1:11:Z36899:11 yields stx1
                        if ':' in name:
                            sample[self.analysistype].delimiter = ':'
                        else:
                            sample[self.analysistype].delimiter = '_'
                        genename = name.split(sample[self.analysistype].delimiter)[0]
                        # Set the best observed percent identity for each unique gene
                        try:
                            # Pull the previous best identity from the dictionary
                            bestidentity = sample[self.analysistype].uniquegenes[genename]
                            # If the current identity is better than the old identity, save it
                            if float(identity) > float(bestidentity):
                                sample[self.analysistype].uniquegenes[genename] = float(identity)
                        # Initialise the dictionary if necessary
                        except KeyError:
                            sample[self.analysistype].uniquegenes[genename] = float(identity)
            except AttributeError:
                raise
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        # Initialise strings to store the results
        data = 'Strain,Gene,Subtype/Allele,Description,Accession,PercentIdentity,FoldCoverage\n'
        with open(os.path.join(self.reportpath, self.analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                try:
                    if sample[self.analysistype].results:
                        # If there are many results for a sample, don't write the sample name in each line of the report
                        for name, identity in sorted(sample[self.analysistype].results.items()):
                            # Check to see which delimiter is used to separate the gene name, allele, accession, and
                            # subtype information in the header
                            if len(name.split(sample[self.analysistype].delimiter)) == 4:
                                # Split the name on the delimiter: stx2A:63:AF500190:d; gene: stx2A, allele: 63,
                                # accession: AF500190, subtype: d
                                genename, allele, accession, subtype = name.split(sample[self.analysistype].delimiter)
                            elif len(name.split(sample[self.analysistype].delimiter)) == 3:
                                # Treat samples without a subtype e.g. icaC:intercellular adhesion protein C: differently.
                                # Extract the allele as the 'subtype', and the gene name, and accession as above
                                genename, subtype, accession = name.split(sample[self.analysistype].delimiter)
                            else:
                                genename = name
                                subtype = ''
                                accession = ''
                            # Retrieve the best identity for each gene
                            percentid = sample[self.analysistype].uniquegenes[genename]
                            # If the percent identity of the current gene matches the best percent identity, add it to
                            # the report - there can be multiple occurrences of genes e.g.
                            # sul1,1,AY224185,100.00,840 and sul1,2,CP002151,100.00,927 are both included because they
                            # have the same 100% percent identity
                            if float(identity) == percentid:
                                # Treat the initial vs subsequent results for each sample slightly differently - instead
                                # of including the sample name, use an empty cell instead
                                try:
                                    description = genedict[genename]
                                except KeyError:
                                    description = 'na'
                                # Populate the results
                                data += '{samplename},{gene},{subtype},{description},{accession},{identity},{depth}\n'\
                                    .format(samplename=sample.name,
                                            gene=genename,
                                            subtype=subtype,
                                            description=description,
                                            accession=accession,
                                            identity=identity,
                                            depth=sample[self.analysistype].avgdepth[name])
                    else:
                        data += sample.name + '\n'
                except (KeyError, AttributeError):
                    data += sample.name + '\n'
            # Write the strings to the file
            report.write(data)
786361
def get_proteome_counts_impute_missing(prots_filtered_feathers, outpath, length_filter_pid=None,
                                       copynum_scale=False, copynum_df=None,
                                       force_rerun=False):
    """"""
    if ssbio.utils.force_rerun(flag=force_rerun, outfile=outpath):
        big_strain_counts_df = pd.DataFrame()
        first = True
        for feather in prots_filtered_feathers:
            loaded = load_feather(protein_feather=feather, length_filter_pid=length_filter_pid,
                                  copynum_scale=copynum_scale,
                                  copynum_df=copynum_df)
            if first:
                big_strain_counts_df = pd.DataFrame(index=_all_counts, columns=loaded.columns)
                first = False

            new_columns = list(set(loaded.columns.tolist()).difference(big_strain_counts_df.columns))
            if new_columns:
                for col in new_columns:
                    big_strain_counts_df[col] = big_strain_counts_df.mean(axis=1)

            not_in_loaded = list(set(big_strain_counts_df.columns).difference(loaded.columns.tolist()))
            if not_in_loaded:
                for col in not_in_loaded:
                    big_strain_counts_df[col] = big_strain_counts_df[col] + loaded.mean(axis=1)

            big_strain_counts_df = big_strain_counts_df.add(loaded, fill_value=0)

        if len(big_strain_counts_df) > 0:
            big_strain_counts_df.astype(float).reset_index().to_feather(outpath)
        return big_strain_counts_df
    else:
        return pd.read_feather(outpath).set_index('index')
1203080
def corewriter(self):
        
        printtime('Creating core allele files', self.start)
        for gene in sorted(self.genesequence):
            self.geneset.add(gene)
            # Set the name of the allele file
            genefile = os.path.join(self.coregenelocation, '{}.fasta'.format(gene))
            # If the file doesn't exist, create it
            if not os.path.isfile(genefile):
                with open(genefile, 'w') as core:
                    for count, sequence in enumerate(self.genesequence[gene]):
                        # The definition line is the gene name, and the allele number (count (+ 1 to compensate for
                        # base zero))
                        definitionline = '{}-{}'.format(gene, count + 1)
                        # Create a sequence record using BioPython
                        fasta = SeqRecord(Seq(sequence),
                                          # Without this, the header will be improperly formatted
                                          description='',
                                          # Use >:definitionline as the header
                                          id=definitionline)
                        # Use the SeqIO module to properly format the new sequence record
                        SeqIO.write(fasta, core, 'fasta')
                        for strain in self.coresequence[sequence]:
                            # Record the strain name, the gene name, and the allele number.
                            # [:-6] removes the contig number: 2014-SEQ-0276_00001 becomes 2014-SEQ-0276
                            try:
                                self.corealleles[strain[:-6]].update({gene: count + 1})
                            except KeyError:
                                self.corealleles[strain[:-6]] = {gene: count + 1}
            else:
                # If the file exists, don't recreate it; only iterate through the dictionary of gene sequences
                for count, sequence in enumerate(self.genesequence[gene]):
                    for strain in self.coresequence[sequence]:
                        # Populate the dictionary as above
                        try:
                            self.corealleles[strain[:-6]].update({gene: count + 1})
                        except KeyError:
                            self.corealleles[strain[:-6]] = {gene: count + 1}
        # Create a combined file of all the core genes to be used in typing strain(s) of interest
        if not os.path.isfile(os.path.join(self.coregenelocation, 'core_combined.fasta')):
            fastafiles = glob(os.path.join(self.coregenelocation, '*.fasta'))
            # Run the method for each allele
            self.combinealleles(fastafiles)
        # Run the profiler
        self.profiler()
1026293
def load_sequences_to_strains(self, joblib=False, cores=1, force_rerun=False):
        """"""
        log.info('Loading sequences to strain GEM-PROs...')
        if joblib:
            result = DictList(Parallel(n_jobs=cores)(delayed(self._load_sequences_to_strain)(s, force_rerun) for s in self.strain_ids))
        else:
            result = []
            for s in tqdm(self.strain_ids):
                result.append(self._load_sequences_to_strain(s, force_rerun))

        for strain_id, gp_seqs_path in result:
            self.strain_infodict[strain_id]['gp_seqs_path'] = gp_seqs_path
1298531
def _process_prb_strain_genotype_view(self, limit=None):
        

        line_counter = 0
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        LOG.info("Getting genotypes for strains")
        raw = '/'.join((self.rawdir, 'prb_strain_genotype_view'))
        with open(raw, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for line in filereader:
                line_counter += 1
                if line_counter == 1:
                    continue

                (strain_key, genotype_key) = line

                if self.test_mode is True:
                    if int(genotype_key) not in self.test_keys.get('genotype') \
                            and int(strain_key) not in self.test_keys.get('strain'):
                        continue

                strain_id = self.idhash['strain'].get(strain_key)
                if strain_id is None:
                    strain_id = self._makeInternalIdentifier(
                        'strain', strain_key)
                genotype_id = self.idhash['genotype'].get(genotype_key)
                if genotype_id is None:
                    genotype_id = self._makeInternalIdentifier(
                        'genotype', genotype_key)

                if strain_id is not None and genotype_id is not None:
                    self.strain_to_genotype_map[strain_id] = genotype_id

                graph.addTriple(strain_id, self.globaltt['has_genotype'], genotype_id)
                # TODO
                # verify if this should be contingent on the exactness or not
                # if qualifier == 'Exact':
                #     gu.addTriple(
                #       graph, strain_id,
                #       self.globaltt['has_genotype'],
                #       genotype_id)
                # else:
                #     gu.addXref(graph, strain_id, genotype_id)

                if not self.test_mode and limit is not None and line_counter > limit:
                    break

        return
989011
def gene_by_id(self, gene_id):
        
        if gene_id not in self._genes:
            field_names = [
                "seqname",
                "start",
                "end",
                "strand",
            ]
            optional_field_names = [
                "gene_name",
                "gene_biotype",
            ]
            # Do not look for gene_name and gene_biotype if they are
            # not in the database.
            field_names.extend([
                name for name in optional_field_names
                if self.db.column_exists("gene", name)
            ])
            result = self.db.query_one(
                field_names,
                filter_column="gene_id",
                filter_value=gene_id,
                feature="gene")
            if not result:
                raise ValueError("Gene not found: %s" % (gene_id,))

            gene_name, gene_biotype = None, None
            assert len(result) >= 4 and len(result) <= 6, \
                "Result is not the expected length: %d" % len(result)
            contig, start, end, strand = result[:4]
            if len(result) == 5:
                if "gene_name" in field_names:
                    gene_name = result[4]
                else:
                    gene_biotype = result[4]
            elif len(result) == 6:
                gene_name, gene_biotype = result[4:]

            self._genes[gene_id] = Gene(
                gene_id=gene_id,
                gene_name=gene_name,
                contig=contig,
                start=start,
                end=end,
                strand=strand,
                biotype=gene_biotype,
                genome=self)

        return self._genes[gene_id]
1018496
def preproc_data(data, gene_subset=False, **kwargs):
    
    import uncurl
    from uncurl.preprocessing import log1p, cell_normalize
    from sklearn.decomposition import TruncatedSVD
    data_subset = data
    if gene_subset:
        gene_subset = uncurl.max_variance_genes(data)
        data_subset = data[gene_subset, :]
    tsvd = TruncatedSVD(min(8, data_subset.shape[0] - 1))
    data_tsvd = tsvd.fit_transform(log1p(cell_normalize(data_subset)).T)
    return data_tsvd
1156988
def reporter(self):
        
        # Initialise variables
        header = ''
        row = ''
        databasedict = dict()
        # Load the database sequence type into a dictionary
        strainprofile = os.path.join(self.profilelocation, 'strainprofiles.txt')
        databaseprofile = DictReader(open(strainprofile))
        # Put the strain profile dictionary into a more easily searchable format
        for data in databaseprofile:
            databasedict[data['Strain']] = data['SequenceType']
        for sample in self.metadata.samples:
            closestmatches = list()
            if sample[self.analysistype].reportdir != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Populate the header with the appropriate data, including all the genes in the list of targets
                    header = 'Strain,SequenceType,Matches,Mismatches,NA,TotalGenes,ClosestDatabaseMatch,{},\n' \
                        .format(','.join(sorted(sample[self.analysistype].allelenames)))
                sortedmatches = sorted(sample[self.analysistype].profilematches.items(),
                                       key=operator.itemgetter(1), reverse=True)[0]
                closestseqtype = sortedmatches[0]
                # Pull out the closest database match
                for strain, seqtype in databasedict.items():
                    if seqtype == closestseqtype:
                        closestmatches.append(strain)
                sample[self.analysistype].closestseqtype = closestseqtype
                nummatches = int(sortedmatches[1])
                numna = 0
                queryallele = list()
                # Get all the alleles into a list
                for gene, allele in sorted(sample[self.analysistype].profiledata[closestseqtype].items()):
                    try:
                        # Extract the allele (anything after the -) from the allele matches
                        query = sample[self.analysistype].allelematches[gene].split('-')[1]
                        if allele == query:
                            queryallele.append(query)
                        else:
                            queryallele.append('{} ({})'.format(query, allele))
                    except KeyError:
                        queryallele.append('NA')
                        numna += 1
                mismatches = len(sample[self.analysistype].alleles) - nummatches - numna
                row += '{},{},{},{},{},{},{},{}'\
                    .format(sample.name, closestseqtype, nummatches, mismatches, numna,
                            len(sample[self.analysistype].alleles), ','.join(closestmatches), ','.join(queryallele))
                row += '\n'
        # Create the report folder
        make_path(self.reportpath)
        # Create the report containing all the data from all samples
        with open(os.path.join(self.reportpath, '{}.csv'.format(self.analysistype)), 'w') as combinedreport:
            # Write the results to this report
            combinedreport.write(header)
            combinedreport.write(row)
405982
def check_genes(self):
        
        gene_dna_set = set([g.dna for g in self.genes])
        assert gene_dna_set == self.dna_choices_set
1091344
def get_gene_names(self):
        """"""
        # Collect all gene names in network
        gene_names = []
        for node in self._nodes:
            members = node['data'].get('members')
            if members:
                gene_names += list(members.keys())
            else:
                if node['data']['name'].startswith('Group'):
                    continue
                gene_names.append(node['data']['name'])
        self._gene_names = gene_names
218926
def get_gene_seqs(database_path, gene):
    
    gene_path = database_path + "/" + gene + ".fsa"
    gene_seq = ""
    # Open fasta file
    with open(gene_path) as gene_file:
        header = gene_file.readline()
        for line in gene_file:
            seq = line.strip()
            gene_seq += seq
    return gene_seq
833600
def database(self):
        
        import sqlite3
        try:
            os.remove('{}/metadatabase.sqlite'.format(self.reportpath))
        except OSError:
            pass
        # Set the name of the database
        db = sqlite3.connect('{}/metadatabase.sqlite'.format(self.reportpath))
        # Create a cursor to allow access to the database
        cursor = db.cursor()
        # Set up the db
        cursor.execute('''
          CREATE TABLE IF NOT EXISTS Samples (
            id     INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT UNIQUE,
            name   TEXT UNIQUE
          )
        ''')
        # Create a variable to store the names of the header values for each individual table
        # This will store a set of all the headers from all the strains, as there can be some variability present, as
        # not all analyses are available for all taxonomic groups
        columns = dict()
        for sample in self.metadata:
            # Create a metadata object to store the new tables
            data = MetadataObject()
            data.name = sample.name
            # Insert each strain name into the Samples table
            cursor.execute('''
              INSERT OR IGNORE INTO Samples (name)
              VALUES ( ? )
            ''', (sample.name, ))
            # Each header in the .json file represents a major category e.g. ARMI, GeneSeekr, commands, etc. and
            # will be made into a separate table
            for header in sample.datastore.items():
                # Allow for certain analyses, such as core genome, not being performed on all strains
                try:

                    # Key and value: data description and data value e.g. targets present: 1012, etc.
                    for key, value in sorted(header[1].datastore.items()):
                        # Only the values consisting of dictionaries are of interest
                        if type(value) == dict:
                            # Clean the column names so there are no issues entering names into the database
                            cleanedcolumn = self.columnclean(key)
                            # Set the table name
                            tablename = '{}_{}'.format(header[0].replace('.', '_'), cleanedcolumn)
                            # Create the table (if it doesn't already exist)
                            cursor.execute('''
                                          CREATE TABLE IF NOT EXISTS {} (
                                            sample_id INTEGER
                                          )
                                          '''.format(tablename))
                            # Add the attributes with the dictionaries (values) to the metadata object
                            setattr(data, tablename, GenObject(value))
                            for gene, result in sorted(value.items()):
                                # Add the data header to the dictionary
                                try:
                                    columns[tablename].add(gene)
                                # Initialise the dictionary the first time a table name is encountered
                                except KeyError:
                                    columns[tablename] = set()
                                    columns[tablename].add(str(gene))
                except (AttributeError, IndexError):
                    pass
            self.tabledata.append(data)
        # Iterate through the dictionary containing all the data headers
        for table, setofheaders in sorted(columns.items()):
            # Each header will be used as a column in the appropriate table
            for cleanedcolumn in sorted(setofheaders):
                # Alter the table by adding each header as a column
                cursor.execute('''
                  ALTER TABLE {}
                  ADD COLUMN {} TEXT
                '''.format(table, cleanedcolumn))
            # Iterate through the samples and pull out the data for each table/column
            # for sample in self.metadata:
            for sample in self.tabledata:
                # Find the id associated with each sample in the Sample table
                cursor.execute('''
                  SELECT id from Samples WHERE name=?
                ''', (sample.name,))
                sampleid = cursor.fetchone()[0]
                # Add the sample_id to the table
                cursor.execute('''
                  INSERT OR IGNORE INTO {}
                  (sample_id) VALUES ("{}")
                  '''.format(table, sampleid))
                # Add the data to the table
                try:
                    # Find the data for each table/column
                    for item in sorted(sample[table].datastore.items()):
                        # Clean the names
                        cleanedcolumn = self.columnclean(str(item[0]))
                        # Add the data to the column of the appropriate table,
                        # where the sample_id matches the current strain
                        cursor.execute('''
                          UPDATE {}
                          SET {} = ?
                          WHERE sample_id = {}
                          '''.format(table, cleanedcolumn, sampleid), (str(item[1]), ))
                except KeyError:
                    pass
        # Commit the changes to the database
        db.commit()
1053895
def dalignbedannot2daggbyquery(cfg):
    
    datatmpd=cfg['datatmpd']

    dalignbedannot=del_Unnamed(pd.read_csv(cfg['dalignbedannotp'],sep='\t',low_memory=False))
    
    daggbyqueryp=f'{datatmpd}/10_daggbyquery.tsv'      
    logging.info(basename(daggbyqueryp))
    if not exists(daggbyqueryp) or cfg['force']:
        dalignbedannot=dfliststr2dflist(dalignbedannot,
                                ['type', 'gene name', 'gene id', 'transcript id', 'protein id', 'exon id'],
                                colfmt='tuple')    
        dalignbedannot['alternate alignments count']=1
        import itertools
        from rohan.dandage.io_sets import unique
        def unique_dropna(l): return unique(l,drop='nan')
        def merge_unique_dropna(l): return unique(list(itertools.chain(*l)),drop='nan')
        cols2aggf={'id':unique_dropna,
         'type':merge_unique_dropna,
         'gene name':merge_unique_dropna,
         'gene id':merge_unique_dropna,
         'transcript id':merge_unique_dropna,
         'protein id':merge_unique_dropna,
         'exon id':merge_unique_dropna,
         'region':unique_dropna,
         'alternate alignments count':sum,
        }
        daggbyquery=dalignbedannot.groupby('query id').agg(cols2aggf)
        daggbyquery.to_csv(daggbyqueryp,sep='\t')
        daggbyquery.to_csv(cfg['dalignannotedp'],sep='\t')
    return cfg
191635
def find_databases(databases):
    
    # 16 ribosomal proteins in their expected order
    proteins = ['L15', 'L18', 'L6', 'S8', 'L5', 'L24', 'L14',
            'S17', 'L16', 'S3', 'L22', 'S19', 'L2', 'L4', 'L3', 'S10']
    # curated databases
    protein_databases = {
                'L14': 'rpL14_JGI_MDM.filtered.faa',
                'L15': 'rpL15_JGI_MDM.filtered.faa',
                'L16': 'rpL16_JGI_MDM.filtered.faa',
                'L18': 'rpL18_JGI_MDM.filtered.faa',
                'L22': 'rpL22_JGI_MDM.filtered.faa',
                'L24': 'rpL24_JGI_MDM.filtered.faa',
                'L2': 'rpL2_JGI_MDM.filtered.faa',
                'L3': 'rpL3_JGI_MDM.filtered.faa',
                'L4': 'rpL4_JGI_MDM.filtered.faa',
                'L5': 'rpL5_JGI_MDM.filtered.faa',
                'L6': 'rpL6_JGI_MDM.filtered.faa',
                'S10': 'rpS10_JGI_MDM.filtered.faa',
                'S17': 'rpS17_JGI_MDM.filtered.faa',
                'S19': 'rpS19_JGI_MDM.filtered.faa',
                'S3': 'rpS3_JGI_MDM.filtered.faa',
                'S8': 'rpS8_JGI_MDM.filtered.faa'}
    protein_databases = {key: '%s/%s' % (databases, database) \
            for key, database in list(protein_databases.items())}
    return proteins, protein_databases
408896
def parse(self, limit=None):
        
        if limit is not None:
            LOG.info("Only parsing first %d rows", limit)

        sgd_file = '/'.join((self.rawdir, self.files['sgd_phenotype']['file']))
        columns = [
            'Feature Name', 'Feature Type', 'Gene Name', 'SGDID', 'Reference',
            'Experiment Type', 'Mutant Type', 'Allele', 'Strain Background',
            'Phenotype', 'Chemical', 'Condition', 'Details', 'Reporter']
        sgd_df = pd.read_csv(sgd_file, sep='\t', names=columns)
        records = sgd_df.to_dict(orient='records')
        for index, assoc in enumerate(records):
            if isinstance(assoc['Gene Name'], str):
                if limit is not None and index > limit:
                    break
                self.make_association(assoc)

        return