Returns BlastResult from input seqs, using blastp.

    Need to add doc string
509737
def blastp(seqs, blast_db="nr", e_value="1e-20", max_hits=200,
           working_dir=tempfile.gettempdir(), blast_mat_root=None,
           extra_params={}):
    

    # set up params to use with blastp
    params = {
        # matrix
        "-M":"BLOSUM62",

        # max procs
        "-a":"1",

        # expectation
        "-e":e_value,

        # max seqs to show
        "-b":max_hits,

        # max one line descriptions
        "-v":max_hits,

        # program
        "-p":"blastp"
    }
    params.update(extra_params)

    # blast
    blast_res =  blast_seqs(seqs,
        Blastall,
        blast_mat_root=blast_mat_root,
        blast_db=blast_db,
        params=params,
        add_seq_names=False,
        WorkingDir=working_dir
        )

    # get prot id map
    if blast_res['StdOut']:
        lines = [x for x in blast_res['StdOut']]
        return BlastResult(lines)

    return None
370664
def convert_block_dicts_to_string(self, block_1st2nd, block_1st, block_2nd, block_3rd):
        """"""
        out = ""
        # We need 1st and 2nd positions
        if self.codon_positions in ['ALL', '1st-2nd']:
            for gene_code, seqs in block_1st2nd.items():
                out += '>{0}_1st-2nd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        elif self.codon_positions == '1st':
            for gene_code, seqs in block_1st.items():
                out += '>{0}_1st\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        elif self.codon_positions == '2nd':
            for gene_code, seqs in block_2nd.items():
                out += '>{0}_2nd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq

        # We also need 3rd positions
        if self.codon_positions in ['ALL', '3rd']:
            for gene_code, seqs in block_3rd.items():
                out += '\n>{0}_3rd\n----\n'.format(gene_code)
                for seq in seqs:
                    out += seq
        return out
756761
def get_seqs_fasta(seqs, names, out_fa):
    """"""
    with open(out_fa, 'w') as fa_handle:
        for s, n in itertools.izip(seqs, names):
            print(">cx{1}-{0}\n{0}".format(s, n), file=fa_handle)
    return out_fa
888172
def blastparser(self, report, sample):
        
        # Open the sequence profile file as a dictionary
        blastdict = DictReader(open(report), fieldnames=self.fieldnames, dialect='excel-tab')
        resultdict = dict()
        # Initialise a dictionary to store all the target sequences
        sample[self.analysistype].targetsequence = dict()
        # Go through each BLAST result
        for row in blastdict:
            # Calculate the percent identity and extract the bitscore from the row
            # Percent identity is the (length of the alignment - number of mismatches) / total subject length
            percentidentity = float('{:0.2f}'.format((float(row['positives']) - float(row['gaps'])) /
                                                     float(row['subject_length']) * 100))
            target = row['subject_id']
            # If the percent identity is greater than the cutoff
            if percentidentity >= self.cutoff:
                # Update the dictionary with the target and percent identity
                resultdict.update({target: percentidentity})
                # Determine if the orientation of the sequence is reversed compared to the reference
                if int(row['subject_end']) < int(row['subject_start']):
                    # Create a sequence object using Biopython
                    seq = Seq(row['query_sequence'], IUPAC.unambiguous_dna)
                    # Calculate the reverse complement of the sequence
                    querysequence = str(seq.reverse_complement())
                # If the sequence is not reversed, use the sequence as it is in the output
                else:
                    querysequence = row['query_sequence']
                # Add the sequence in the correct orientation to the sample
                sample[self.analysistype].targetsequence[target] = querysequence
            # Add the percent identity to the object
            sample[self.analysistype].blastresults = resultdict
        # Populate missing results with 'NA' values
        if len(resultdict) == 0:
            sample[self.analysistype].blastresults = 'NA'
1113405
def blast(args):
    
    p = OptionParser(blast.__doc__)

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    btabfile, = args
    btab = Btab(btabfile)
    for b in btab:
        print(b.blastline)
