1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
from cgelib.alignment.aligner import BlastNAligner, KMAAligner
import matplotlib.pyplot as plt
import os
import json
def runblast_finder(data, output, db):
""" Function for running the BlastNAligner. It takes:
- Data: the assembled fasta file.
- Output: Path where to save the output files.
- Db: List of files that serve as 'template' for blast. For example:
['/path/to/aminoglycoside.fsa','/path/to/beta-lactam.fsa']
"""
# Declare an instance of the blast aligner, like:
# aligner = BlastNAligner(result_file="XML")
# Set the parameters for blast
# Run the aligner
# Fit the alignment
# Blast do not report directly the coverage. Calculate it with the following
# formula: perc_coverage = ((((1+int(query_end)-int(query_start)) - int(gaps))
# / float(alignment_length)) * 100)
# If you wanna save the results in a json file (recommended if you are
# starting), read the results and save the results
def runkma_finder(data, output, db):
""" Function for running the KMAAligner. It takes:
- Data: If single-end reads, a string (example: 'path/to/rawread.fq').
If paired-end reads, a list (example: ['path/to/rawread1.fq',
'path/to/rawread2.fq'])
- Output: Path where to save the output files.
- Db: List of files (without extension) that serve as 'template' for
kma (indexed files). For example: ['/path/to/aminoglycoside',
'/path/to/beta-lactam']
"""
# Declare an instance of the kma aligner, like:
# aligner = KMAAligner(result_file=["Result", "Mapstat"])
# Set the parameters for kma
# Run the aligner
# Fit the alignment
# If you wanna save the results in a json file (recommended if you are
# starting and do not want to run the aligner several times), read the
# results and save the results
def graph_fragmentsfound(json_file, output, x, y, diameter=None):
""" Example of how to analyze the alignment from the json file. It creates
a graph that show two/three features of each hit:
- Json_file: Path to the json created by an aligner
(example: '/path/to/json')
- Output: path and name of the file where to save the graph
(example: '/path/to/graph.png')
- x: feature of the hit on the x axis
- y: feature of the hit on the y axis
- diameter: feature of the hit on the diameter of the dot
representing the hit. Not necessary
"""
with open(json_file, "r") as read_file:
results = json.load(read_file)
fig, ax = plt.subplots()
for n in results["aln_hits"]:
if diameter is None:
ax.scatter(x=float(results["aln_hits"][n][x]),
y=float(results["aln_hits"][n][y]),
alpha=0.7,
label=results["aln_hits"][n]["templateID"])
else:
ax.scatter(x=float(results["aln_hits"][n][x]),
y=float(results["aln_hits"][n][y]),
s=float(results["aln_hits"][n][diameter]),
alpha=0.7,
label=results["aln_hits"][n]["templateID"])
plt.xlabel(x)
plt.ylabel(y)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.savefig(output, bbox_inches="tight")
def gene_to_phenotype(aln_hits, db):
""" """
pass
# In order to run the finders/graphs functions, fill and uncomment some of the
# next lines:
#runblast_finder(data=?, output=?, db=?)
#runkma_finder(data=?, output=?, db=?)
#graph_fragmentsfound(json_file=?, output=?, x=?, y=?, diameter=?)
|
