-Passes all test cases for 'align' mode (#8)

Author: kwade4

poplars/ref_genomes/K03455-protein.fasta

@@ -65,16 +65,16 @@ IGVTRQRRARNGASRS
 >Tat|HIVHXB2CG
 MEPVDPRLEPWKHPGSQPKTACTNCYCKKCCFHCQVCFITKALGISYGRKKRRQRRRAHQNSQTHQASLSKQPTSQPRGD
 PTGPKESKKKVERETETDPFD
->Tat(exon 1)|HIVHXB2CG
+>Tat(exon1)|HIVHXB2CG
 MEPVDPRLEPWKHPGSQPKTACTNCYCKKCCFHCQVCFITKALGISYGRKKRRQRRRAHQNSQTHQASLSK
->Tat(exon 2)|HIVHXB2CG
+>Tat(exon2)|HIVHXB2CG
 PTSQPRGDPTGPKESKKKVERETETDPFD
 >Rev|HIVHXB2CG
 MAGRSGDSDEELIRTVRLIKLLYQSNPPPNPEGTRQARRNRRRRWRERQRQIHSISERILGTYLGRSAEPVPLQLPPLER
 LTLDCNEDCGTSGTQGVGSPQILVESPTVLESGTKE
->Rev(exon 1)|HIVHXB2CG
+>Rev(exon1)|HIVHXB2CG
 MAGRSGDSDEELIRTVRLIKLLYQS
->Rev(exon 2)|HIVHXB2CG
+>Rev(exon2)|HIVHXB2CG
 PPPNPEGTRQARRNRRRRWRERQRQIHSISERILGTYLGRSAEPVPLQLPPLERLTLDCNEDCGTSGTQGVGSPQILVES
 PTVLESGTKE
 >Vpu|HIVHXB2CG

poplars/ref_genomes/K03455_genome_coordinates.csv

@@ -14,10 +14,8 @@ RNase,3870,4229
 Integrase,4230,5096
 Vif,5041,5619
 Vpr,5559,5850
-Tat,5831,8469
 Tat(exon1),5831,6045
 Tat(exon2),8379,8469
-Rev,5970,8653
 Rev(exon1),5970,6045
 Rev(exon2),8739,8653
 Vpu,6062,6310

poplars/ref_genomes/K03455_protein_coordinates.csv

@@ -12,10 +12,8 @@ RNase,1096,1215
 Integrase,1216,1503
 Vif,1504,1695
 Vpr,1696,1791
-Tat,1792,1892
 Tat(exon1),1792,1862
 Tat(exon2),1863,1892
-Rev(with intron),1893,2008
 Rev(exon1),1893,1918
 Rev(exon2),1919,2008
 Vpu,2009,2090

poplars/ref_genomes/M33262-protein.fasta

@@ -65,20 +65,20 @@ MHCKKGCRCLGEGHGAGGWRPGPPPPPPPGLA
 >Vpr|SIVMM239
 MEERPPENEGPQREPWDEWVVEVLEELKEEALKHFDPRLLTALGNHIYNRHGDTLEGAGELIRILQRALFMHFRGGCIHS
 RIGQPGGGNPLSAIPPSRSML
->Tat(full protein)|SIVMM239
+>Tat|SIVMM239
 METPLREQENSLESSNERSSCISEADASTPESANLGEEILSQLYRPLEACYNTCYCKKCCYHCQFCFLKKGLGICYEQSR
 KRRRTPKKAKANTSSASNKPISNRTRHCQPEKAKKETVEKAVATAPGLGR
->Tat(exon 1)|SIVMM239
+>Tat(exon1)|SIVMM239
 METPLREQENSLESSNERSSCISEADASTPESANLGEEILSQLYRPLEACYNTCYCKKCCYHCQFCFLKKGLGICYEQSR
 KRRRTPKKAKANTSSASN
->Tat(exon 2)|SIVMM239
+>Tat(exon2)|SIVMM239
 PISNRTRHCQPEKAKKETVEKAVATAPGLGR
 >Rev|SIVMM239
 MSNHEREEELRKRLRLIHLLHQTNPYPTGPGTANQRRQRKRRWRRRWQQLLALADRIYSFPDPPTDTPLDLAIQQLQNLA
 IESIPDPPTNTPEALCDPTEDSRSPQD
->Rev(exon 1)|SIVMM239
+>Rev(exon1)|SIVMM239
 MSNHEREEELRKRLRLIHLLHQT
->Rev(exon 2)|SIVMM239
+>Rev(exon2)|SIVMM239
 PYPTGPGTANQRRQRKRRWRRRWQQLLALADRIYSFPDPPTDTPLDLAIQQLQNLAIESIPDPPTNTPEALCDPTEDSRS
 PQD
 >Env(gp160)|SIVMM239

poplars/ref_genomes/M33262_genome_coordinates.csv

@@ -15,10 +15,8 @@ Integrase,4785,5666
 Vif,5596,6240
 Vpx,6068,6406
 Vpr,6407,6712
-Tat,6558,9158
 Tat(exon1),6558,6853
 Tat(exon2),9062,9158
-Rev,6784,9315
 Rev(exon1),6784,6853
 Rev(exon2),9062,9315
 Env(gp160),6860,9499

poplars/ref_genomes/M33262_protein_coordinates.csv

@@ -13,10 +13,8 @@ Integrase,1237,1529
 Vif,1530,1743
 Vpx,1744,1855
 Vpr,1856,1956
-Tat,1957,2086
 Tat(exon1),1957,2055
 Tat(exon2),2056,2086
-Rev,2087,2193
 Rev(exon1),2087,2110
 Rev(exon2),2111,2193
 Env(gp160),2194,3072

poplars/sequence_locator.py

@@ -128,12 +128,12 @@ def make_codon_aln(self):
         if self.nt_seq is not None and self.aa_seq is not None:
             codon_aln = []
             codons = [''.join(t) for t in zip(*[iter(self.nt_seq)] * 3)]
-            for i in range(len(codons)):
+            for aa, codon in zip(self.aa_seq, codons):
                 # Check if stop codon
-                if codons[i] == 'TAA' or codons[i] == 'TGA' or codons[i] == 'TAG':
+                if codon == 'TAA' or codon == 'TGA' or codon == 'TAG':
                     codon_aln.append('-*-')
                 else:
-                    codon_aln.append('-{}-'.format(self.aa_seq[i]))
+                    codon_aln.append('-{}-'.format(aa))
             self.codon_aln = ''.join(codon_aln)
         return self.codon_aln
 
@@ -194,10 +194,10 @@ def set_regions(virus, base, nt_coords, nt_seq, aa_coords, aa_seq):
     :param base: The base of the sequence
     :param nt_coords: Path to the csv file containing the global coordinates of the nucleotide region.
             The file stream has one genomic entry per line and has the following format: region_name,start,end
-    :param nt_seq: The file stream containing the reference nucleotide sequence in read mode
+    :param nt_seq: The nucleotide sequence
     :param aa_coords: Path to the csv file containing the global coordinates of the protein region.
             The file stream has one genomic entry per line and has the following format: region_name,start,end
-    :param aa_seqe: The file stream containing the reference nucleotide sequence in read mode
+    :param aa_seq: A list of lists containing the amino acid sequences
     :return: A list of GenomeRegions
     """
 
@@ -225,26 +225,24 @@ def set_regions(virus, base, nt_coords, nt_seq, aa_coords, aa_seq):
             genome_regions.append(seq_region)
 
         # Parse protein coordinates file
-        prot_names = []
         prot_coords = []
         for aa_line in aa_coords:
             aa_line = aa_line.strip()
             aa_line = aa_line.split(',')
             prot_coords.append([int(aa_line[1]), int(aa_line[2])])
-            prot_names.append(aa_line[0])
 
-        for i in range(len(prot_names)):
+        for i, coords in enumerate(prot_coords):
             for seq_region in genome_regions:
-                if prot_names[i] in seq_region.region_name:
+                if aa_seq[i][0].startswith(seq_region.region_name):
                     # Set global and local protein coordinates
-                    seq_region.set_coords(prot_coords[i], 'prot')
-                    seq_region.set_seq_from_ref(aa_seq, 'prot')
-
+                    seq_region.set_coords(coords, 'prot')
+                    seq_region.set_sequence(aa_seq[i][1], 'prot')
                     seq_region.set_pos_from_pstart(virus)
+                    seq_region.make_codon_aln()
 
     else:
         # Parse protein region coordinates file
-        for aa_line in aa_coords:
+        for i, aa_line in enumerate(aa_coords):
             aa_line = aa_line.strip()
             aa_line = aa_line.split(',')
             prot_coords = [int(aa_line[1]), int(aa_line[2])]
@@ -253,7 +251,7 @@ def set_regions(virus, base, nt_coords, nt_seq, aa_coords, aa_seq):
 
             # Set global and local nucleotide coordinates
             seq_region.set_coords(prot_coords, 'prot')
-            seq_region.set_seq_from_ref(aa_seq, 'prot')
+            seq_region.set_sequence(aa_seq[i][1], 'prot')
 
             # Set relative positions
             seq_region.set_pos_from_cds(virus)
@@ -268,17 +266,17 @@ def set_regions(virus, base, nt_coords, nt_seq, aa_coords, aa_seq):
         for nt_line in nt_coords:
             nt_line = nt_line.strip()
             nt_line = nt_line.split(',')
-            nucl_coords.append([int(nt_line[1]), int(nt_line[2])])
             nucl_names.append(nt_line[0])
+            nucl_coords.append([int(nt_line[1]), int(nt_line[2])])
 
-        for i in range(len(nucl_names)):
+        for i, name in enumerate(nucl_names):
             for seq_region in genome_regions:
-                if nucl_names[i] in seq_region.region_name:
+                if name.startswith(seq_region.region_name):
                     # Set global and local protein coordinates
                     seq_region.set_coords(nucl_coords[i], 'nucl')
                     seq_region.set_seq_from_ref(nt_seq, 'nucl')
-
                     seq_region.set_pos_from_pstart(virus)
+                    seq_region.make_codon_aln()
 
     return genome_regions
 
@@ -557,7 +555,7 @@ def output_retrieved_region(region, outfile=None):
         print("\tNucleotide Sequence:")
         seq_lines = [region.nt_seq[i:i + 60] for i in range(0, len(region.nt_seq), 60)]
         for line in seq_lines:
-            print('\t\t{}\n'.format(line))
+            print('\t\t{}'.format(line))
 
         if region.aa_seq:
             print("\tProtein Sequence:")
@@ -739,8 +737,8 @@ def retrieve(virus, base, ref_regions, region, qstart=1, qend='end'):
             query_region = GenomeRegion(region)
 
             # Set local and global coordinates
-            query_region.set_coords([qstart, qend], base)
-            global_coords = query_region.local_to_global_index(ref_region, [qstart, qend], base)
+            # query_region.set_coords([qstart, qend], base)
+            global_coords = GenomeRegion.local_to_global_index(ref_region, [qstart, qend], base)
             query_region.set_coords(global_coords, base)
 
             # Set sequences protein and nucleotide sequences
@@ -878,7 +876,7 @@ def main():
     # Ensure proper configuration files are set
     configs = handle_args(args.virus, args.base, args.ref_nt, args.nt_coords, args.ref_aa, args.aa_coords)
     ref_nt_seq = configs[0][0][1]
-    ref_aa_seq = configs[1][0][1]
+    ref_aa_seq = configs[1]
     nt_coords = configs[2]
     aa_coords = configs[3]
     reference_sequence = configs[4]

poplars/tests/fixtures/siv-test-proteins.fasta

@@ -1,8 +1,4 @@
->Rev|SIVMM239
-MSNHEREEELRKRLRLIHLLHQTNPYPTGPGTANQRRQRKRRWRRRWQQLLALADRIYSFPDPPTDTPLDLAIQQLQNLA
-IESIPDPPTNTPEALCDPTEDSRSPQD
->Rev(exon 1)|SIVMM239
+>Rev(exon1)|SIVMM239
 MSNHEREEELRKRLRLIHLLHQT
->Rev(exon 2)|SIVMM239
-PYPTGPGTANQRRQRKRRWRRRWQQLLALADRIYSFPDPPTDTPLDLAIQQLQNLAIESIPDPPTNTPEALCDPTEDSRS
-PQD
+>Rev(exon2)|SIVMM239
+PYPTGPGTANQRRQRKRRWRRRWQQLLALADRIYSFPDPPTDTPLDLAIQQLQNLAIESIPDPPTNTPEALCDPTEDSRSPQD

poplars/tests/fixtures/siv_test_nt_coords.csv

@@ -1,3 +1,2 @@
-Rev(with intron),6784,9315
 Rev(exon1),6784,6853
 Rev(exon2),9062,9315

poplars/tests/fixtures/siv_test_prot_coords.csv

@@ -1,3 +1,2 @@
-Rev,2087,2193
 Rev(exon1),2087,2110
 Rev(exon2),2111,2193