Runs through a protein table and adds quant data from ANOTHER protein
    table that contains that data.
170108
def add_isoquant_data(proteins, quantproteins, quantacc, quantfields):
    """"""
    for protein in base_add_isoquant_data(proteins, quantproteins,
                                          prottabledata.HEADER_PROTEIN,
                                          quantacc, quantfields):
        yield protein
139906
def add_all_exchange_reactions(model, compartment, allow_duplicates=False):
    

    all_reactions = {}
    if not allow_duplicates:
        # TODO: Avoid adding reactions that already exist in the database.
        # This should be integrated in the database.
        for rxnid in model.database.reactions:
            rx = model.database.get_reaction(rxnid)
            all_reactions[rx] = rxnid

    added = set()
    added_compounds = set()
    initial_compounds = set(model.compounds)
    reactions = set(model.database.reactions)
    for model_compound in initial_compounds:
        compound = model_compound.in_compartment(compartment)
        if compound in added_compounds:
            continue

        rxnid_ex = create_exchange_id(reactions, compound)

        reaction_ex = Reaction(Direction.Both, {compound: -1})
        if reaction_ex not in all_reactions:
            model.database.set_reaction(rxnid_ex, reaction_ex)
            reactions.add(rxnid_ex)
        else:
            rxnid_ex = all_reactions[reaction_ex]

        if not model.has_reaction(rxnid_ex):
            added.add(rxnid_ex)
        model.add_reaction(rxnid_ex)
        added_compounds.add(compound)

    return added
448697
def addRnaQuantificationSet(self):
        
        self._openRepo()
        dataset = self._repo.getDatasetByName(self._args.datasetName)
        if self._args.name is None:
            name = getNameFromPath(self._args.filePath)
        else:
            name = self._args.name
        rnaQuantificationSet = rna_quantification.SqliteRnaQuantificationSet(
            dataset, name)
        referenceSetName = self._args.referenceSetName
        if referenceSetName is None:
            raise exceptions.RepoManagerException(
                "A reference set name must be provided")
        referenceSet = self._repo.getReferenceSetByName(referenceSetName)
        rnaQuantificationSet.setReferenceSet(referenceSet)
        rnaQuantificationSet.populateFromFile(self._args.filePath)
        rnaQuantificationSet.setAttributes(json.loads(self._args.attributes))
        self._updateRepo(
            self._repo.insertRnaQuantificationSet, rnaQuantificationSet)
492294
def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'):
    """"""
    if not as_string:
        startmessage = '\nStarting analysis of %s\n' % pdbfile.split('/')[-1]
    else:
        startmessage = "Starting analysis from stdin.\n"
    write_message(startmessage)
    write_message('='*len(startmessage)+'\n')
    mol = PDBComplex()
    mol.output_path = outpath
    mol.load_pdb(pdbfile, as_string=as_string)
    # #@todo Offers possibility for filter function from command line (by ligand chain, position, hetid)
    for ligand in mol.ligands:
        mol.characterize_complex(ligand)

    create_folder_if_not_exists(outpath)

    # Generate the report files
    streport = StructureReport(mol, outputprefix=outputprefix)

    config.MAXTHREADS = min(config.MAXTHREADS, len(mol.interaction_sets))

    ######################################
    # PyMOL Visualization (parallelized) #
    ######################################

    if config.PYMOL or config.PICS:
        try:
            from plip.modules.visualize import visualize_in_pymol
        except ImportError:
            from modules.visualize import visualize_in_pymol
        complexes = [VisualizerData(mol, site) for site in sorted(mol.interaction_sets)
                     if not len(mol.interaction_sets[site].interacting_res) == 0]
        if config.MAXTHREADS > 1:
            write_message('\nGenerating visualizations in parallel on %i cores ...' % config.MAXTHREADS)
            parfn = parallel_fn(visualize_in_pymol)
            parfn(complexes, processes=config.MAXTHREADS)
        else:
            [visualize_in_pymol(plcomplex) for plcomplex in complexes]

    if config.XML:  # Generate report in xml format
        streport.write_xml(as_string=config.STDOUT)

    if config.TXT:  # Generate report in txt (rst) format
        streport.write_txt(as_string=config.STDOUT)
