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)
Datasets
Models
