Skip to content

Commit b29b9f4

Browse files
yangao07yangao07
authored
Merge pull request #93 from dialvarezs/python-use-qscores
feat: support quality scores in python api
2 parents 7acdac7 + 7a86871 commit b29b9f4

3 files changed

Lines changed: 159 additions & 107 deletions

File tree

‎python/README.md‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -76,9 +76,9 @@ This constructs a multiple sequence alignment handler of pyabpoa, it accepts the
7676
* **extra_f**: second adaptive banding paremete; the number of extra bases added on both sites of the band is *b+f\*L*, where *L* is the length of the aligned sequence; default : **0.01**
7777
* **cons_algrm**: consensus calling algorithm. 'HB': heaviest bunlding, 'MF': most frequent bases; default: **'HB'**
7878

79-
The `msa_aligner` handler provides one method `msa` which performs multiple sequence alignment and takes four arguments:
79+
The `msa_aligner` handler provides one method `msa` which performs multiple sequence alignment and accepts the following arguments:
8080
```
81-
pyabpoa.msa_aligner.msa(seqs, out_cons, out_msa, out_pog='', incr_fn='')
81+
pyabpoa.msa_aligner.msa(seqs, out_cons, out_msa, out_pog='', incr_fn='', qscores=None)
8282
```
8383

8484
* **seqs**: a list variable containing a set of input sequences; **positional**
@@ -88,12 +88,15 @@ pyabpoa.msa_aligner.msa(seqs, out_cons, out_msa, out_pog='', incr_fn='')
8888
* **min_freq**: minimum frequency of each consensus to output (effective when **max_n_cons** > 1); default: **0.3**
8989
* **out_pog**: name of a file (`.png` or `.pdf`) to store the plot of the final alignment graph; default: **''**
9090
* **incr_fn**: name of an existing graph (GFA) or MSA (FASTA) file, incrementally align sequence to this graph/MSA; default: **''**
91+
* **qscores**: optional per-sequence quality information used to weight the consensus graph. Each entry must match the corresponding sequence length and should be a list of integer Phred scores, e.g. `record.letter_annotations["phred_quality"]` from Biopython; default: **None**
92+
93+
When `qscores` is provided, pyabpoa enables quality-weighted consensus (`use_qv`) for that run. This is effective for heaviest-bundling consensus (`cons_algrm='HB'`), matching the CLI `-Q/--use-qual-weight` behavior.
9194

9295
`msa_aligner` also provides three methods for incrementally adding sequences to graph/MSA:
9396

9497
```
95-
pyabpoa.msa_aligner.msa_align(seqs, out_cons, out_msa, max_n_cons=1, min_freq=0.25, incr_fn=b'')
96-
pyabpoa.msa_aligner.msa_add(new_seqs)
98+
pyabpoa.msa_aligner.msa_align(seqs, out_cons, out_msa, max_n_cons=1, min_freq=0.25, incr_fn=b'', qscores=None)
99+
pyabpoa.msa_aligner.msa_add(new_seqs, qscores=None)
97100
pyabpoa.msa_aligner.msa_output()
98101
```
99102

@@ -113,6 +116,7 @@ This class describes the information of the generated consensus sequence and the
113116
* **cons_len**: an array of consensus sequence length(s)
114117
* **cons_seq**: an array of consensus sequence(s)
115118
* **cons_cov**: an array of consensus sequence coverage for each base
119+
* **cons_qv**: an array of consensus quality strings in FASTQ encoding
116120
* **msa_len**: size of each row in the RC-MSA
117121
* **msa_seq**: an array containing `n_seq`+`n_cons` strings that demonstrates the RC-MSA, each consisting of one input sequence and several `-` indicating the alignment gaps.
118122

‎python/example.py‎

Lines changed: 65 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
import pyabpoa as pa
22

3-
#@parameters of msa_aligner:
3+
# @parameters of msa_aligner:
44
# aln_mode='g' # g: global, l: local, e: extension
55
# is_aa=False # set as True if input is amino acid sequence
66
# score_matrix='' # file of score matrix, e.g. HOXD70.mtx/BLOSUM62.mtx
@@ -20,71 +20,95 @@
2020

2121
print("==== First exmaple: 2 consensus sequences ====\n")
2222
# for multiple consensus
23-
seqs=[
24-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT',
25-
'CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT',
26-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT',
27-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT',
28-
'CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT',
29-
'CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT',
30-
'CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT',
31-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT',
32-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT',
33-
'CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT'
34-
]
35-
36-
#@parameters of msa
37-
#seqs: multiple sequences
38-
out_cons=True # generate consensus sequence, set as False to disable
39-
out_msa=True # generate row-column multiple sequence alignment, set as False to disable
40-
#out_pog="example1.png" # generate plot of alignment graph, set None to disable, require `dot` to be installed
23+
seqs = [
24+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
25+
"CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT",
26+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
27+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
28+
"CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT",
29+
"CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT",
30+
"CGATCGATCGATAAAAAAAAAAAAAAAAAAACGATGCATGCATCGATGCATCGATCGATGCATGCAT",
31+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
32+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
33+
"CGATCGATCGATCGATGCATGCATCGATGCATCGATCGATGCATGCAT",
34+
]
35+
36+
# @parameters of msa
37+
# seqs: multiple sequences
38+
out_cons = True # generate consensus sequence, set as False to disable
39+
out_msa = True # generate row-column multiple sequence alignment, set as False to disable
40+
# out_pog="example1.png" # generate plot of alignment graph, set None to disable, require `dot` to be installed
4141
max_n_cons = 2
4242