511859
def run(bam_file, data, out_dir):
    """"""
    m = {"base": None, "secondary": []}
    m.update(_mirbase_stats(data, out_dir))
    m["secondary"].append(_seqcluster_stats(data, out_dir))
531185
def updateQTable( self, reward, current_state ):
        """"""
        
        # retrieve the old action value from the Q-table (indexed by the previous state and the previous action)
        old_q = self.q_table[self.prev_s][self.prev_a]
        
        # TODO: what should the new action value be?
        new_q = old_q
        
        # assign the new action value to the Q-table
        self.q_table[self.prev_s][self.prev_a] = new_q
916640
def add_ms1_quant_from_top3_mzidtsv(proteins, psms, headerfields, protcol):
    """"""
    if not protcol:
        protcol = mzidtsvdata.HEADER_MASTER_PROT
    top_ms1_psms = generate_top_psms(psms, protcol)
    for protein in proteins:
        prot_acc = protein[prottabledata.HEADER_PROTEIN]
        prec_area = calculate_protein_precursor_quant(top_ms1_psms, prot_acc)
        outprotein = {k: v for k, v in protein.items()}
        outprotein[headerfields['precursorquant'][
            prottabledata.HEADER_AREA][None]] = str(prec_area)
        yield outprotein
1085323
def add_all_database_reactions(model, compartments):
    

    added = set()
    for rxnid in model.database.reactions:
        reaction = model.database.get_reaction(rxnid)
        if all(compound.compartment in compartments
               for compound, _ in reaction.compounds):
            if not model.has_reaction(rxnid):
                added.add(rxnid)
            model.add_reaction(rxnid)

    return added
107285
def run_rnaseq(bam_file, data, out_dir):
    
    strandedness = {"firststrand": "strand-specific-reverse",
                    "secondstrand": "strand-specific-forward",
                    "unstranded": "non-strand-specific"}

    # Qualimap results should be saved to a directory named after sample.
    # MultiQC (for parsing additional data) picks the sample name after the dir as follows:
    #   <sample name>/raw_data_qualimapReport/insert_size_histogram.txt
    results_dir = os.path.join(out_dir, dd.get_sample_name(data))
    results_file = os.path.join(results_dir, "rnaseq_qc_results.txt")
    report_file = os.path.join(results_dir, "qualimapReport.html")
    config = data["config"]
    gtf_file = dd.get_gtf_file(data)
    library = strandedness[dd.get_strandedness(data)]
    if not utils.file_exists(results_file):
        with file_transaction(data, results_dir) as tx_results_dir:
            utils.safe_makedir(tx_results_dir)
            bam.index(bam_file, config)
            cmd = _rnaseq_qualimap_cmd(data, bam_file, tx_results_dir, gtf_file, library)
            do.run(cmd, "Qualimap for {}".format(dd.get_sample_name(data)))
            tx_results_file = os.path.join(tx_results_dir, "rnaseq_qc_results.txt")
            cmd = "sed -i 's/bam file = .*/bam file = %s.bam/' %s" % (dd.get_sample_name(data), tx_results_file)
            do.run(cmd, "Fix Name Qualimap for {}".format(dd.get_sample_name(data)))
    metrics = _parse_rnaseq_qualimap_metrics(report_file)
    metrics.update(_detect_duplicates(bam_file, results_dir, data))
    metrics.update(_detect_rRNA(data, results_dir))
    metrics.update({"Average_insert_size": salmon.estimate_fragment_size(data)})
    metrics = _parse_metrics(metrics)
    # Qualimap output folder (results_dir) needs to be named after the sample (see comments above). However, in order
    # to keep its name after upload, we need to put  the base QC file (results_file) into the root directory (out_dir):
    base_results_file = os.path.join(out_dir, os.path.basename(results_file))
    shutil.copyfile(results_file, base_results_file)
    return {"base": base_results_file,
            "secondary": _find_qualimap_secondary_files(results_dir, base_results_file),
            "metrics": metrics}
