5190565251

pending completion

Build # 5190565251

Build Type

push

github

Committed by

veghp

Commit Message

Bump to v3.2.11

Run Details

1 of 1 new or added line in 1 file covered. (100.0%)

2966 of 3299 relevant lines covered (89.91%)

0.9 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.33

/dnachisel/biotools/blast_sequence.py

from Bio.Blast import NCBIXML
import os
import tempfile
import subprocess
import time


def blast_sequence(
    sequence,
    blast_db=None,
    subject_sequences=None,
    subject=None,
    word_size=4,
    perc_identity=80,
    num_alignments=1000,
    ungapped=False,
    num_threads=3,
    culling_limit=None,
    e_value=None,
    task=None,
    dust="no",
):
    """Return a Biopython BLAST record of the given sequence BLASTed
    against the provided database.

    Parameters
    ----------

    sequence
      An ATGC sequence

    Examples
    --------

    >>> blast_record = blast_sequence("ATTGTGCGTGTGTGCGT", "blastdb/ecoli")
    >>> for alignment in blast_record.alignments:
    >>>     for hit in alignment.hsps:
    >>>         print (hit.identities)
    """

    xml_file, xml_name = tempfile.mkstemp(".xml")
    fasta_file, fasta_name = tempfile.mkstemp(".fa")
    with open(fasta_name, "w+") as f:
        f.write(">seq\n" + sequence)

    close_subject = False
    remove_subject = False

    if subject is not None:
        close_subject = True

    if subject_sequences is not None:
        close_subject = True
        remove_subject = True
        subject_file, subject = tempfile.mkstemp(".fa")
        if isinstance(subject_sequences[0], str):
            subject_sequences = [
                ("%06d" % i, seq) for i, seq in enumerate(subject_sequences)
            ]
        fasta_content = "\n".join(
            [">%s\n%s" % name_sequence for name_sequence in subject_sequences]
        )
        with open(subject, "w+") as f:
            f.write(fasta_content)

    def parameter_if_not_none(label, param):
        return [label, str(param)] if param else []

    command = [
        "blastn",
        "-out",
        xml_name,
        "-outfmt",
        "5",
        "-max_target_seqs",
        str(num_alignments),
        "-query",
        fasta_name,
        "-word_size",
        str(word_size),
        "-num_threads",
        str(num_threads),
        "-perc_identity",
        str(perc_identity),
    ]
    command += (
        (["-db", blast_db] if subject is None else ["-subject", subject])
        + parameter_if_not_none("-dust", dust)
        + parameter_if_not_none("-evalue", e_value)
        + parameter_if_not_none("-culling_limit", culling_limit)
        + parameter_if_not_none("-task", task)
    )
    if ungapped:
        command += ["-ungapped"]
    
    process = subprocess.run(command, stderr=subprocess.PIPE, close_fds=True)
    if process.returncode:
        raise OSError("BLAST failed: %s" % process.stderr)
    with open(xml_name, "r") as f:
        result = list(NCBIXML.parse(f))
    os.fdopen(xml_file, "w").close()
    os.fdopen(fasta_file, "w").close()
    os.remove(xml_name)
    os.remove(fasta_name)
    if close_subject:
        open(subject, "w").close()
        if remove_subject:
            os.remove(subject)
    if len(result) == 1:
        return result[0]
    else:
        return result

1	from Bio.Blast import NCBIXML	1✔
2	import os	1✔
3	import tempfile	1✔
4	import subprocess	1✔
5	import time	1✔
6
7
8	def blast_sequence(	1✔
9	sequence,
10	blast_db=None,
11	subject_sequences=None,
12	subject=None,
13	word_size=4,
14	perc_identity=80,
15	num_alignments=1000,
16	ungapped=False,
17	num_threads=3,
18	culling_limit=None,
19	e_value=None,
20	task=None,
21	dust="no",
22	):
23	"""Return a Biopython BLAST record of the given sequence BLASTed
24	against the provided database.
25
26	Parameters
27	----------
28
29	sequence
30	An ATGC sequence
31
32	Examples
33	--------
34
35	>>> blast_record = blast_sequence("ATTGTGCGTGTGTGCGT", "blastdb/ecoli")
36	>>> for alignment in blast_record.alignments:
37	>>> for hit in alignment.hsps:
38	>>> print (hit.identities)
39	"""
40
41	xml_file, xml_name = tempfile.mkstemp(".xml")	1✔
42	fasta_file, fasta_name = tempfile.mkstemp(".fa")	1✔
43	with open(fasta_name, "w+") as f:	1✔
44	f.write(">seq\n" + sequence)	1✔
45
46	close_subject = False	1✔
47	remove_subject = False	1✔
48
49	if subject is not None:	1✔
50	close_subject = True	×
51
52	if subject_sequences is not None:	1✔
53	close_subject = True	1✔
54	remove_subject = True	1✔
55	subject_file, subject = tempfile.mkstemp(".fa")	1✔
56	if isinstance(subject_sequences[0], str):	1✔
57	subject_sequences = [	1✔
58	("%06d" % i, seq) for i, seq in enumerate(subject_sequences)
59	]
60	fasta_content = "\n".join(	1✔
61	[">%s\n%s" % name_sequence for name_sequence in subject_sequences]
62	)
63	with open(subject, "w+") as f:	1✔
64	f.write(fasta_content)	1✔
65
66	def parameter_if_not_none(label, param):	1✔
67	return [label, str(param)] if param else []	1✔
68
69	command = [	1✔
70	"blastn",
71	"-out",
72	xml_name,
73	"-outfmt",
74	"5",
75	"-max_target_seqs",
76	str(num_alignments),
77	"-query",
78	fasta_name,
79	"-word_size",
80	str(word_size),
81	"-num_threads",
82	str(num_threads),
83	"-perc_identity",
84	str(perc_identity),
85	]
86	command += (	1✔
87	(["-db", blast_db] if subject is None else ["-subject", subject])
88	+ parameter_if_not_none("-dust", dust)
89	+ parameter_if_not_none("-evalue", e_value)
90	+ parameter_if_not_none("-culling_limit", culling_limit)
91	+ parameter_if_not_none("-task", task)
92	)
93	if ungapped:	1✔
94	command += ["-ungapped"]	1✔
95
96	process = subprocess.run(command, stderr=subprocess.PIPE, close_fds=True)	1✔
97	if process.returncode:	1✔
98	raise OSError("BLAST failed: %s" % process.stderr)	×
99	with open(xml_name, "r") as f:	1✔
100	result = list(NCBIXML.parse(f))	1✔
101	os.fdopen(xml_file, "w").close()	1✔
102	os.fdopen(fasta_file, "w").close()	1✔
103	os.remove(xml_name)	1✔
104	os.remove(fasta_name)	1✔
105	if close_subject:	1✔
106	open(subject, "w").close()	1✔
107	if remove_subject:	1✔
108	os.remove(subject)	1✔
109	if len(result) == 1:	1✔
110	return result[0]	1✔
111	else:
112	return result	×

Edinburgh-Genome-Foundry / DnaChisel / 5190565251

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous