Added securities in case inputted alignments are incompatible with the graph

Author: Roland Faure

graphunzip.py

@@ -320,9 +320,11 @@ def main():
 
             print("WARNING: all interaction matrices are empty, GraphUnzip does not do anything")
 
+        print("\nDone untangling, now merging all contigs that can be merged")
         merge_adjacent_contigs(segments)
 
-        # now exporting the output        
+        # now exporting the output  
+        print("Now exporting the result")
         io.export_to_GFA(segments, gfaFile, exportFile=outFile, merge_adjacent_contigs=merge)
 
         if fastaFile != "None":

input_output.py

@@ -12,6 +12,7 @@
 import pickle #for writing files and reading them
 import re #to find all numbers in a mixed number/letters string (such as 31M1D4M), to split on several characters (<> in longReads_interactionMatrix)
 import shutil #to remove directories
+import sys
 
 from segment import Segment
 from segment import compute_copiesNumber
@@ -150,9 +151,8 @@ def read_TSV(tsv_file, names, lines):
                     alignment += contig[:-1]
 
                     if contig[:-1] not in names :
-                        print("Problem: ", line)
-                        while True :
-                            rien = 0
+                        print("ERROR: while reading the .tsv, I am coming across a contig that was not in the .gfa, namely ", contig[:-1], ". I recommend you check that you are using the same GFA that you aligned the long reads on.")
+                        sys.exit()
                 lines += [alignment]
 
 

solve_with_long_reads.py

@@ -12,6 +12,7 @@
 from determine_multiplicity import determine_multiplicity
 from input_output import read_GAF
 from input_output import read_TSV
+from segment import find_this_link
 
 import segment
 
