#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="", default=None, type=str, help="URI to the forward reads") group.add_argument('-r', '--reverse-reads', action="store", dest="alignpairedend:reverse", metavar="", default=None, type=str, help="URI to the reverse reads") # TODO # group.add_argument('--index-file', # action="store", dest="alignpairedend:index", # metavar="", # 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"]) if finput is None: raise Exception("Could not open forward reads") forward = finput[1] rinput = open_uri(config["alignpairedend"]["reverse"]) 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"]: # TODO index = open_uri(config["alignpairedend"]["index"]) else: index = None # Open the output output = open_uri(config['obi']['outputURI'], 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()