961494
def findCodons(seq, codons):
    
    seqLen = len(seq)
    start = 0
    while start < seqLen:
        triplet = str(seq[start:start + 3])
        if triplet in codons:
            yield start
        start = start + 3
879079
function() {
            while (i < tb.seqCount) {
                var ns = ba[p];
                if (p + ns + 2 + (tb.version == 0 ? 4 : 8) >= ba.length) {
                    headerBlocksFetched += headerBlockSize;
                    headerBlockSize = Math.max(HEADER_BLOCK_SIZE,Math.floor(headerBlocksFetched*tb.seqCount/i));
                    return tb.data.slice(o + p, headerBlockSize).fetch(function (r) {
                        o += p;
                        p = 0;
                        ba = new Uint8Array(r);
                        parseSeqInfo();
                    });
                } else {
                    ++p;
                    var name = '';
                    for (var j = 1; j <= ns; ++j) {
                        name += String.fromCharCode(ba[p++]);
                    }
                    if (tb.version == 0) {
                        var offset = tb.readInt(ba, p);
                        p += 4;
                    } else {
                        var offset = tb.readInt64(ba, p);
                        p += 8;
                    }
                    tb.seqDict[name] = new TwoBitSeq(tb, offset);
                    ++i;
                }
            }
            return cnt(tb);
        }
655228
def runblast(self, assembly, allele, sample):
        
        genome = os.path.split(assembly)[1].split('.')[0]
        # Run the BioPython BLASTn module with the genome as query, fasta(target gene) as db.
        # Do not re-perform the BLAST search each time
        make_path(sample[self.analysistype].reportdir)
        try:
            report = glob('{}{}*rawresults*'.format(sample[self.analysistype].reportdir, genome))[0]
            size = os.path.getsize(report)
            if size == 0:
                os.remove(report)
                report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                          time.strftime("%Y.%m.%d.%H.%M.%S"))
        except IndexError:
            report = '{}{}_rawresults_{:}.csv'.format(sample[self.analysistype].reportdir, genome,
                                                      time.strftime("%Y.%m.%d.%H.%M.%S"))
        db = allele.split('.')[0]
        # BLAST command line call. Note the mildly restrictive evalue, and the high number of alignments.
        # Due to the fact that all the targets are combined into one database, this is to ensure that all potential
        # alignments are reported. Also note the custom outfmt: the doubled quotes are necessary to get it work
        blastn = NcbiblastnCommandline(query=assembly, db=db, evalue='1E-20', num_alignments=1000000,
                                       num_threads=12,
                                       outfmt="'6 qseqid sseqid positive mismatch gaps "
                                              "evalue bitscore slen length qstart qend qseq sstart send'",
                                       out=report)
        # Save the blast command in the metadata
        sample[self.analysistype].blastcommand = str(blastn)
        sample[self.analysistype].blastreport = report
        if not os.path.isfile(report):
            # Run BLAST
            blastn()
        # Run the blast parsing module
        self.blastparser(report, sample)