835887
def merge(self, other_roc):
        
        if other_roc.thresholds.size == self.thresholds.size and np.all(other_roc.thresholds == self.thresholds):
            self.contingency_tables += other_roc.contingency_tables
        else:
            print("Input table thresholds do not match.")
153269
def get_peptide_quant(quantdata, quanttype):
    """"""
    parsefnx = {'precur': max}
    quantfloats = []
    for q in quantdata:
        try:
            quantfloats.append(float(q))
        except(TypeError, ValueError):
            pass
    if not quantfloats:
        return 'NA'
    return str(parsefnx[quanttype](quantfloats))
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:]}
1261722
def create_tables(self, tables):
        """"""
        cursor = self.get_cursor()
        for table in tables:
            columns = mslookup_tables[table]
            try:
                cursor.execute('CREATE TABLE {0}({1})'.format(
                    table, ', '.join(columns)))
            except sqlite3.OperationalError as error:
                print(error)
                print('Warning: Table {} already exists in database, will '
                      'add to existing tables instead of creating '
                      'new.'.format(table))
            else:
                self.conn.commit()
1294689
def new_table(self, name, add_id=True, **kwargs):
        '''
        '''
        from . import Table
        from .exc import NotFoundError

        try:
            table = self.table(name)
            extant = True
        except NotFoundError:

            extant = False

            if 'sequence_id' not in kwargs:
                kwargs['sequence_id'] = self._database.next_sequence_id(Dataset, self.vid, Table)

            table = Table(name=name, d_vid=self.vid, **kwargs)

            table.update_id()

        # Update possibly extant data
        table.data = dict(
            (list(table.data.items()) if table.data else []) + list(kwargs.get('data', {}).items()))

        for key, value in list(kwargs.items()):

            if not key:
                continue
            if key[0] != '_' and key not in ['vid', 'id', 'id_', 'd_id', 'name', 'sequence_id', 'table', 'column', 'data']:
                setattr(table, key, value)

        if add_id:
            table.add_id_column()

        if not extant:
            self.tables.append(table)

        return table
2808
def quantitate(data):
    
    data = to_single_data(to_single_data(data))
    data = generate_transcript_counts(data)[0][0]
    data["quant"] = {}
    if "sailfish" in dd.get_expression_caller(data):
        data = to_single_data(sailfish.run_sailfish(data)[0])
        data["quant"]["tsv"] = data["sailfish"]
        data["quant"]["hdf5"] = os.path.join(os.path.dirname(data["sailfish"]), "abundance.h5")
    if ("kallisto" in dd.get_expression_caller(data) or "pizzly" in dd.get_fusion_caller(data, [])):
        data = to_single_data(kallisto.run_kallisto_rnaseq(data)[0])
        data["quant"]["tsv"] = os.path.join(data["kallisto_quant"], "abundance.tsv")
        data["quant"]["hdf5"] = os.path.join(data["kallisto_quant"], "abundance.h5")
    if (os.path.exists(os.path.join(data["kallisto_quant"], "fusion.txt"))):
        data["quant"]["fusion"] = os.path.join(data["kallisto_quant"], "fusion.txt")
    else:
        data["quant"]["fusion"] = None
    if "salmon" in dd.get_expression_caller(data):
        data = to_single_data(salmon.run_salmon_reads(data)[0])
        data["quant"]["tsv"] = data["salmon"]
        data["quant"]["hdf5"] = os.path.join(os.path.dirname(data["salmon"]), "abundance.h5")
    return [[data]]
918776
def run(self, data):
        
        if self.normalize_data:
            data = cell_normalize(data)
        M, W, ll = poisson_estimate_state(data, **self.params)
        outputs = []
        if self.return_w:
            outputs.append(W)
        if self.return_m:
            outputs.append(M)
        if self.return_mw:
            outputs.append(M.dot(W))
        if self.return_mds:
            X = dim_reduce(M, W, 2)
            outputs.append(X.T.dot(W))
        return outputs, ll