@@ -57,9 +58,9 @@ def bridge_with_long_reads(segments, names, copiesnumber, gafFile, supported_lin
         read_GAF(gafFile, 0.7, 0.1, lines)
         print("Finished going through the gaf file.")
     elif '.tsv' in gafFile :
-        print("Reading the gpa file...")
+        print("Reading the tsv file...")
         read_TSV(gafFile, names, lines)
-        print("Finished going through the gpa file.")
+        print("Finished going through the tsv file.")
     else :
         print("ERROR: input format of mapped read not recognized. It should be .gfa or .gpa")
         sys.exit()
@@ -71,7 +72,7 @@ def bridge_with_long_reads(segments, names, copiesnumber, gafFile, supported_lin
     longContigs = [True for i in range(len(names))] #then all contigs that are in the middle of a read will be marked as False
     bridges = [[[],[]] for i in range(len(haploidContigs))] #bridges is a list inventoring at index haploidCOntigsNames[seg.names[0]] all the links left and right of the contig, supported by the gaf
     minimum_supported_links = sparse.lil_matrix((len(names)*2, len(names)*2)) #minimum_supported links is the list of all links between different contigs found at least once in the gaf file
-    inventoriate_bridges(lines, bridges, minimum_supported_links, haploidContigsNames, longContigs, names)
+    inventoriate_bridges(lines, bridges, minimum_supported_links, haploidContigsNames, longContigs, names, segments)
 
     #now, from all the bridges, build consensus bridges
     consensus_bridges = [['',''] for i in range(len(haploidContigs))] #consensus bridge is essentially the same as bridges, except there is only one bridge left at each side for each contig
@@ -124,74 +125,92 @@ def bridge_with_long_reads(segments, names, copiesnumber, gafFile, supported_lin
 
 #input : a list of alignments of a gaf file
 #output : the completed bridges list, with for each haploid contig a list of what was found left and right of the contig
-def inventoriate_bridges(lines, bridges, minimum_supported_links, haploidContigsNames, longContigs, names) :
+def inventoriate_bridges(lines, bridges, minimum_supported_links, haploidContigsNames, longContigs, names, segments) :
 
 
     for l, line in enumerate(lines) :      
 
         if (l+1) % 1000 == 0 :
-            print("Inventoried ", l+1, " long reads over ", len(lines))
+            print("Inventoried ", l+1, " long reads over ", len(lines), end = '\r')
 
         contigs = re.split('[><]' , line)
         orientations = "".join(re.findall("[<>]", line))
         del contigs[0] #because the first element is always ''
-    
+        
+        #first go through the alignment to make sure it is possible on the gfa
+        possible = True
         for c, contig in enumerate(contigs) :
-            
             if c>0 :
-                minimum_supported_links[2*names[contigs[c-1]] + '<>'.index(orientations[c-1]) , 2*names[contigs[c]] + '><'.index(orientations[c])] = 1
-                minimum_supported_links[2*names[contigs[c]] + '><'.index(orientations[c]), 2*names[contigs[c-1]] + '<>'.index(orientations[c-1])] = 1
+                or1 = '<>'.index(orientations[c-1])
+                or2 = '><'.index(orientations[c])
+                #check if the link actually exists (it should, if the aligner did its job correctly, but apparently sometimes SPAligner behaves strangely)
+                if -1 == find_this_link(segments[names[contig]], or2, segments[names[contigs[c-1]]].links[or1], segments[names[contigs[c-1]]].otherEndOfLinks[or1]) :
+                    print ("WARNING: discrepancy between what's found in the alignment files and the inputted GFA graph. Link ", contigs[c-1:c+1], orientations[c-1:c+1], " not found in the gfa")
+                    possible = False
+                    
+        #then, only inventoriate the bridge if it is possible with respect to the graph     
+        if  possible :
 
-            if c > 0 and c < len(contigs) - 1 :
-                longContigs[names[contig]] = False
-                                        
-            if contig in haploidContigsNames :
-                
+            for c, contig in enumerate(contigs) :
 
-                if orientations[c] == ">" :
-                    r = 0
-                    #first look at what contigs are left of the contig of interest
-                    bridges[haploidContigsNames[contig]][1] +=  [""]
-                    for c2 in range(c+1, len(contigs)) :
-                        
-                        bridges[haploidContigsNames[contig]][1][-1] += orientations[c2] + contigs[c2]
-                        
-                        # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
-                        #     break
-                        
-                    #then look at what's left of the contig of interest (so mirror the orientations)
-                    bridges[haploidContigsNames[contig]][0] +=  [""]
-                    for c2 in range(c-1 , -1, -1) :
-                        
-                        if orientations[c2] == '>' :
-                            bridges[haploidContigsNames[contig]][0][-1] += '<' + contigs[c2]
-                        else :
-                            bridges[haploidContigsNames[contig]][0][-1] += '>' + contigs[c2]
+                if c>0 :
+                    
+                    or1 = '<>'.index(orientations[c-1])
+                    or2 = '><'.index(orientations[c])
+                    
+                    minimum_supported_links[2*names[contigs[c-1]] + or1, 2*names[contigs[c]] + or2] = 1
+                    minimum_supported_links[2*names[contigs[c]] + or2, 2*names[contigs[c-1]] + or1] = 1
+        
+                if c > 0 and c < len(contigs) - 1 :
+                    longContigs[names[contig]] = False
+                                            
+                if contig in haploidContigsNames :
+                    
+                    
+                    if orientations[c] == ">" :
+                        r = 0
+                        #first look at what contigs are left of the contig of interest
+                        bridges[haploidContigsNames[contig]][1] +=  [""]
+                        for c2 in range(c+1, len(contigs)) :
 
-                        # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
-                        #     break
+                            bridges[haploidContigsNames[contig]][1][-1] += orientations[c2] + contigs[c2]
 
-                else :
-                    
-                    #first look at what contigs are left of the contig of interest
-                    bridges[haploidContigsNames[contig]][0] +=  [""]
-                    for c2 in range(c+1, len(contigs)) :
-                        bridges[haploidContigsNames[contig]][0][-1] += orientations[c2] + contigs[c2]
-                        # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
-                        #     break
-                        
-                    #then look at what's left of the contig of interest (so mirror the orientations)
-                    bridges[haploidContigsNames[contig]][1] +=  [""]
-                    for c2 in range(c-1 , -1, -1 ) :
+                            # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
+                            #     break
+                            
+                        #then look at what's left of the contig of interest (so mirror the orientations)
+                        bridges[haploidContigsNames[contig]][0] +=  [""]
+                        for c2 in range(c-1 , -1, -1) :
+                            
+                            if orientations[c2] == '>' :
+                                bridges[haploidContigsNames[contig]][0][-1] += '<' + contigs[c2]
+                            else :
+                                bridges[haploidContigsNames[contig]][0][-1] += '>' + contigs[c2]
+                                
+                            # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
+                            #     break
+                                
+                    else :
 
-                        if orientations[c2] == '>' :
-                            bridges[haploidContigsNames[contig]][1][-1] += '<' + contigs[c2]
-                        else :
-                            bridges[haploidContigsNames[contig]][1][-1] += '>' + contigs[c2]
+                        #first look at what contigs are left of the contig of interest
+                        bridges[haploidContigsNames[contig]][0] +=  [""]
+                        for c2 in range(c+1, len(contigs)) :
+                            bridges[haploidContigsNames[contig]][0][-1] += orientations[c2] + contigs[c2]
+                            # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
+                            #     break
 
-                        # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
-                        #     break
-                                                 
+                        #then look at what's left of the contig of interest (so mirror the orientations)
+                        bridges[haploidContigsNames[contig]][1] +=  [""]
+                        for c2 in range(c-1 , -1, -1 ) :
+                            
+                            if orientations[c2] == '>' :
+                                bridges[haploidContigsNames[contig]][1][-1] += '<' + contigs[c2]
+                            else :
+                                bridges[haploidContigsNames[contig]][1][-1] += '>' + contigs[c2]
+                                
+                            # if contigs[c2] in haploidContigsNames : #you can stop the bridge, you've reached the other side
+                            #     break
+                                                     
 #input : list of bridges for each haploid contig
 #output : completed consensus_bridges, where there is max one bridge at each end of contig
 def build_consensus_bridges(consensus_bridges, bridges, names, haploidContigs, haploidContigsNames):
@@ -206,7 +225,7 @@ def build_consensus_bridges(consensus_bridges, bridges, names, haploidContigs, h
     for c in range(len(bridges)) :
 
         if (c)%100 == 0 :
-            print("consensused ", c, " bridges out of ", len(consensus_bridges))
+            print("consensused ", c, " bridges out of ", len(consensus_bridges), end = '\r')
 
         localContigs = [ [ re.split('[><]' , bridges[c][j][k])[1:] for k in range(len(bridges[c][j])) ] for j in range(2)]
         localOrientations = [ [ "".join(re.findall("[<>]", bridges[c][j][k])) for k in range(len(bridges[c][j])) ] for j in range(2)]
@@ -410,7 +429,7 @@ def unzip_graph_with_bridges(segments, non_overlapping_bridges, copiesnumber, ha
     for se in range(l) :
 
         if (se)%1000 == 0 :
-            print("Processed ", se, " contigs out of ", l, ", while untangling with long reads")
+            print("Processed ", se, " contigs out of ", l, ", while untangling with long reads", end = '\r')
         s = segments[se]
 
         if s.names[0] in haploidContigsNames :
@@ -524,7 +543,7 @@ def unzip_graph_with_bridges(segments, non_overlapping_bridges, copiesnumber, ha
 
                                 for alreadyDuplicatedContig in segments[oldContigsIndices[c-1]].links[end0] :
                                     if alreadyDuplicatedContig.names[0] == contigs[c]:
-                                        segment.add_link(segments[newContigsIndices[-1]] , end0, alreadyDuplicatedContig, end1)
+                                        segment.add_link(segments[newContigsIndices[-1]] , end0, alreadyDuplicatedContig, end1, CIGAR)
 
 
                         else :