--- a
+++ b/bin/DeepMod_scripts/myGetFeatureBasedPos.py
@@ -0,0 +1,782 @@
+
+import os;
+import sys;
+import string;
+import glob;
+import time
+import copy
+
+import h5py
+import numpy as np
+import multiprocessing
+
+from collections import defaultdict
+from distutils.version import LooseVersion
+
+import tempfile
+import subprocess
+
+import re;
+
+from . import myCom
+from . import myDetect
+
+#
+# map long reads
+# then call another function to get feature for each base of interest
+#
+def mGetFeature1(moptions, sp_options, f5files):
+   # associate signals to events
+   f5data = myDetect.get_Event_Signals(moptions, sp_options, f5files)
+
+   # save all sequences information
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG: start_time = time.time();
+   temp_fa = tempfile.NamedTemporaryFile(suffix='.fa', mode='w')
+   f5keys = sorted(f5data.keys()); #f5keys.sort()
+   for f5k in f5keys:
+      temp_fa.write(''.join(['>', f5k, '\n', f5data[f5k][0], '\n']))
+   temp_fa.flush();
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG:
+      end_time = time.time();
+      print ("Write consuming time %d" % (end_time-start_time))
+
+   # alignment using bwa-mem or minimap2
+   temp_sam = tempfile.NamedTemporaryFile()
+   if moptions['alignStr']=='bwa':
+      cmd_opt = ['mem', '-x', 'ont2d', '-v', '1', '-t', '1', moptions['Ref'], temp_fa.name]
+   else:
+      cmd_opt = ['-ax', 'map-ont', moptions['Ref'], temp_fa.name]
+   returncode = subprocess.call([moptions['alignStr'],]+cmd_opt, stdout=temp_sam)
+   if not returncode==0:
+      print ('Fatal Error!!! returncode is non-zero(%d) for "%s"' % (returncode, curcmd))
+      errkey = "Cannot running aligment"
+      for f5k in f5keys:
+         sp_options["Error"][errkey].append(f5data[f5k][3])
+      return;
+
+   temp_fa.close();
+   temp_sam.seek(0);
+   # get sam information
+   align_info = temp_sam.readlines()
+   align_info = [str(align_info[i], 'utf-8').strip() for i in range(len(align_info))]
+   temp_sam.close();
+
+   sp_param = defaultdict();
+   sp_param['f5data'] = f5data
+
+   # for alignment
+   f5align = defaultdict()
+   f5keydict = defaultdict();
+   sp_param['ref_info'] = defaultdict()
+
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG:start_time = time.time();
+   ilid = 0;
+   # for each record in sam, get alignment information
+   while ilid < len(align_info):
+      if len(align_info[ilid])==0 or align_info[ilid][0]=='@':
+         ilid += 1
+         continue;
+
+      sp_param['f5status'] = "";
+      sp_param['line'] = align_info[ilid]
+      qname = handle_line(moptions, sp_param, f5align)
+      if sp_param['f5status'] == "":
+         f5keydict[qname] = True;
+      ilid += 1
+
+   # for unmapped reads
+   for f5k in f5keys:
+      if f5k not in f5keydict:
+         sp_options["Error"]["Not in alignment sam"].append(f5data[f5k][3])
+
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG:
+      end_time = time.time();
+      print ("Get BAM consuming time %d" % (end_time-start_time))
+
+   sp_param['f5status']= ""
+   sp_param['line'] = ""
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG:start_time = time.time();
+   # handle each alignment
+   handle_record(moptions, sp_options, sp_param, f5align, f5data)
+   if moptions['outLevel']<=myCom.OUTPUT_DEBUG:
+      end_time = time.time();
+      print ("Analyze & annotate & save consuming time %d" % (end_time-start_time))
+
+#
+# get mapping information
+# then call another function to get feature of each base in a long read
+#
+def handle_record(moptions, sp_options, sp_param, f5align, f5data):
+   alignkeys = list(f5align.keys());
+   # alignment detail
+   numreg = re.compile('\d+')
+   mdireg = re.compile('[MIDNSHPX=]{1}')
+
+   feat_file_ind_dict = []
+   feat_list = None; feat_file_ind = 0
+   start_c_time = time.time();
+
+   for readk in alignkeys:
+     if len(feat_file_ind_dict)>0 and feat_list.nbytes > moptions['size_per_batch']:
+        # save features when the size is larger than the defined size
+        cur_feat_file_base = sp_options['ctfolder'] + '/'+str(feat_file_ind)
+        np.savetxt(cur_feat_file_base+'.xy.gz', feat_list, fmt='%.3f')
+        with open(cur_feat_file_base+'.xy.ind', 'w') as ind_mw:
+            for f_ind in feat_file_ind_dict:
+               ind_mw.write('%d %s\n' % (f_ind[1], f_ind[0]))
+        print ("\t%s-%d Total consuming time %d" % (sp_options['ctfolder'][sp_options['ctfolder'].rfind('/'):], feat_file_ind, time.time()-start_c_time)); sys.stdout.flush()
+        feat_file_ind_dict = []
+        feat_list = None;
+        feat_file_ind += 1
+
+     # get alignment detail
+     mapq, flag, rname, pos, cigar, readseq = f5align[readk]
+
+     if not ( (rname in moptions['fulmodlist'] and len(moptions['fulmodlist'][rname])>0) or \
+        ((not moptions['anymodlist']==None) and rname in moptions['anymodlist'] and len(moptions['anymodlist'][rname])>0) or \
+        ((not moptions['nomodlist']==None) and rname in moptions['nomodlist'] and len(moptions['nomodlist'][rname])>0) ):
+        continue;
+
+     # get reference sequece
+     if rname not in sp_param['ref_info']:
+        myDetect.getRefSeq(moptions, sp_param, rname)
+     refseq = sp_param['ref_info'][rname]
+
+     # mapped starting position and strand
+     pos = pos - 1
+     forward_reverse = '-' if flag&0x10 else '+'
+
+     numinfo = numreg.findall(cigar);
+     mdiinfo = mdireg.findall(cigar)
+     numinfo = [int(numinfo[i]) for i in range(len(numinfo))] #map(int, numinfo)
+
+     # remove clip from both tails
+     leftclip = 0; rightclip = 0;
+     while mdiinfo[0] in ['I', 'D', 'N', 'S', 'H', 'P', 'X']:
+         if mdiinfo[0] in ['I', 'S', 'X']:
+            leftclip += numinfo[0];  readseq = readseq[numinfo[0]:]
+         if mdiinfo[0] in ['H']: leftclip += numinfo[0]
+         if mdiinfo[0] in ['D', 'N', 'X']:
+            pos += numinfo[0]
+         numinfo = numinfo[1:];  mdiinfo = mdiinfo[1:]
+     while mdiinfo[-1] in ['I', 'D', 'N', 'S', 'H', 'P', 'X']:
+         if mdiinfo[-1] in ['I', 'S', 'X']:
+            rightclip += numinfo[-1]; readseq = readseq[:-numinfo[-1]]
+         if mdiinfo[-1] in ['H']: rightclip += numinfo[-1]
+         numinfo = numinfo[:-1]; mdiinfo = mdiinfo[:-1]
+     if forward_reverse=='+':
+        if rightclip>0: m_event = f5data[readk][1][leftclip:-rightclip]
+        else: m_event = f5data[readk][1][leftclip:]
+     else:
+        if leftclip>0: m_event = f5data[readk][1][rightclip:-leftclip]
+        else: m_event = f5data[readk][1][rightclip:]
+
+     # is in region of interest if provided
+     isinreg = False;
+     if (moptions['region'][0] in ['', None, rname]) and \
+        (moptions['region'][1] in ['', None] or pos>moptions['region'][1]) and \
+        (moptions['region'][2] in ['', None] or pos+len(m_event)<moptions['region'][2]):
+        isinreg = True;
+     if not isinreg:
+        continue;
+
+     # associate mapped reference positions with read positions
+     lastmatch = None; firstmatch = None;
+     first_match_pos = None; last_match_pos = None
+     last_al_match = None;  first_al_match = None
+     lasmtind = 0;
+     base_map_info = []; #indel_groups = defaultdict()
+     nummismatch = 0; numinsert = 0; numdel = 0;
+     read_ind = 0;
+     for n1ind in range(len(numinfo)):
+        mdi = mdiinfo[n1ind];
+        # for each mapped types
+        for n1i in range(numinfo[n1ind]):
+           if mdi=='M':
+              base_map_info.append((refseq[pos], readseq[read_ind], pos, read_ind))
+              if refseq[pos]==readseq[read_ind]:
+                 if firstmatch==None: firstmatch = read_ind
+                 if lastmatch==None or lastmatch<read_ind: lastmatch = read_ind; lasmtind=n1ind
+                 if first_al_match==None: first_al_match=len(base_map_info)-1
+                 if last_al_match==None or last_al_match<len(base_map_info): last_al_match=len(base_map_info)-1
+                 if first_match_pos==None: first_match_pos = pos
+                 if last_match_pos==None or last_match_pos<pos: last_match_pos = pos
+              else: nummismatch += 1
+              pos += 1; read_ind += 1;
+           elif mdi =='I':
+              base_map_info.append(('-', readseq[read_ind], pos, read_ind))
+              read_ind += 1;
+              numinsert += 1
+           elif mdi == 'D':
+              base_map_info.append((refseq[pos], '-', pos, read_ind))
+              pos += 1;
+              numdel += 1
+           elif mdi == 'N':
+              base_map_info.append((refseq[pos], '-', pos, read_ind))
+              pos += 1;
+              if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+                 print ('CIGAR-Error N exist', f5data[readk][3])
+           elif mdi == 'S':
+              read_ind += 1;
+              if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+                 print ('CIGAR-Error!!! S in the middle of the sequence', f5data[readk][3])
+           elif mdi == 'H':
+              if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+                 print ('CIGAR-Error!!! H in the middle of the sequence', f5data[readk][3])
+           elif mdi == 'P':
+              if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+                 print ('CIGAR-Error!!! P exist', f5data[readk][3])
+           elif mdi == '=':
+             base_map_info.append((refseq[pos], readseq[read_ind], pos, read_ind))
+             if first_match_pos==None: first_match_pos  = pos
+             if last_match_pos==None or last_match_pos<pos: last_match_pos = pos
+             pos += 1; read_ind += 1;
+             if firstmatch==None: firstmatch = read_ind - 1
+             if lastmatch==None or lastmatch<read_ind-1: lastmatch = read_ind - 1; lasmtind=n1ind
+             if last_al_match==None or last_al_match<len(base_map_info): last_al_match=len(base_map_info)-1
+             if first_al_match==None: first_al_match=len(base_map_info)-1
+           elif mdi == 'X':
+             base_map_info.append((refseq[pos], readseq[read_ind], pos, read_ind))
+             pos += 1; read_ind += 1;
+             nummismatch += 1
+           else:
+             if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+                print ('CIGAR-Error!!!', 'Warning unknow CIGAR element ' + str(numinfo[n1ind]) + ' ' + mdi, f5data[readk][3])
+     if firstmatch==None or lastmatch==None or firstmatch<0 or lastmatch<0:
+        if moptions['outLevel']<=myCom.OUTPUT_WARNING:
+           print ("Errorfast5 "+f5data[readk][3])
+           print('match-Error!!! no first and/or last match',f5data[readk][3],('firstmatch=%d' % firstmatch) if not (firstmatch==None) else "N", ('lastmatch%d' % lastmatch) if not (lastmatch==None) else "N", str(flag), rname, str(pos));
+           print('\tf=%d, chr=%s, p=%d, c=%s, s=%s' % (flag, rname, pos, cigar, readseq))
+           continue;
+
+     # remove unmatch in both tails
+     if not firstmatch==None: leftclip += firstmatch
+     if (not lastmatch==None) and len(m_event)-lastmatch>1: rightclip += len(m_event)-lastmatch-1
+     # remove events whose bases are not mapped.
+     if forward_reverse=='+':
+        if len(m_event)-lastmatch>1:
+           m_event = m_event[firstmatch:(lastmatch+1-len(m_event))]
+        elif firstmatch>0: m_event = m_event[firstmatch:]
+     else:
+        if firstmatch>0: m_event = m_event[(len(m_event)-1-lastmatch):-firstmatch]
+        elif len(m_event)-lastmatch>1: m_event = m_event[(len(m_event)-1-lastmatch):]
+     # print detail if unexpected error occurs
+     if firstmatch>0 or len(base_map_info)-last_al_match>1:
+        if moptions['outLevel']<=myCom.OUTPUT_WARNING and ((firstmatch>0) or (len(base_map_info)-last_al_match>1 and refseq[last_match_pos+1] not in ['N'])):
+           print ("Errorfast5"+f5data[readk][3])
+           print ('Warning!!! first not match', firstmatch, lastmatch, first_al_match, last_al_match, len(base_map_info), numinfo[lasmtind-2:(lasmtind+5)], mdiinfo[lasmtind-2:(lasmtind+5)], lasmtind, len(numinfo))
+           print('\tref='+refseq[last_match_pos:last_match_pos+20]+"\n\tred="+readseq[lastmatch:lastmatch+20])
+           if firstmatch>0:
+              print('\tref='+refseq[(first_match_pos-20 if first_match_pos-20>0 else 0):first_match_pos]+"\n\tred="+readseq[(firstmatch-20 if firstmatch-20>0 else 0):firstmatch])
+           print('\tf=%d, chr=%s, p=%d, c=%s, s=%s' % (flag, rname, pos, cigar, readseq)) # flag, rname, pos, cigar, readseq
+
+        if len(base_map_info)-last_al_match>1:
+           base_map_info = base_map_info[first_al_match:(last_al_match+1-len(base_map_info))]
+        elif first_al_match>0:
+           base_map_info = base_map_info[first_al_match:]
+
+     # post-process mapping information
+     base_map_info = np.array(base_map_info, dtype=[('refbase', 'U1'), ('readbase', 'U1'), ('refbasei', np.uint64), ('readbasei', np.uint64)])
+     if forward_reverse=='-':
+        base_map_info = np.flipud(base_map_info)
+        for bmii in range(len(base_map_info)):
+            base_map_info['refbase'][bmii]  = get_complement(base_map_info['refbase'][bmii])
+            base_map_info['readbase'][bmii] = get_complement(base_map_info['readbase'][bmii])
+        leftclip, rightclip = rightclip, leftclip
+     if False: #True: # for test base_map_info  ### for check consistency
+        ref_align_key = '/Analyses/NanomoCorrected_000/BaseCalled_template/Alignment/genome_alignment'
+        read_align_key = '/Analyses/NanomoCorrected_000/BaseCalled_template/Alignment/read_alignment'
+        with h5py.File(f5data[readk][3], 'r') as mf5:
+           read_align_list = [bt.decode(encoding="utf-8") for bt in mf5[read_align_key]]
+           ref_align_list = [bt.decode(encoding="utf-8") for bt in mf5[ref_align_key]]
+           for rali in range(len(read_align_list)):
+              if not read_align_list[rali]==base_map_info['readbase'][rali]:
+                 print ("Error not equal1! %s %s %d %s" % (read_align_list[rali], base_map_info['readbase'][rali], rali, f5data[readk][3]))
+              if not ref_align_list[rali]==base_map_info['refbase'][rali]:
+                 print ("Error not equal2! %s %s %d %s" % (ref_align_list[rali], base_map_info['refbase'][rali], rali, f5data[readk][3]))
+     #
+     # handle map like
+     # CCG    or CGG
+     # C-G       C-G
+     #
+     if 'motif' in moptions and moptions['motif'][0]=='CG':
+        for ali in range(len(base_map_info)):
+           if base_map_info['refbase'][ali]=='C' and base_map_info['readbase'][ali]=='C':
+              if ali+1<len(base_map_info) and base_map_info['readbase'][ali+1]=='-' and base_map_info['refbase'][ali+1]=='G':
+                 addali = 2;
+                 while ali + addali < len(base_map_info):
+                     if base_map_info['readbase'][ali+addali]=='-' and base_map_info['refbase'][ali+addali]=='G': addali += 1;
+                     else: break;
+                 if ali + addali < len(base_map_info) and base_map_info['readbase'][ali+addali]=='G' and base_map_info['refbase'][ali+addali]=='G':
+                    base_map_info['readbase'][ali+1], base_map_info['readbase'][ali+addali] = base_map_info['readbase'][ali+addali], base_map_info['readbase'][ali+1]
+           if base_map_info['refbase'][ali]=='G' and base_map_info['readbase'][ali]=='G':
+              if ali-1>-1 and base_map_info['readbase'][ali-1]=='-' and base_map_info['refbase'][ali-1]=='C':
+                 addali = 2;
+                 while ali - addali >-1:
+                     if base_map_info['readbase'][ali-addali]=='-' and base_map_info['refbase'][ali-addali]=='C': addali += 1;
+                     else: break;
+                 if ali - addali>-1 and base_map_info['readbase'][ali-addali]=='C' and base_map_info['refbase'][ali-addali]=='C':
+                     base_map_info['readbase'][ali-1], base_map_info['readbase'][ali-addali] = base_map_info['readbase'][ali-addali], base_map_info['readbase'][ali-1]
+     # too short reads
+     if len(m_event)<500:
+         sp_options["Error"]["Less(<500) events"].append(f5data[readk][3])
+         continue;
+
+     # get feautre
+     mfeatures,isdif = get_Feature(moptions, sp_options, sp_param, f5align, f5data, readk, leftclip, rightclip, base_map_info, forward_reverse, rname, first_match_pos, numinsert, numdel)
+     if isdif and moptions['outLevel']<=myCom.OUTPUT_WARNING:
+        print("Dif is true")
+        print([lastmatch, firstmatch, first_match_pos, last_match_pos, first_al_match, last_al_match, lasmtind, len(base_map_info), nummismatch, numinsert, numdel, len(base_map_info)-nummismatch-numinsert-numdel])
+
+     # merge to previously handled features of other fast5 files
+     if len(mfeatures)>0:
+        if len(feat_file_ind_dict)==0:
+           feat_file_ind_dict.append((f5data[readk][3], 0));
+           feat_list = mfeatures
+        else:
+           feat_file_ind_dict.append((f5data[readk][3], len(feat_list)))
+           feat_list = np.concatenate((feat_list, mfeatures), axis=0)
+
+   # store the last feature data.
+   if len(feat_file_ind_dict)>0:
+      cur_feat_file_base = sp_options['ctfolder'] + '/'+str(feat_file_ind)
+      np.savetxt(cur_feat_file_base+'.xy.gz', feat_list, fmt='%.3f')
+      with open(cur_feat_file_base+'.xy.ind', 'w') as ind_mw:
+          for f_ind in feat_file_ind_dict:
+             ind_mw.write('%d %s\n' % (f_ind[1], f_ind[0]))
+      print ("\t%s-%d Total consuming time %d" % (sp_options['ctfolder'][sp_options['ctfolder'].rfind('/'):], feat_file_ind, time.time()-start_c_time)); sys.stdout.flush()
+      feat_file_ind_dict = []
+      feat_list = None;
+      feat_file_ind += 1
+
+#
+# get feature for each base of interest in long reads according to raw signals and mapping information
+#
+def get_Feature(moptions, sp_options, sp_param, f5align, f5data, readk, start_clip, end_clip, base_map_info, forward_reverse, rname, mapped_start_pos, num_insertions, num_deletions):
+   # event information
+   modevents = sp_param['f5data'][readk][1]
+   # class number, bin num and bin length
+   clnum = 2; binnum = 50; binlen = 0.2;
+   if forward_reverse=='+':
+      align_ref_pos = mapped_start_pos
+   else:
+      align_ref_pos = mapped_start_pos + len(base_map_info) - num_insertions - 1
+
+   # initialize feature matrix for all events.
+   if moptions['fnum']==57:
+      #mfeatures = np.zeros((len(modevents)-end_clip+100-(start_clip-100), (binnum+3+3+4)));
+      mfeatures = np.zeros((len(modevents)-end_clip+100-(start_clip-100), (binnum+3+3+4)));
+   else: mfeatures = np.zeros((len(modevents)-end_clip+100-(start_clip-100), (3+3+4)));
+
+   # filter poor alignment
+   checkneighbornums = [3,6]
+   checkratios = {3:[6,5,4,2], 6:[11,10,9,3]}
+   checkratios = {3:[6,5,4,2], 6:[12,10,9,3]}
+   cgpos = [[], []]
+   affectneighbor = 1; # 2;
+   for aligni in range(len(base_map_info)):
+      # for methylated positions and not-used adjacent positions
+      if 'motif' in moptions and base_map_info['readbase'][aligni]==moptions['motif'][0][moptions['motif'][1]]:
+         m_a_st = aligni-moptions['motif'][1]; m_a_end = aligni+len(moptions['motif'][0])-moptions['motif'][1]
+         if m_a_st>-1 and m_a_end<=len(base_map_info) and ''.join(base_map_info['readbase'][m_a_st:m_a_end])==moptions['motif'][0] and (not ''.join(base_map_info['refbase'][m_a_st:m_a_end])==moptions['motif'][0]):
+            cgpos[1].extend([(forward_reverse, base_map_info['refbasei'][addi]) for addi in range(aligni-affectneighbor if aligni-affectneighbor>-1 else 0, aligni+affectneighbor+1 if aligni+affectneighbor+1<len(base_map_info) else len(base_map_info))])
+      if (not base_map_info['refbase'][aligni]=='-') and \
+         (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['fulmodlist'][rname]:
+         if not base_map_info['readbase'][aligni]=='-':
+            nextnogap = aligni + 1;
+            while nextnogap<len(base_map_info):
+               if not base_map_info['refbase'][nextnogap]=='-': break;
+               nextnogap += 1
+            iscg = False;
+            # find gaps
+            for checkneighbornum in checkneighbornums:
+               if not nextnogap<len(base_map_info): continue;
+               matchnum = 0; gapnum = 0;
+               # get gaps for two window sizes
+               for checki in range(aligni-checkneighbornum, aligni+checkneighbornum+1):
+                  if checki>-1 and checki<len(base_map_info):
+                     if base_map_info['refbase'][checki]==base_map_info['readbase'][checki]: matchnum += 1
+                     if base_map_info['refbase'][checki]=='-' or base_map_info['readbase'][checki]=='-': gapnum += 1
+               if gapnum<=checkratios[checkneighbornum][3]:
+                  for addi in range(aligni-affectneighbor if aligni-affectneighbor>-1 else 0, nextnogap+affectneighbor if nextnogap+affectneighbor<len(base_map_info) else len(base_map_info)):
+                     if addi==aligni: # for methylated positions
+                        cgpos[0].append((forward_reverse, base_map_info['refbasei'][addi]))
+                     else: # for non-used positions
+                        cgpos[1].append((forward_reverse, base_map_info['refbasei'][addi]))
+                  iscg = True; break;
+            if iscg: continue;
+         # add more not-used positions if more gaps exist
+         if not base_map_info['readbase'][aligni]=='-':
+            nextnogap = aligni
+            for _ in range(affectneighbor):
+               nextnogap += 1;
+               while nextnogap<len(base_map_info['refbase']):
+                 if not base_map_info['refbase'][nextnogap]=='-': break;
+                 nextnogap += 1
+            prenogap = aligni
+            for _ in range(affectneighbor):
+               prenogap -= 1;
+               while prenogap>-1:
+                  if not base_map_info['refbase'][prenogap]=='-': break;
+                  prenogap -= 1
+
+            read0 = aligni; read1 = aligni
+            for _ in range(affectneighbor):
+               read0 -= 1
+               while read0>-1:
+                  if base_map_info['readbase'][read0]=='-': read0 -= 1
+                  else: break;
+               read1 += 1
+               while read1<len(base_map_info['readbase']):
+                  if base_map_info['readbase'][read1]=='-': read1 += 1
+                  else: break;
+
+            if read0<prenogap:
+               if read0>-1: prenogap = read0
+               else: prenogap = 0
+            if read1>nextnogap:
+               if read1<len(base_map_info['readbase']): nextnogap = read1
+               else: nextnogap = len(base_map_info['readbase'])-1
+            if prenogap<0: prenogap = 0
+            if not nextnogap<len(base_map_info['readbase']): nextnogap=len(base_map_info['readbase'])-1
+            if not prenogap<len(base_map_info['readbase']): prenogap=len(base_map_info['readbase'])-1
+            for excldi in range(prenogap, nextnogap+1):
+               cgpos[1].append((forward_reverse, base_map_info['refbasei'][excldi]))
+
+   print ('%s%s %d, %d >> %d %d, %d-%d=%d' % (forward_reverse, f5data[readk][3], len(cgpos[0]), len(cgpos[1]), len(modevents)-end_clip-start_clip, start_clip, len(modevents), end_clip, len(modevents)-end_clip))
+
+   aligni = 0; isdif = False;
+   for ie in range(start_clip-100, len(modevents)-end_clip+100):
+      cur_row_num = ie - (start_clip-100); cur_base = ''
+      # for aligned bases
+      if ie>=start_clip and ie<len(modevents)-end_clip:
+         if align_ref_pos<mapped_start_pos:
+            print ('ERRRR align_ref_pos(%d)<mapped_start_pos(%d)' % (align_ref_pos, mapped_start_pos))
+         while base_map_info['readbase'][aligni]=='-':
+            if not align_ref_pos==base_map_info['refbasei'][aligni]:
+               print ('ERRRR align_ref_pos(%d) not equal to %d' % (align_ref_pos, base_map_info['refbasei'][aligni] ))
+            if not base_map_info['refbase'][aligni]=='-':
+               if forward_reverse=='+': align_ref_pos += 1
+               else: align_ref_pos -= 1
+            aligni += 1
+         if not base_map_info['readbase'][aligni] == modevents['model_state'][ie][2]:
+            print ('Error Does not match', base_map_info['readbase'][aligni], modevents['model_state'][ie][2], aligni, ie)
+            isdif = True;
+         # the first column is the aligned reference position
+         mfeatures[cur_row_num][0] = align_ref_pos
+         cur_base = base_map_info['refbase'][aligni]
+         # the second/third column is for negative/positive labels of methylation
+         if moptions['posneg'] == 0: # for a data without any modification
+            if ( (not moptions['anymodlist']==None) and rname in moptions['nomodlist'] and (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['nomodlist'][rname] ):
+                mfeatures[cur_row_num][1] = 1; mfeatures[cur_row_num][2] = 0
+            elif (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['fulmodlist'][rname]:
+                mfeatures[cur_row_num][1] = 1; mfeatures[cur_row_num][2] = 0
+            elif ((not moptions['anymodlist']==None) and rname in moptions['anymodlist'] and (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['anymodlist'][rname] ):
+                mfeatures[cur_row_num][1] = 1; mfeatures[cur_row_num][2] = 0
+         else: # for a data with both modified and un-modified positions
+            if (forward_reverse, base_map_info['refbasei'][aligni]) in cgpos[0] and (not base_map_info['refbase'][aligni]=='-'):
+               mfeatures[cur_row_num][1] = 0; mfeatures[cur_row_num][2] = 1
+            else:
+               if (forward_reverse, base_map_info['refbasei'][aligni]) not in cgpos[1]:
+                  if moptions['anymodlist']==None:
+                      if moptions['nomodlist']==None or ( rname in moptions['nomodlist'] and (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['nomodlist'][rname] ):
+                         mfeatures[cur_row_num][1] = 1; mfeatures[cur_row_num][2] = 0
+                  elif rname in moptions['anymodlist'] and (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['anymodlist'][rname]:
+                      pass
+                  else:
+                      if moptions['nomodlist']==None or ( rname in moptions['nomodlist'] and (forward_reverse, base_map_info['refbasei'][aligni]) in moptions['nomodlist'][rname] ):
+                         mfeatures[cur_row_num][1] = 1; mfeatures[cur_row_num][2] = 0
+         if not base_map_info['refbase'][aligni]=='-':
+            if forward_reverse=='+': align_ref_pos += 1
+            else: align_ref_pos -= 1
+         aligni += 1
+
+      # for bin features
+      if ie>=0 and ie<len(modevents) and moptions['fnum']==57:
+         for currs in sp_param['f5data'][readk][2][modevents['start'][ie]:int(modevents['start'][ie]+int(modevents['length'][ie]+0.5))]:
+             if currs>10 or currs<-10: print ('Error raw signal', currs, ie, modevents['start'][ie], modevents['length'][ie])
+             curbin = int((currs+5)/binlen)
+             if curbin<0: curbin = 0
+             elif not curbin<binnum: curbin = binnum-1
+             mfeatures[cur_row_num][curbin+3] += 1
+      if ie>=0 and ie<len(modevents):
+         # for reference base type feature
+         if cur_base in myCom.g_ACGT:
+            mfeatures[cur_row_num][moptions['fnum']-3+3-4+myCom.g_ACGT.index(cur_base)] = 1
+         cur_index_add = moptions['fnum'] - 3 + 3
+         # for signal mean std and length.
+         mfeatures[cur_row_num][cur_index_add + 0] = modevents["mean"][ie]
+         mfeatures[cur_row_num][cur_index_add + 1] = modevents["stdv"][ie]
+         mfeatures[cur_row_num][cur_index_add + 2] = modevents["length"][ie]
+
+   # truncated too much not-used positions
+   nbkeys = defaultdict();
+   for mfi in range(len(mfeatures)):
+      if mfeatures[mfi][1] + mfeatures[mfi][2] > 0.9:
+         for ini in range(mfi-25, mfi+26):
+            if ini<0 or ini>len(mfeatures)-1:
+               print("Warning wrong del mfeatures id %d for %s" % (ini, f5data[readk][3]))
+            else:
+               nbkeys[ini] = True;
+   keepInd = sorted(list(nbkeys.keys()));
+   if len(keepInd)>0:
+      if not len(keepInd)>len(mfeatures)*0.9:
+         mfeatures = mfeatures[np.array(keepInd)]
+   else:
+      mfeatures = []
+
+   return (mfeatures, isdif)
+
+
+#
+# get the complementary bases
+#
+def get_complement(na):
+   if na in myCom.acgt: return myCom.na_bp[na]
+   else: return na;
+
+#
+# get required information for reach mapping records.
+#
+def handle_line(moptions, sp_param, f5align):
+   lsp = sp_param['line'].split('\t')
+   qname, flag, rname, pos, mapq, cigar, _, _, _, seq, _ = lsp[:11]
+   # checked query name
+   if qname=='*': sp_param['f5status'] = "qname is *"
+   # check mapping quality
+   elif int(mapq)==255: sp_param['f5status'] = "mapq is 255"
+   # check mapped positions
+   elif int(pos)==0: sp_param['f5status'] = "pos is 0"
+   # check mapped string
+   elif cigar=='*': sp_param['f5status'] = "cigar is *"
+   # check reference name
+   elif rname=='*': sp_param['f5status'] = "rname is *"
+   if not sp_param['f5status']=="": return qname
+
+   if (qname not in f5align) or f5align[qname][0]<int(mapq):
+      f5align[qname] = (int(mapq), int(flag), rname, int(pos), cigar, seq)
+
+   return qname
+
+#
+# feature handler/workder for multiprocessing
+#
+def getFeature_handler(moptions, h5files_Q, failed_Q, version_Q):
+   while not h5files_Q.empty():
+      try:
+         # get a list of files
+         f5files, ctfolderid = h5files_Q.get(block=False)
+      except:
+         break;
+
+      sp_options = defaultdict();
+      sp_options['ctfolder'] = moptions['outFolder']+str(ctfolderid)
+      if not os.path.isdir(sp_options['ctfolder']):
+         os.system('mkdir '+sp_options['ctfolder'])
+      # get features
+      mGetFeature1(moptions, sp_options, f5files)
+      # output errors
+      for errtype, errfiles in sp_options["Error"].items():
+         failed_Q.put((errtype, errfiles));
+      # double check albacore version
+      for vk in sp_options["get_albacore_version"]:
+         version_Q.put((vk, sp_options["get_albacore_version"][vk]))
+
+#
+# read sequence information from a reference genome
+#
+def readFA(mfa, t_chr=None):
+   fadict = defaultdict();
+   with open(mfa, 'r') as mr:
+      cur_chr = None;
+      line = mr.readline();
+      while line:
+         # remove empty spaces
+         line = line.strip();
+         if len(line)>0:
+            if line[0]=='>': # for each chromosome line
+               if (not cur_chr==None) and (t_chr in [None, cur_chr]):
+                  fadict[cur_chr] = ''.join(fadict[cur_chr])
+               cur_chr = line[1:].split()[0]
+               if t_chr in [None, cur_chr]:
+                  fadict[cur_chr] = []
+            else: # for sub-sequence line in a reference file
+               if t_chr in [None, cur_chr]:
+                  fadict[cur_chr].append(line.upper())
+         line = mr.readline();
+      # for the last chromosome in the file
+      if (not cur_chr==None) and (t_chr in [None, cur_chr]):
+         fadict[cur_chr] = ''.join(fadict[cur_chr])
+   return fadict
+
+#
+# get reference positions for motif-based modifications
+#
+def readMotifMod(fadict, mpat='Cg', mposinpat=0, t_chr=None, t_start=None, t_end=None):
+   pos_dict = defaultdict(int)
+
+   # get motif and complementary motif
+   pat3 = copy.deepcopy(mpat.upper())
+   comp_pat3 = ''.join([get_complement(curna) for curna in pat3][::-1])
+   comp_mposinpat = len(comp_pat3)-1-mposinpat
+
+   fakeys = fadict.keys();
+   cpgdict = defaultdict(int);
+   all_a = defaultdict()
+   for fak in fakeys:
+       cpgnum = [0, 0]
+       # motif-based reference positions
+       cpgdict[fak] = defaultdict()
+       # position of bases of interest
+       all_a[fak] = defaultdict()
+       for i in range(len(fadict[fak])):
+          if (t_start==None or i>=t_start) and (t_end==None or i<=t_end):
+             if fadict[fak][i]==mpat[mposinpat]: # for forward strand
+                all_a[fak][('+', i)] = True;
+             elif get_complement(fadict[fak][i])==mpat[mposinpat]: # for reverse strand
+                all_a[fak][('-', i)] = True;
+
+             # check motif in forward strand
+             if i-mposinpat>=0 and i+len(comp_pat3)-1-mposinpat<len(fadict[fak]) and ''.join(fadict[fak][i-mposinpat:(i+len(comp_pat3)-1-mposinpat+1)])==pat3:
+                cpgdict[fak][('+', i)] = [1, fadict[fak][i]]; cpgnum[0] += 1
+             elif i-comp_mposinpat>=0 and i+len(comp_pat3)-1-comp_mposinpat<len(fadict[fak]) and ''.join(fadict[fak][i-comp_mposinpat:(i+len(comp_pat3)-1-comp_mposinpat+1)])==comp_pat3: # check motif in reverse strand
+                cpgdict[fak][('-', i)] = [1, fadict[fak][i]]; cpgnum[1] += 1
+             else:
+                pass
+       print('%s%d site: %d(+) %d(-) for %s' % (pat3, mposinpat, cpgnum[0], cpgnum[1], fak))
+   return (cpgdict, all_a)
+
+
+#
+# a multiprocessing manager to get features from all long reads.
+#
+def getFeature_manager(moptions):
+   start_time = time.time();
+   # multipprocessing manager
+   pmanager = multiprocessing.Manager();
+
+   # prepare output folder
+   if os.path.isdir(moptions['outFolder']):
+      os.system('rm -dr '+moptions['outFolder'])
+   if not os.path.isdir(moptions['outFolder']):
+      os.system('mkdir '+moptions['outFolder'])
+
+   moptions['size_per_batch'] = moptions['size_per_batch'] * (10**7)
+
+   # read reference information
+   fadict = readFA(moptions['Ref'],moptions['region'][0])
+   if moptions['motifORPos']==1: # get motif-based positions for modifications
+      moptions['fulmodlist'], moptions['nomodlist'] = readMotifMod(fadict, moptions['motif'][0], moptions['motif'][1], moptions['region'][0], moptions['region'][1], moptions['region'][2])
+      moptions['anymodlist'] = None
+      moptions['nomodlist'] = None; # add for simple process
+   elif moptions['motifORPos']==2: # modification position is specified by the files
+      fuldfiles = glob.glob(moptions["fulmod"]);
+      moptions['fulmodlist'] = defaultdict(lambda: defaultdict());
+      if not moptions["anymod"]==None: # partially modified positions
+         anydfiles = glob.glob(moptions["anymod"])
+         moptions['anymodlist'] = defaultdict(lambda: defaultdict());
+      else:
+         moptions['anymodlist'] = None
+      if not moptions["nomod"]==None: # completely un-modified positions
+         nodfiles = glob.glob(moptions["nomod"])
+         moptions['nomodlist'] = defaultdict(lambda: defaultdict());
+      else:
+         moptions['nomodlist'] = None
+      mthreadin = [moptions['fulmodlist'], moptions['anymodlist'], moptions['nomodlist']]
+      mthfiles = [fuldfiles, anydfiles, nodfiles]
+      # read completely modified positions, partially modified positions, completely un-modified positions from files
+      for mthi in range(len(mthreadin)):
+         curmeth = mthreadin[mthi]; curfilelist = mthfiles[mthi]
+         if curmeth==None or curfilelist==None: continue;
+         for curmthf in curfilelist:
+             with open(curmthf, 'r') as mreader:
+                line = mreader.readline();
+                while line:
+                   if len(line)>0:
+                      tchr, tstrand, tpos = line.split()[:3]
+                      curmeth[tchr][(tstrand, int(tpos))] = [1-mthi, fadict[tchr][int(tpos)]];
+                   line = mreader.readline();
+   for tchr in moptions['fulmodlist'] if moptions['anymodlist']==None else moptions['anymodlist']:
+      if len(moptions['fulmodlist'][tchr])>0 or ((not moptions['anymodlist']==None) and len(moptions['anymodlist'][tchr])>0):
+          print ('%s fulmod=%d anymod=%d nomod=%d' % (tchr, len(moptions['fulmodlist'][tchr]), len(moptions['anymodlist'][tchr]) if (not moptions['anymodlist']==None) else -1, len(moptions['nomodlist'][tchr]) if (not moptions['nomodlist']==None) else -1))
+
+   if True: #False:
+      # get all input fast5 files
+      f5files = glob.glob(os.path.join(moptions['wrkBase'],"*.fast5" ))
+      if moptions['recursive']==1:
+         f5files.extend(glob.glob(os.path.join(moptions['wrkBase'],"*/*.fast5" )))
+         f5files.extend(glob.glob(os.path.join(moptions['wrkBase'],"*/*/*.fast5" )))
+         f5files.extend(glob.glob(os.path.join(moptions['wrkBase'],"*/*/*/*.fast5" )))
+
+
+   print('Total files=%d' % len(f5files))
+   h5files_Q = pmanager.Queue();
+   failed_Q = pmanager.Queue()
+   version_Q = pmanager.Queue()
+
+   # split input fast5 files into different batch
+   h5_batch = []; h5batchind = 0;
+   for f5f in f5files:
+      h5_batch.append(f5f);
+      if len(h5_batch)==moptions['files_per_thread']:
+         h5files_Q.put((h5_batch, h5batchind))
+         h5batchind += 1
+         h5_batch = []; #break; ### feature500
+   if len(h5_batch)>0:
+      h5files_Q.put((h5_batch, h5batchind))
+
+   # each thread handle a batch a time and repeat for all batches.
+   share_var = (moptions, h5files_Q, failed_Q, version_Q)
+   handlers = []
+   for hid in range(moptions['threads']):
+      p = multiprocessing.Process(target=getFeature_handler, args=share_var);
+      p.start();
+      handlers.append(p);
+
+   # get failed files.
+   failed_files = defaultdict(list);
+   version_default = defaultdict(lambda: defaultdict(int));
+   while any(p.is_alive() for p in handlers):
+      try:
+         errk, fns = failed_Q.get(block=False);
+         failed_files[errk].extend(fns)
+         curv, curv_num = version_Q.get(block=False);
+         version_default[curv] += curv_num
+      except:
+         time.sleep(1);
+         continue;
+
+   # output failure information
+   if len(failed_files)>0:
+      print ('Error information for different fast5 files:')
+      for errtype, errfiles in failed_files.items():
+         print ('\t%s %d' % (errtype, len(errfiles)))
+   print("abversion info {}".format(str(version_default)))
+   sys.stdout.flush()
+   end_time = time.time();
+   print ("Total consuming time %d" % (end_time-start_time))
+
+
+
+# for indepdent testing of code
+if __name__=='__main__':
+#   if len(sys.argv)>4:
+      moptions = {}
+      moptions['basecall_1d'] = 'Basecall_1D_000'
+      moptions['basecall_1d'] = ['Basecall_1D_000']
+      moptions['basecall_2strand'] = 'BaseCalled_template'
+
+      moptions['outLevel'] = myCom.OUTPUT_WARNING
+      moptions['outLevel'] = myCom.OUTPUT_INFO
+
+      moptions['modfile'] = '../../mod_output/train1/2/mod_train'
+
+      moptions['fnum'] = 53;
+      moptions['hidden'] = 100;
+      moptions['windowsize'] = 21;
+
+      moptions['threads'] = 8
+      moptions['threads'] = 1
+      moptions['files_per_thread'] = 500
+
+      mDetect_manager(moptions)