62596
def alignprotein(self, sample, target):
        
        # Initialise dictionaries
        sample[self.analysistype].dnaseq = dict()
        sample[self.analysistype].protseq = dict()
        sample[self.analysistype].ntindex = dict()
        sample[self.analysistype].aaindex = dict()
        sample[self.analysistype].ntalign = dict()
        sample[self.analysistype].aaalign = dict()
        sample[self.analysistype].aaidentity = dict()
        # Remove any gaps incorporated into the sequence
        sample[self.analysistype].targetsequence[target] = \
            sample[self.analysistype].targetsequence[target].replace('-', '')
        # In order to properly translate the nucleotide sequence, BioPython requests that the sequence is a multiple of
        # three - not partial codons. Trim the sequence accordingly
        remainder = 0 - len(sample[self.analysistype].targetsequence[target]) % 3
        seq = sample[self.analysistype].targetsequence[target] if remainder == 0 \
            else sample[self.analysistype].targetsequence[target][:remainder]
        # Set the DNA and protein sequences of the target in the sample
        sample[self.analysistype].dnaseq[target] = Seq(seq, IUPAC.unambiguous_dna)
        # Translate the nucleotide sequence
        sample[self.analysistype].protseq[target] = str(sample[self.analysistype].dnaseq[target].translate())
        for targetfile in self.targetfiles:
            # Trim the reference sequence to multiples of three
            refremainder = 0 - len(self.records[targetfile][target].seq) % 3
            refseq = str(self.records[targetfile][target].seq) if refremainder % 3 == 0 \
                else str(self.records[targetfile][target].seq)[:refremainder]
            # Translate the nucleotide sequence of the reference sequence
            refdna = Seq(refseq, IUPAC.unambiguous_dna)
            refprot = str(refdna.translate())
            # Use pairwise2 to perform a local alignment with the following parameters:
            # x     No match parameters. Identical characters have score of 1, otherwise 0.
            # s     Same open (-1)  and extend (-.1) gap penalties for both sequences
            ntalignments = pairwise2.align.localxs(seq, refseq, -1, -.1)
            # Use format_alignment to create a formatted alignment that is subsequently split on newlines e.g.
            '''
            ACCGT
            | ||
            A-CG-
            Score=3
            '''
            ntformat = (str(format_alignment(*ntalignments[0])).split('\n'))
            # Align the nucleotide sequence of the reference (ntalignments[2]) to the sample (ntalignments[0]).
            # If the corresponding bases match, add a |, otherwise a space
            ntalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(ntformat[0], ntformat[2])))
            # Create the nucleotide alignment: the sample sequence, the (mis)matches, and the reference sequence
            sample[self.analysistype].ntalign[target] = self.interleaveblastresults(ntformat[0], ntformat[2])
            # Regex to determine location of mismatches in the sequences
            count = 0
            sample[self.analysistype].ntindex[target] = str()
            for snp in re.finditer(' ', ntalignment):
                # If there are many SNPs, then insert line breaks for every 10 SNPs
                if count <= 10:
                    sample[self.analysistype].ntindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].ntindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1
            # Perform the same steps, except for the amino acid sequence
            aaalignments = pairwise2.align.localxs(sample[self.analysistype].protseq[target], refprot, -1, -.1)
            aaformat = (str(format_alignment(*aaalignments[0])).split('\n'))
            aaalignment = ''.join(map(lambda x: '|' if len(set(x)) == 1 else ' ',
                                      zip(aaformat[0], aaformat[2])))
            sample[self.analysistype].aaidentity[target] = '{:.2f}'\
                .format(float(aaalignment.count('|')) / float(len(aaalignment)) * 100)
            sample[self.analysistype].aaalign[target] = self.interleaveblastresults(aaformat[0], aaformat[2])
            count = 0
            sample[self.analysistype].aaindex[target] = str()
            for snp in re.finditer(' ', aaalignment):
                if count <= 10:
                    sample[self.analysistype].aaindex[target] += str(snp.start()) + ';'
                else:
                    sample[self.analysistype].aaindex[target] += '\n' + str(snp.start()) + ';'
                    count = 0
                count += 1