105746
def _generate_docstring_return_section(self, return_vals, header,
                                           return_element_name,
                                           return_element_type,
                                           placeholder, indent):
        """"""
        # If all return values are None, return none
        non_none_vals = [return_val for return_val in return_vals
                         if return_val and return_val != 'None']
        if not non_none_vals:
            return header + indent + 'None.'

        # Get only values with matching brackets that can be cleaned up
        non_none_vals = [return_val.strip(' ()\t\n').rstrip(',')
                         for return_val in non_none_vals]
        non_none_vals = [re.sub('([\"\'])(?:(?=(\\\\?))\\2.)*?\\1',
                                '"string"', return_val)
                         for return_val in non_none_vals]
        unambiguous_vals = []
        for return_val in non_none_vals:
            try:
                cleaned_val = self.find_top_level_bracket_locations(return_val)
            except IndexError:
                continue
            unambiguous_vals.append(cleaned_val)
        if not unambiguous_vals:
            return header + placeholder

        # If remaining are a mix of tuples and not, return single placeholder
        single_vals, tuple_vals = [], []
        for return_val in unambiguous_vals:
            (tuple_vals.append(return_val) if ',' in return_val
             else single_vals.append(return_val))
        if single_vals and tuple_vals:
            return header + placeholder

        # If return values are tuples of different length, return a placeholder
        if tuple_vals:
            num_elements = [return_val.count(',') + 1
                            for return_val in tuple_vals]
            if num_elements.count(num_elements[0]) != len(num_elements):
                return header + placeholder
            num_elements = num_elements[0]
        else:
            num_elements = 1

        # If all have the same len but some ambiguous return that placeholders
        if len(unambiguous_vals) != len(non_none_vals):
            return header + '\n'.join(
                [placeholder for __ in range(num_elements)])

        # Handle tuple (or single) values position by position
        return_vals_grouped = zip(*[
            [return_element.strip() for return_element in
             return_val.split(',')]
            for return_val in unambiguous_vals])
        return_elements_out = []
        for return_vals_group in return_vals_grouped:
            return_elements_out.append(
                self.parse_return_elements(return_vals_group,
                                           return_element_name,
                                           return_element_type,
                                           placeholder))

        return header + '\n'.join(return_elements_out)
53163
def make_arrays(blastfile, qpadbed, spadbed, qpadnames, spadnames):
    
    m, n = len(qpadnames), len(spadnames)
    qpadorder, spadorder = qpadbed.order, spadbed.order
    qpadid = dict((a, i) for i, a in enumerate(qpadnames))
    spadid = dict((a, i) for i, a in enumerate(spadnames))
    qpadlen = dict((a, len(b)) for a, b in qpadbed.sub_beds())
    spadlen = dict((a, len(b)) for a, b in spadbed.sub_beds())

    qsize, ssize = len(qpadbed), len(spadbed)

    assert sum(qpadlen.values()) == qsize
    assert sum(spadlen.values()) == ssize

    # Populate arrays of observed counts and expected counts
    logging.debug("Initialize array of size ({0} x {1})".format(m, n))
    observed = np.zeros((m, n))
    fp = open(blastfile)
    all_dots = 0
    for row in fp:
        b = BlastLine(row)
        qi, q = qpadorder[b.query]
        si, s = spadorder[b.subject]
        qseqid, sseqid = q.seqid, s.seqid
        qsi, ssi = qpadid[qseqid], spadid[sseqid]
        observed[qsi, ssi] += 1
        all_dots += 1

    assert int(round(observed.sum())) == all_dots

    logging.debug("Total area: {0} x {1}".format(qsize, ssize))
    S = qsize * ssize
    expected = np.zeros((m, n))
    qsum = 0
    for i, a in enumerate(qpadnames):
        alen = qpadlen[a]
        qsum += alen
        for j, b in enumerate(spadnames):
            blen = spadlen[b]
            expected[i, j] = all_dots * alen * blen * 1. / S

    assert int(round(expected.sum())) == all_dots

    # Calculate the statistical significance for each cell
    from scipy.stats.distributions import poisson
    M = m * n  # multiple testing
    logmp = np.zeros((m, n))
    for i in xrange(m):
        for j in xrange(n):
            obs, exp = observed[i, j], expected[i, j]
            pois = max(poisson.pmf(obs, exp), 1e-250)  # Underflow
            logmp[i, j] = max(- log(pois), 0)

    return logmp
1013855
def get_sequences(options, fasta_file_contents):
    ''' 
    '''
    errors = []
    fasta_files_str = ", ".join(fasta_file_contents.keys())
    fasta_records = None
    reverse_mapping = {}

    try:
        fasta_records, reverse_mapping = parse_FASTA_files(options, fasta_file_contents)
        if not fasta_records:
            errors.append("No protein sequences found in the FASTA file(s) %s." % fasta_files_str)
    except Exception, e:
        e = '\n'.join([l for l in traceback.format_exc(), str('e') if l.strip()])
        errors.append("Error parsing FASTA file(s) %s:\n%s" % (fasta_files_str, str(e)))

    if not fasta_records:
        return None, {}, errors

    colorprinter.message('Found %d protein sequence(s).' % len(fasta_records))
    return fasta_records, reverse_mapping, errors
