primalutils/quick_align.diff.txt at master · joshquick/primalutils · GitHub

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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
diff --git a/ncov-random-scripts/quick_align.py b/Users/Quickj-admin/bioinfo/tyson_scheme/Feb-2022/multi_fasta_mask/ncov-random-scripts/quick_align.py
index ee6d4a0..5937c7b 100644
--- a/ncov-random-scripts/quick_align.py
+++ b/Users/Quickj-admin/bioinfo/tyson_scheme/Feb-2022/multi_fasta_mask/ncov-random-scripts/quick_align.py
@@ -7,41 +7,41 @@ import parasail
 import argparse
 import textwrap as tw

+
 def get_sequence(file):
     fasta = pysam.FastxFile(file)
     reference = None
     for record in fasta:
         # For this narrow application there should be exactly one entry in the reference file
-        assert(reference is None)
+        assert reference is None
         reference = record
     return reference

+
 def get_alignment_parasail(reference_genome, input_genome):
-
+
     # the dna full matrix supports ambiguity codes, although "N"s are not given free mismatches as we might like
     # the alignments appear good enough for our purpose however
-    # do not penalise gaps at end of database
-    result = parasail.sg_dx_trace_striped_32(input_genome.sequence, reference_genome.sequence, 10, 1, parasail.dnafull)
+    result = parasail.nw_trace_striped_32(
+        input_genome.sequence, reference_genome.sequence, 10, 1, parasail.dnafull
+    )
     traceback = result.traceback
+
     return traceback.ref, traceback.comp, traceback.query

+
 def alignment2vcf(reference_name, reference_aligned, query_aligned):

+    reference_position = 0
+    query_position = 0
+
     print("##fileformat=VCFv4.2")
     print('##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">')
     print("##contig=<ID=%s>" % (reference_name))
     print("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tsample")

-
-    # only write out out records for overlapping region
-    q_start = len(query_aligned) - len(query_aligned.lstrip('-'))
-    q_end = len(query_aligned) - (len(query_aligned) - len(query_aligned.rstrip('-')))
-
-    i = q_start
-    query_position = q_start
-    reference_position = q_start
-    n = q_end
-
+    i = 0
+    n = len(reference_aligned)
     while i < n:
         if query_aligned[i] == reference_aligned[i]:
             reference_position += 1
@@ -63,16 +63,24 @@ def alignment2vcf(reference_name, reference_aligned, query_aligned):
             offset = 0
             if q_gaps > 0 or r_gaps > 0:
                 # append a single base to the start
-                q_sub = query_aligned[i-1] + q_sub.replace("-", "")
-                r_sub = reference_aligned[i-1] + r_sub.replace("-", "")
+                q_sub = query_aligned[i - 1] + q_sub.replace("-", "")
+                r_sub = reference_aligned[i - 1] + r_sub.replace("-", "")
                 offset = 1

             # Record the difference
-            print("%s\t%d\t.\t%s\t%s\t.\t.\t.\tGT\t1" % (reference_name, reference_position - offset + 1, r_sub.upper(), q_sub.upper()))
+            print(
+                "%s\t%d\t.\t%s\t%s\t.\t.\t.\tGT\t1"
+                % (
+                    reference_name,
+                    reference_position - offset + 1,
+                    r_sub.upper(),
+                    q_sub.upper(),
+                )
+            )

             # update counters
-            reference_position += (j - i - r_gaps)
-            query_position += (j - i - q_gaps)
+            reference_position += j - i - r_gaps
+            query_position += j - i - q_gaps
             i = j


@@ -80,11 +88,13 @@ def main():
     """
     Main method for script
     """
-    description = 'Align a pair of genomes and write the results in various formats'
+    description = "Align a pair of genomes and write the results in various formats"
     parser = argparse.ArgumentParser(description=description)
-    parser.add_argument('-g', '--genome', help='consensus genome FASTA file to process')
-    parser.add_argument('-r', '--reference-genome', help='fasta file containing the reference genome')
-    parser.add_argument('-o', '--output-mode', default="differences")
+    parser.add_argument("-g", "--genome", help="consensus genome FASTA file to process")
+    parser.add_argument(
+        "-r", "--reference-genome", help="fasta file containing the reference genome"
+    )
+    parser.add_argument("-o", "--output-mode", default="differences")

     if len(sys.argv) <= 1:
         parser.print_help(sys.stderr)
@@ -93,15 +103,17 @@ def main():

     reference_genome = get_sequence(file=args.reference_genome)
     input_genome = get_sequence(file=args.genome)
-    (reference_aligned, comparison_aligned, query_aligned) = get_alignment_parasail(reference_genome, input_genome)
+    (reference_aligned, comparison_aligned, query_aligned) = get_alignment_parasail(
+        reference_genome, input_genome
+    )

     if args.output_mode == "differences":
         columns = 120
         for start in range(0, len(reference_aligned), columns):

-            r = reference_aligned[start:start+columns]
-            c = comparison_aligned[start:start+columns]
-            q = query_aligned[start:start+columns]
+            r = reference_aligned[start : start + columns]
+            c = comparison_aligned[start : start + columns]
+            q = query_aligned[start : start + columns]

             if r != q:
                 print(start, r)
@@ -114,5 +126,6 @@ def main():
     elif args.output_mode == "vcf":
         alignment2vcf(reference_genome.name, reference_aligned, query_aligned)

-if __name__ == '__main__':
+
+if __name__ == "__main__":
     main()