New command: obi alignpairedend

python/obitools3/commands/alignpairedend.pyx

@@ -0,0 +1,170 @@
+#cython: language_level=3
+
+from obitools3.apps.progress cimport ProgressBar  # @UnresolvedImport
+from obitools3.dms import DMS
+from obitools3.dms.view.view cimport View
+from obitools3.dms.view import RollbackException
+from obitools3.dms.view.typed_view.view_NUC_SEQS cimport View_NUC_SEQS
+from obitools3.dms.column.column cimport Column, Column_line
+from obitools3.dms.capi.obiview cimport QUALITY_COLUMN, COUNT_COLUMN, NUC_SEQUENCE_COLUMN, ID_COLUMN
+from obitools3.dms.capi.obitypes cimport OBI_INT, OBI_STR, index_t, OBI_QUAL
+from obitools3.apps.optiongroups import addMinimalOutputOption
+from obitools3.uri.decode import open_uri
+from obitools3.apps.config import logger
+from obitools3.align._qsassemble import QSolexaReverseAssemble
+from obitools3.align._qsrassemble import QSolexaRightReverseAssemble
+from obitools3.align._solexapairend import buildConsensus, buildJoinedSequence
+from obitools3.dms.obiseq cimport Nuc_Seq_Stored
+
+
+__title__="Aligns paired-ended reads"
+
+
+def addOptions(parser):
+
+    addMinimalOutputOption(parser)
+ 
+    group = parser.add_argument_group('obi alignpairedend specific options')
+ 
+    group.add_argument('-f', '--forward-reads',
+                     action="store", dest="alignpairedend:forward",
+                     metavar="<FILENAME>",
+                     default=None,
+                     type=str,
+                     help="URI to the forward reads")
+
+    group.add_argument('-r', '--reverse-reads',
+                     action="store", dest="alignpairedend:reverse",
+                     metavar="<FILENAME>",
+                     default=None,
+                     type=str,
+                     help="URI to the reverse reads")
+
+# TODO
+#    group.add_argument('--index-file',
+#                     action="store", dest="alignpairedend:index",
+#                     metavar="<FILENAME>",
+#                     default=None,
+#                     type=str,
+#                     help="URI to the index reads")
+
+    group.add_argument('--score-min',
+                     action="store", dest="alignpairedend:smin",
+                     metavar="#.###",
+                     default=None,
+                     type=float,
+                     help="Minimum score for keeping alignments")
+
+
+la = QSolexaReverseAssemble()
+ra = QSolexaRightReverseAssemble()
+def buildAlignment(Nuc_Seq_Stored direct, Nuc_Seq_Stored reverse):
+    
+    if len(direct)==0 or len(reverse)==0:
+        return None
+    
+    la.seqA = direct
+    la.seqB = reverse
+    
+    ali=la()
+    ali.direction='left'
+
+    ra.seqA = direct
+    ra.seqB = reverse
+
+    rali=ra()
+    rali.direction='right'
+     
+    if ali.score < rali.score:
+        ali = rali
+    
+    return ali
+
+
+def alignmentIterator(View_NUC_SEQS forward, View_NUC_SEQS reverse):    
+    for i in range(len(forward)):
+        ali = buildAlignment(forward[i], reverse[i])
+        if ali is None:
+            continue
+        yield ali
+
+
+def run(config):
+        
+    DMS.obi_atexit()
+    
+    logger("info", "obi alignpairedend")
+    
+    # Open the input
+    finput = open_uri(config["alignpairedend"]["forward"], config)
+    if finput is None:
+        raise Exception("Could not open forward reads")
+    forward = finput[1]
+    rinput = open_uri(config["alignpairedend"]["reverse"], config)
+    if rinput is None:
+        raise Exception("Could not open reverse reads")
+    reverse = rinput[1]
+
+    if len(forward) != len(reverse):
+        raise Exception("Error: the number of forward and reverse reads are different")
+    
+    if "index" in config["alignpairedend"]:
+        index = open_uri(config["alignpairedend"]["index"], config)
+    else:
+        index = None
+    
+    # Open the output
+    output = open_uri(config['obi']['outputURI'],
+                      config,
+                      input=False,
+                      newviewtype=View_NUC_SEQS)
+    if output is None:
+        raise Exception("Could not create output view")
+
+    view = output[1]
+    Column.new_column(view, b"QUALITY", OBI_QUAL)   #TODO output URI quality option
+
+    if 'smin' in config['alignpairedend']:
+        config['alignpairedend']['sminL'] = config['alignpairedend']['smin']
+        config['alignpairedend']['sminR'] = config['alignpairedend']['smin']
+        sminL = config['alignpairedend']['sminL']
+        sminR = config['alignpairedend']['sminR']
+    else:
+        sminL = 0
+        sminR = 0
+
+    # Initialize the progress bar
+    pb = ProgressBar(len(forward), config, seconde=5)
+    
+    ba = alignmentIterator(forward, reverse)
+            
+    i = 0
+    for ali in ba:
+        
+        pb(i)
+        
+        consensus = view[i]
+        
+        if sminL > 0:
+            if ((ali.direction=='left' and ali.score > sminL)
+                or (ali.score > sminR)):
+                buildConsensus(ali, consensus)
+            else:
+                buildJoinedSequence(ali, reverse[i], consensus)
+                
+            consensus[b"sminL"] = sminL
+            consensus[b"sminR"] = sminR
+        else:
+            buildConsensus(ali, consensus)
+
+# TODO            
+#        if "index" in config['alignpairedend'] and config['alignpairedend']['index'] is not None:
+#            idx = str(config['alignpairedend']['index'].next())
+#            consensus[b"illumina_index"] = idx
+    
+        i+=1
+    
+    finput[0].close()
+    rinput[0].close()
+    output[0].close()
+