785750
def _read_function(schema):
    
    def func(
        filename,
        seq_label='sequence',
        alphabet=None,
        use_uids=True,
        **kwargs):
        # Use generic write class to write data.
        return _read(
            filename=filename,
            schema=schema,
            seq_label=seq_label,
            alphabet=alphabet,
            use_uids=use_uids,
            **kwargs
        )
    # Update docs
    func.__doc__ = _read_doc_template(schema)
    return func
861033
def chain(args):
    
    p = OptionParser(chain.__doc__)
    p.add_option("--dist", dest="dist",
            default=100, type="int",
            help="extent of flanking regions to search [default: %default]")

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    blastfile, = args
    dist = opts.dist
    assert dist > 0

    blast = BlastSlow(blastfile)
    logging.debug("A total of {} records imported".format(len(blast)))
    chained_hsps = chain_HSPs(blast, xdist=dist, ydist=dist)
    logging.debug("A total of {} records after chaining".format(len(chained_hsps)))

    for b in chained_hsps:
        print(b)
890402
def get_doc(additional_doc=False,
            field_prefix='$',
            field_suffix=':',
            indent=4):
    
    if additional_doc:
        f = fields.copy()
        f.update(additional_doc)
    else:
        f = fields
    field_length = get_max_field_length(f)
    field_length = field_length + len(field_prefix) + len(field_suffix) + 4
    description_indent = ' ' * (indent + field_length)
    output = ''
    for field, description in sorted(f.items()):
        description = description['description']
        field = ' ' * indent + field_prefix + field + ':'
        output += field.ljust(field_length) + \
            textwrap.fill(
                description,
                width=78,
                initial_indent=description_indent,
                subsequent_indent=description_indent
            )[field_length:] + '\n\n\n'

    return output