4343
# multiple sequence alignment for 'seqs'
44-
res=a.msa(seqs, out_cons=out_cons, out_msa=out_msa, max_n_cons=max_n_cons) #, out_pog=out_pog)
44+
res = a.msa(seqs, out_cons=out_cons, out_msa=out_msa, max_n_cons=max_n_cons) # , out_pog=out_pog)
4545

4646
# output result
4747
if out_cons:
4848
for i in range(res.n_cons):
49-
print(">Consensus_sequence_{}".format(i+1))
49+
print(f">Consensus_sequence_{i + 1}")
5050
print(res.cons_seq[i])
5151
if out_msa:
5252
res.print_msa()
5353

5454

5555
print("\n\n==== Second exmaple: 1 consensus sequence ====\n")
56-
seqs=[
57-
'CGTCAATCTATCGAAGCATACGCGGGCAGAGCCGAAGACCTCGGCAATCCA',
58-
'CCACGTCAATCTATCGAAGCATACGCGGCAGCCGAACTCGACCTCGGCAATCAC',
59-
'CGTCAATCTATCGAAGCATACGCGGCAGAGCCCGGAAGACCTCGGCAATCAC',
60-
'CGTCAATGCTAGTCGAAGCAGCTGCGGCAGAGCCGAAGACCTCGGCAATCAC',
61-
'CGTCAATCTATCGAAGCATTCTACGCGGCAGAGCCGACCTCGGCAATCAC',
62-
'CGTCAATCTAGAAGCATACGCGGCAAGAGCCGAAGACCTCGGCCAATCAC',
63-
'CGTCAATCTATCGGTAAAGCATACGCTCTGTAGCCGAAGACCTCGGCAATCAC',
64-
'CGTCAATCTATCTTCAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATC',
65-
'CGTCAATGGATCGAGTACGCGGCAGAGCCGAAGACCTCGGCAATCAC',
66-
'CGTCAATCTAATCGAAGCATACGCGGCAGAGCCGTCTACCTCGGCAATCACGT'
56+
seqs = [
57+
"CGTCAATCTATCGAAGCATACGCGGGCAGAGCCGAAGACCTCGGCAATCCA",
58+
"CCACGTCAATCTATCGAAGCATACGCGGCAGCCGAACTCGACCTCGGCAATCAC",
59+
"CGTCAATCTATCGAAGCATACGCGGCAGAGCCCGGAAGACCTCGGCAATCAC",
60+
"CGTCAATGCTAGTCGAAGCAGCTGCGGCAGAGCCGAAGACCTCGGCAATCAC",
61+
"CGTCAATCTATCGAAGCATTCTACGCGGCAGAGCCGACCTCGGCAATCAC",
62+
"CGTCAATCTAGAAGCATACGCGGCAAGAGCCGAAGACCTCGGCCAATCAC",
63+
"CGTCAATCTATCGGTAAAGCATACGCTCTGTAGCCGAAGACCTCGGCAATCAC",
64+
"CGTCAATCTATCTTCAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATC",
65+
"CGTCAATGGATCGAGTACGCGGCAGAGCCGAAGACCTCGGCAATCAC",
66+
"CGTCAATCTAATCGAAGCATACGCGGCAGAGCCGTCTACCTCGGCAATCACGT",
6767
]
6868

69-
#@parameters of msa
70-
#seqs: multiple sequences
71-
out_cons=True # generate consensus sequence, set as False to disable
72-
out_msa=True # generate row-column multiple sequence alignment, set as False to disable
69+
# @parameters of msa
70+
# seqs: multiple sequences
71+
out_cons = True # generate consensus sequence, set as False to disable
72+
out_msa = True # generate row-column multiple sequence alignment, set as False to disable
7373
# out_pog="example2.png" # generate plot of alignment graph, set None to disable
7474
max_n_cons = 1
7575

7676
# multiple sequence alignment for 'seqs'
77-
res=a.msa(seqs, out_cons=out_cons, out_msa=out_msa, max_n_cons=max_n_cons) #, out_pog=out_pog)
77+
res = a.msa(seqs, out_cons=out_cons, out_msa=out_msa, max_n_cons=max_n_cons) # , out_pog=out_pog)
7878

7979
# output result
8080
if out_cons:
8181
for i in range(res.n_cons):
82-
print(">Consensus_sequence_{}".format(i+1))
82+
print(f">Consensus_sequence_{i + 1}")
8383
print(res.cons_seq[i])
8484
if out_msa:
8585
for i in range(res.n_seq):
86-
print(">Seq_{}".format(i+1))
86+
print(f">Seq_{i + 1}")
8787
print(res.msa_seq[i])
8888
for i in range(res.n_cons):
89-
print(">Consensus_sequence_{}".format(i+1))
90-
print(res.msa_seq[res.n_seq+i])
89+
print(f">Consensus_sequence_{i + 1}")
90+
print(res.msa_seq[res.n_seq + i])
91+
92+
93+
print("\n\n==== Third exmaple: quality-weighted consensus ====\n")
94+
seqs = [
95+
"ACGT",
96+
"ACGT",
97+
"ACGA",
98+
]
99+
100+
# Per-base Phred qualities. This matches the shape returned by
101+
# Biopython's record.letter_annotations["phred_quality"] for FASTQ input.
102+
qscores = [
103+
[40, 40, 40, 40],
104+
[35, 35, 35, 35],
105+
[30, 30, 30, 5],
106+
]
107+
108+
res = a.msa(seqs, out_cons=True, out_msa=False, qscores=qscores)
109+
110+
for i in range(res.n_cons):
111+
print(f"@Consensus_sequence_{i + 1}")
112+
print(res.cons_seq[i])
113+
print("+")
114+
print(res.cons_qv[i])

0 commit comments

Comments
 (0)