1262700
def parseblast(self):
        
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA' and sample[self.analysistype].assemblyfile != 'NA':
                # Initialise variables
                sample[self.analysistype].blastresults = dict()
                sample[self.analysistype].contigs = dict()
                sample[self.analysistype].hits = dict()
                sample[self.analysistype].mismatches = dict()
                sample[self.analysistype].blastrecords = list()
                sample[self.analysistype].range = dict()
                sample[self.analysistype].genespresent = dict()
                # Open blast output csv file
                csvfile = open(sample[self.analysistype].report)
                # Skip header
                csvfile.readline()
                # Open the sequence profile file as a dictionary
                blastdict = DictReader(csvfile, fieldnames=self.fieldnames, dialect='excel-tab')

                # Go through each BLAST result
                for row in blastdict:
                    # Ensure that the hit is full-length, and that the number of mismatches is equal to or lesser
                    # than the supplied cutoff value
                    if int(row['alignment_length']) == self.faidict[row['subject_id']] and \
                                    int(row['mismatches']) <= self.mismatches:
                        # Add the current row to the list for future work
                        sample[self.analysistype].blastrecords.append(row)
                        # Populate the dictionaries with the contig name (e.g. CA_CFIA-515_NODE_1_length_1791),
                        # the gene name (e.g. vtx2a), and the primer name (e.g. vtx2a-R3_1) as required
                        # accounts for primer names with "-" in addition to the terminal "-F" or "-R"
                        try:
                            sample[self.analysistype].blastresults[row['query_id']].add(row['subject_id'])
                            sample[self.analysistype].contigs[row['query_id']].add('-'.join(row['subject_id']
                                                                                            .split('-')[:-1]))
                        except KeyError:
                            sample[self.analysistype].blastresults[row['query_id']] = set()
                            sample[self.analysistype].blastresults[row['query_id']].add(row['subject_id'])
                            sample[self.analysistype].contigs[row['query_id']] = set()
                            sample[self.analysistype].contigs[row['query_id']].add('-'.join(row['subject_id']
                                                                                            .split('-')[:-1]))
                # Check to see if both forward and reverse primers are present for a particular gene within a contig
                for contig, genes in sample[self.analysistype].contigs.items():
                    # Split off the primer details (e.g. vtx2a-R3_1 -> vtx2a-R) from the blast results dictionary in
                    # order to create a searchable list of primers
                    # accounts for primer names with "-" in addition to the terminal "-F" or "-R"
                    reformatted = {'-'.join(['-'.join(x.split('-')[:-1]), x.split('-')[-1][0]])
                                   for x in sample[self.analysistype].blastresults[contig]}
                    # Iterate through the list of genes to check if primers are present
                    for gene in genes:
                        # Add -F and -R to the gene, and ensure that both options are in the reformatted list of genes
                        # e.g. vtx2a-F and vtx2a-R in [vtx1a-R, vtx2c-F ,vtx2a-F, vtx2a-R]
                        if gene + '-F' in reformatted and gene + '-R' in reformatted:
                            # Extract the precise primers from the dictionary e.g. vtx2a use to
                            # find vtx2a-F2_4 (forward) and vtx2a-R3_1 (reverse)
                            forwardprimers = [primer for primer in sample[self.analysistype].blastresults[contig]
                                              if gene == primer.split('-F')[0]]
                            reverseprimers = [primer for primer in sample[self.analysistype].blastresults[contig]
                                              if gene == primer.split('-R')[0]]
                            # Concatenate the lists
                            primers = forwardprimers + reverseprimers
                            # Populate the dictionary with the primers
                            try:
                                sample[self.analysistype].hits[contig].append(primers)
                            except KeyError:
                                sample[self.analysistype].hits[contig] = list()
                                sample[self.analysistype].hits[contig].append(primers)

                            for record in sample[self.analysistype].blastrecords:
                                for primer in primers:
                                    if record['query_id'] == contig and record['subject_id'] == primer:
                                        # Populate the dictionary with the primers
                                        try:
                                            sample[self.analysistype].mismatches[contig][gene]\
                                                .update({primer: int(record['mismatches'])})
                                        except KeyError:
                                            try:
                                                sample[self.analysistype].mismatches[contig][gene] = dict()
                                                sample[self.analysistype].mismatches[contig][gene] = \
                                                    {primer: int(record['mismatches'])}
                                            except KeyError:
                                                sample[self.analysistype].mismatches[contig] = dict()
                                                sample[self.analysistype].mismatches[contig][gene] = dict()
                                                sample[self.analysistype].mismatches[contig][gene] = \
                                                    {primer: int(record['mismatches'])}
                # Use query the stored blast dictionary to find the location of the hits
                for row in sample[self.analysistype].blastrecords:
                    try:
                        # Extract the primers corresponding to the contig
                        for primers in sample[self.analysistype].hits[row['query_id']]:
                            # Extract the name of the contig
                            contig = row['query_id']
                            # Iterate through the forward and reverse primers
                            for primer in primers:
                                # If the primer is present in the current row, then this is the row of interest
                                if row['subject_id'] == primer:
                                    # Split off the primer direction and numbering
                                    # accounts for primer names with "-" in addition to the terminal "-F" or "-R"
                                    gene = '-'.join(primer.split('-')[:-1])
                                    # Populate a dictionary for storing the genes present - will be used in creating
                                    # the report
                                    try:
                                        sample[self.analysistype].genespresent[contig].add(gene)
                                    except KeyError:
                                        sample[self.analysistype].genespresent[contig] = set()
                                        sample[self.analysistype].genespresent[contig].add(gene)
                                    # Populate the range of the hit - the forward primer will have a -F an the name
                                    if '-F' in primer:
                                        # Determine if the sequence is the reverse complement - based on the fact that
                                        # this is the forward primer, if the contig is reversed, then the primer
                                        # (subject) will be reversed.
                                        if int(row['subject_start']) > int(row['subject_end']):
                                            # For reversed sequences, take the larger value of the start and stop
                                            data = max(int(row['query_start']), int(row['query_end']))
                                        else:
                                            # Otherwise take the smaller value
                                            data = min(int(row['query_start']), int(row['query_end']))
                                        # Add the appropriately calculated value to the range dictionary
                                        try:
                                            sample[self.analysistype].range[contig][gene].add(data)
                                        except KeyError:
                                            try:
                                                sample[self.analysistype].range[contig][gene] = set()
                                                sample[self.analysistype].range[contig][gene].add(data)
                                            except KeyError:
                                                sample[self.analysistype].range[contig] = dict()
                                                sample[self.analysistype].range[contig][gene] = set()
                                                sample[self.analysistype].range[contig][gene].add(data)

                                    # Similar to the forward primer, except reverse the min() and max()
                                    elif '-R' in primer:
                                        if int(row['subject_start']) < int(row['subject_end']):
                                            data = min(int(row['query_start']), int(row['query_end']))
                                        else:
                                            data = max(int(row['query_start']), int(row['query_end']))
                                        # Add the appropriately calculated value to the range dictionary
                                        try:
                                            sample[self.analysistype].range[contig][gene].add(data)
                                        except KeyError:
                                            try:
                                                sample[self.analysistype].range[contig][gene] = set()
                                                sample[self.analysistype].range[contig][gene].add(data)
                                            except KeyError:
                                                sample[self.analysistype].range[contig] = dict()
                                                sample[self.analysistype].range[contig][gene] = set()
                                                sample[self.analysistype].range[contig][gene].add(data)
                    except KeyError:
                        pass
638248
function _get_result(Type) {
	return function(rows) {
		if(!rows) { throw new TypeError("failed to parse result"); }
		var doc = rows.shift();
		if(!doc) { return; }

		if(doc instanceof Type) {
			return doc;
		}

		var obj = {};
		ARRAY(Object.keys(doc)).forEach(function(key) {

			if(key === 'documents') {
				obj['$'+key] = {};
				ARRAY(Object.keys(doc[key])).forEach(function(k) {
					if(is.uuid(k)) {
						obj['$'+key][k] = _get_result(NoPg.Document)([doc[key][k]]);
					} else {
						obj['$'+key][k] = doc[key][k];
					}
				});
				return;
			}

			obj['$'+key] = doc[key];
		});

		_parse_object_expressions(obj);

		return new Type(obj);
	};
}
924478
def _get_result_paths(self, data):
        

        result = {}

        # get the file extension of the reads file (sortmerna
        # internally outputs all results with this extension)
        fileExtension = splitext(self.Parameters['--reads'].Value)[1]

        # at this point the parameter --aligned should be set as
        # sortmerna will not run without it
        if self.Parameters['--aligned'].isOff():
            raise ValueError("Error: the --aligned parameter must be set.")

        # file base name for aligned reads
        output_base = self.Parameters['--aligned'].Value

        # Blast alignments
        result['BlastAlignments'] =\
            ResultPath(Path=output_base + '.blast',
                       IsWritten=self.Parameters['--blast'].isOn())

        # SAM alignments
        result['SAMAlignments'] =\
            ResultPath(Path=output_base + '.sam',
                       IsWritten=self.Parameters['--sam'].isOn())

        # OTU map (mandatory output)
        result['OtuMap'] =\
            ResultPath(Path=output_base + '_otus.txt',
                       IsWritten=self.Parameters['--otu_map'].isOn())

        # FASTA file of sequences in the OTU map (madatory output)
        result['FastaMatches'] =\
            ResultPath(Path=output_base + fileExtension,
                       IsWritten=self.Parameters['--fastx'].isOn())

        # FASTA file of sequences not in the OTU map (mandatory output)
        result['FastaForDenovo'] =\
            ResultPath(Path=output_base + '_denovo' +
                       fileExtension,
                       IsWritten=self.Parameters['--de_novo_otu'].isOn())
        # Log file
        result['LogFile'] =\
            ResultPath(Path=output_base + '.log',
                       IsWritten=self.Parameters['--log'].isOn())

        return result