diff --git a/misc/training/README.md b/misc/training/README.md
new file mode 100644
index 0000000..10efa7c
--- /dev/null
+++ b/misc/training/README.md
@@ -0,0 +1,5 @@
+Use generate_SNP_pileups.py to generate training features for SNP calling, and then run model_run.py to training a model.
+Use generate_indel_pileups.py or generate_indel_pileups_hifi.py (for PacBio HiFi)  to generate training features for indel calling, and then run model_run_indels.py to training a model.
+
+You can print help for each script using `--help`, e.g. `python generate_SNP_pileups.py --help`
+This code uses tensorflow 1.13 for training a model.
diff --git a/misc/training/generate_SNP_pileups.py b/misc/training/generate_SNP_pileups.py
new file mode 100644
index 0000000..28e2b6a
--- /dev/null
+++ b/misc/training/generate_SNP_pileups.py
@@ -0,0 +1,480 @@
+import sys, pysam, time, os, copy, argparse, random, datetime
+from collections import Counter
+import numpy as np
+import multiprocessing as mp
+from pysam import VariantFile
+from intervaltree import Interval, IntervalTree
+
+base_to_num_map={'*':4,'A':0,'G':1,'T':2,'C':3,'N':4}
+
+def in_bed(tree,pos):            
+        return tree.overlaps(pos)
+    
+def get_cnd_pos(v_pos,cnd_pos, seq='ont'):
+    if seq=='ont':
+        ls=cnd_pos[abs(cnd_pos-v_pos)<50000] 
+
+        ls1_0= [p for p in ls if (p>=v_pos-2000) &  (p<v_pos)][:2]
+        ls1_1= [p for p in ls if (p>=v_pos-5000) &  (p<v_pos-2000)][-3:]
+        ls1_2= [p for p in ls if (p>=v_pos-10000) & (p<v_pos-5000)][-4:]
+        ls1_3= [p for p in ls if (p>=v_pos-20000) & (p<v_pos-10000)][-5:]
+        ls1_4= [p for p in ls if                  (p<v_pos-20000)][-6:]
+
+        ls2_0= [p for p in ls if (p>v_pos) & (p<=v_pos+2000)][-2:]
+        ls2_1= [p for p in ls if (p>v_pos+2000) & (p<=v_pos+5000)][:3]
+        ls2_2= [p for p in ls if (p>v_pos+5000) & (p<=v_pos+10000)][:4]
+        ls2_3= [p for p in ls if (p>v_pos+10000) & (p<=v_pos+20000)][:5]
+        ls2_4= [p for p in ls if (p>v_pos+20000)][:6]
+
+        ls_total_1=sorted(ls1_0+ls1_1+ls1_2+ls1_3+ls1_4)
+        ls_total_2=sorted(ls2_0+ls2_1+ls2_2+ls2_3+ls2_4)
+
+    elif seq=='ul_ont':
+        ls=cnd_pos[abs(cnd_pos-v_pos)<100000] 
+
+        ls1_0= [p for p in ls if (p>=v_pos-2000) &  (p<v_pos)][:2]
+        ls1_1= [p for p in ls if (p>=v_pos-5000) &  (p<v_pos-2000)][-2:]
+        ls1_2= [p for p in ls if (p>=v_pos-10000) & (p<v_pos-5000)][-3:]
+        ls1_3= [p for p in ls if (p>=v_pos-20000) & (p<v_pos-10000)][-3:]
+        ls1_4= [p for p in ls if (p>=v_pos-40000) & (p<v_pos-20000)][-4:]
+        ls1_5= [p for p in ls if (p>=v_pos-50000) & (p<v_pos-40000)][-3:]
+        ls1_6= [p for p in ls if                  (p<v_pos-50000)][-3:]
+
+        ls2_0= [p for p in ls if (p>v_pos) & (p<=v_pos+2000)][-2:]
+        ls2_1= [p for p in ls if (p>v_pos+2000) & (p<=v_pos+5000)][:2]
+        ls2_2= [p for p in ls if (p>v_pos+5000) & (p<=v_pos+10000)][:3]
+        ls2_3= [p for p in ls if (p>v_pos+10000) & (p<=v_pos+20000)][:3]
+        ls2_4= [p for p in ls if (p>v_pos+20000) & (p<=v_pos+40000)][:4]
+        ls2_5= [p for p in ls if (p>v_pos+40000) & (p<=v_pos+50000)][:3]
+        ls2_6= [p for p in ls if (p>v_pos+50000)][:3]
+
+        ls_total_1=sorted(ls1_0+ls1_1+ls1_2+ls1_3+ls1_4+ls1_5+ls1_6)
+        ls_total_2=sorted(ls2_0+ls2_1+ls2_2+ls2_3+ls2_4+ls2_5+ls2_6)
+    
+    elif seq=='ul_ont_extreme':
+        ls=cnd_pos[abs(cnd_pos-v_pos)<300000] 
+
+        ls1_0= [p for p in ls if (p>=v_pos-10000) &  (p<v_pos)][:2]
+        ls1_1= [p for p in ls if (p>=v_pos-20000) & (p<v_pos-10000)][-2:]
+        ls1_2= [p for p in ls if (p>=v_pos-50000) & (p<v_pos-20000)][-3:]
+        ls1_3= [p for p in ls if (p>=v_pos-75000) & (p<v_pos-50000)][-3:]
+        ls1_4= [p for p in ls if (p>=v_pos-100000) & (p<v_pos-75000)][-4:]
+        ls1_5= [p for p in ls if (p>=v_pos-200000) & (p<v_pos-100000)][-4:]
+        ls1_6= [p for p in ls if                  (p<v_pos-200000)][-2:]
+
+        ls2_0= [p for p in ls if (p>v_pos) & (p<=v_pos+10000)][-2:]
+        ls2_1= [p for p in ls if (p>v_pos+10000) & (p<=v_pos+20000)][:2]
+        ls2_2= [p for p in ls if (p>v_pos+20000) & (p<=v_pos+50000)][:3]
+        ls2_3= [p for p in ls if (p>v_pos+50000) & (p<=v_pos+75000)][:3]
+        ls2_4= [p for p in ls if (p>v_pos+75000) & (p<=v_pos+100000)][:4]
+        ls2_5= [p for p in ls if (p>v_pos+100000) & (p<=v_pos+200000)][:4]
+        ls2_6= [p for p in ls if (p>v_pos+200000)][:2]
+
+        ls_total_1=sorted(ls1_0+ls1_1+ls1_2+ls1_3+ls1_4+ls1_5+ls1_6)
+        ls_total_2=sorted(ls2_0+ls2_1+ls2_2+ls2_3+ls2_4+ls2_5+ls2_6)
+        
+    elif seq=='new_pcb':
+        ls=cnd_pos[abs(cnd_pos-v_pos)<20000]
+   
+        ls1_0= [p for p in ls if (p>=v_pos-2000) &  (p<v_pos)][:4]
+        ls1_1= [p for p in ls if (p>=v_pos-5000) &  (p<v_pos-2000)][-5:]
+        ls1_2= [p for p in ls if (p>=v_pos-10000) & (p<v_pos-5000)][-5:]
+        ls1_3= [p for p in ls if (p>=v_pos-20000) & (p<v_pos-10000)][-6:]
+        
+        ls2_0= [p for p in ls if (p>v_pos) & (p<=v_pos+2000)][-4:]
+        ls2_1= [p for p in ls if (p>v_pos+2000) & (p<=v_pos+5000)][:5]
+        ls2_2= [p for p in ls if (p>v_pos+5000) & (p<=v_pos+10000)][:5]
+        ls2_3= [p for p in ls if (p>v_pos+10000) & (p<=v_pos+20000)][:6]
+        
+        ls_total_1=sorted(ls1_0+ls1_1+ls1_2+ls1_3)
+        ls_total_2=sorted(ls2_0+ls2_1+ls2_2+ls2_3)
+        
+    elif seq=='old_pcb':
+        ls=cnd_pos[abs(cnd_pos-v_pos)<20000] 
+
+        ls_total_1= [p for p in ls if (p>=v_pos-20000) &  (p<v_pos)][-20:]
+        ls_total_2= [p for p in ls if (p>v_pos) & (p<=v_pos+20000)][:20]
+    
+    return ls_total_1, ls_total_2
+    
+def get_nbr(dct, nbr_type='freq'):
+    chrom=dct['chrom']
+    start=max(dct['start']-50000,1)
+    end=dct['end']+50000
+
+    sam_path=dct['sam_path']
+    fasta_path=dct['fasta_path']
+    samfile = pysam.Samfile(sam_path, "rb")
+    fastafile=pysam.FastaFile(fasta_path)
+
+    tbx = pysam.TabixFile(dct['exclude_bed'])
+    exclude_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+    
+    
+    output_seq={}
+    
+    flag=0x4|0x100|0x200|0x400|0x800
+        
+    if nbr_type=='freq':
+        rlist=[s for s in fastafile.fetch(chrom,start-1,end-1)]
+        
+        for pcol in samfile.pileup(chrom,start-1,end-1,min_base_quality=0, flag_filter=flag,truncate=True):
+                
+                r=rlist[pcol.pos+1-start]                
+                if r in 'AGTC' and not in_bed(exclude_intervals,pcol.pos+1):
+                    n=pcol.get_num_aligned()
+                    seq=''.join([x[0] for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False,add_indels=True)]).upper()
+                    alt_freq=max([x[1] for x in Counter(seq).items() if (x[0]!=r and x[0] in 'AGTC')]+[0])/n
+
+                    if dct['threshold'][0]<=alt_freq and alt_freq<dct['threshold'][1] and n>=dct['mincov']:
+                        name=pcol.get_query_names()
+                        output_seq[pcol.pos+1]={n:base_to_num_map[s] for (n,s) in zip(name,seq)}
+                        output_seq[pcol.pos+1]['ref']=base_to_num_map[r]
+    
+    elif nbr_type=='gtruth':
+        gt_map={(0,0):0, (1,1):0, (2,2):0, (1,2):1, (2,1):1, (0,1):1, (1,0):1, (0,2):1, (2,0):1, (1,None):0,(None,1):0}
+        
+        bcf_in = VariantFile(dct['vcf_path'])
+        
+        ground_truth={}
+        for rec in bcf_in.fetch(chrom,start,end+1):
+            gt=rec.samples.items()[0][1].get('GT')
+            if gt_map[gt]:
+                if base_to_num_map[rec.alleles[gt[0]]]<4 and base_to_num_map[rec.alleles[gt[1]]]<4:
+                    ground_truth[rec.pos]=rec.ref
+                         
+        for pcol in samfile.pileup(chrom,start-1,end-1,min_base_quality=0, flag_filter=flag,truncate=True):
+            
+            if pcol.pos+1 in ground_truth:
+                n=pcol.get_num_aligned()
+                r=ground_truth[pcol.pos+1]
+
+
+                if n>=dct['mincov']:
+                    name=pcol.get_query_names()
+                    seq=''.join([x[0] for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False,add_indels=True)]).upper()
+                    output_seq[pcol.pos+1]={n:base_to_num_map[s] for (n,s) in zip(name,seq)}
+                    output_seq[pcol.pos+1]['ref']=base_to_num_map[r]
+    
+    return output_seq
+
+def get_snp_training_pileups(dct):
+    
+    chrom=dct['chrom']
+    start=dct['start']
+    end=dct['end']
+    
+    include_intervals = None
+    
+    tbx = pysam.TabixFile(dct['include_bed'])
+    include_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+
+    
+    sam_path=dct['sam_path']
+    fasta_path=dct['fasta_path']
+    vcf_path=dct['vcf_path']
+
+    bcf_in = VariantFile(vcf_path)
+    
+    gt_map={(0,0):0, (1,1):0, (2,2):0, (1,2):1, (2,1):1, (0,1):1, (1,0):1, (0,2):1,(2,0):1}
+    tr_pos={}
+    for rec in bcf_in.fetch(chrom,start,end+1):
+        gt=rec.samples.items()[0][1].get('GT')
+        if gt in gt_map:
+            if base_to_num_map[rec.alleles[gt[0]]]<4 and base_to_num_map[rec.alleles[gt[1]]]<4:
+                tr_pos[rec.pos]=(gt_map[gt],base_to_num_map[rec.alleles[gt[0]]],base_to_num_map[rec.alleles[gt[1]]])
+        
+        
+    nbr_size=20
+    
+    cnd_seq_freq=get_nbr(dct, nbr_type='freq')
+    cnd_pos_freq=np.array(list(cnd_seq_freq.keys()))
+    
+    cnd_seq_gtruth=get_nbr(dct, nbr_type='gtruth')
+    cnd_pos_gtruth=np.array(list(cnd_seq_gtruth.keys()))
+    
+    samfile = pysam.Samfile(sam_path, "rb")
+    fastafile=pysam.FastaFile(fasta_path)
+
+    ref_dict={j:s.upper() if s in 'AGTC' else '*' for j,s in zip(range(max(1,start-40),end+40+1),fastafile.fetch(chrom,max(1,start-40)-1,end+40)) }
+    
+    pileup_dict={}
+    
+    output={'pos':[],0:[],5:[],10:[],15:[],20:[],25:[]}
+    
+    flag=0x4|0x100|0x200|0x400|0x800
+        
+    for pcol in samfile.pileup(chrom,max(0,start-1),end,min_base_quality=0,\
+                                           flag_filter=flag,truncate=True):
+            
+            r=ref_dict[pcol.pos+1]
+            if in_bed(include_intervals, pcol.pos+1) and r in 'AGTC':
+                n=pcol.get_num_aligned()
+                
+                if n<=dct['maxcov'] and n>=dct['mincov'] and pcol.pos+1>=start and pcol.pos+1<=end:
+                    
+                    
+                    if pcol.pos+1 in tr_pos:
+                        seq=''.join([x[0] for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False,add_indels=True)]).upper()
+                        name=pcol.get_query_names()
+                        pileup_dict[pcol.pos+1]={n:base_to_num_map[s] for (n,s) in zip(name,seq)}
+                        output['pos'].append(pcol.pos+1)
+                    
+                    else:
+                        seq=''.join([x[0] for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False,add_indels=True)]).upper()
+                        alt_freq=max([x[1] for x in Counter(seq).items() if (x[0]!=r and x[0] in 'AGTC')]+[0])/n
+                        
+                        if alt_freq>=0.10:
+                            name=pcol.get_query_names()
+                            pileup_dict[pcol.pos+1]={n:base_to_num_map[s] for (n,s) in zip(name,seq)}
+
+                            if 0.10<=alt_freq<0.15:
+                                output[10].append(pcol.pos+1)
+                            elif 0.15<=alt_freq<0.20:
+                                output[15].append(pcol.pos+1)
+                            elif 0.20<=alt_freq<0.25:
+                                output[20].append(pcol.pos+1)
+                            elif 0.25<=alt_freq:
+                                output[25].append(pcol.pos+1)
+                    
+                        elif np.random.randint(2):
+                            seq=''.join([x[0] for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False,add_indels=True)]).upper()
+                            name=pcol.get_query_names()
+
+                            pileup_dict[pcol.pos+1]={n:base_to_num_map[s] for (n,s) in zip(name,seq)}
+
+                            if alt_freq<0.05:
+                                output[0].append(pcol.pos+1)
+                            
+                            elif 0.05<=alt_freq<0.10:
+                                output[5].append(pcol.pos+1)
+                            
+    
+    pileup_list={'pos':[],'neg':[]}
+
+    if output['pos']:
+        tr_len=len(output['pos'])    
+    else:
+        tr_len=1e16
+        
+    sizes={0:tr_len, 5:tr_len//3,10:tr_len//3,15:tr_len//3, 20:tr_len, 25:tr_len}
+
+    
+    for instance in ['pos',0,5,10,15,20,25]:
+        pos_list=output[instance]
+        
+        if pos_list:
+            if instance!='pos':
+                if sizes[instance]<len(output[instance]):
+                    perm=np.random.permutation(sizes[instance])
+                    pos_list=np.take(pos_list,perm,axis=0)
+        
+            for v_pos in pos_list:
+                for nbr_types in [(cnd_pos_freq, cnd_seq_freq),(cnd_pos_gtruth, cnd_seq_gtruth)]:
+                    cnd_pos, cnd_seq = nbr_types
+                    ls_total_1, ls_total_2 = get_cnd_pos(v_pos, cnd_pos, dct['seq'])
+
+                    nbr_dict={}
+
+                    rlist=[(v_pos,base_to_num_map[ref_dict[v_pos]]) ]
+
+                    sample=pileup_dict[v_pos].keys()
+
+                    nbr_dict={}
+
+                    tmp_mat=np.array([[4]*(len(ls_total_1)+1+len(ls_total_2))]*len(sample))
+
+                    for i,name in enumerate(sample):
+                        for j,nb_pos in enumerate(ls_total_1):
+                            try:    
+                                tmp_mat[i][j]=cnd_seq[nb_pos][name]
+                            except KeyError:
+                                pass
+
+                        tmp_mat[i][len(ls_total_1)]=pileup_dict[v_pos][name]
+
+                        for j,nb_pos in enumerate(ls_total_2):
+
+                            try:
+                                tmp_mat[i][j +len(ls_total_1)+1]=cnd_seq[nb_pos][name]
+                            except KeyError:
+                                pass
+                    for nb_pos in ls_total_1+ls_total_2:
+                                rlist.append((nb_pos, cnd_seq[nb_pos]['ref']))
+
+                    rlist=sorted(rlist, key=lambda x:x[0])
+                    total_rlist=np.array([x[1] for x in rlist])
+
+                    mat=np.dstack([np.sum(np.eye(5)[tmp_mat[tmp_mat[:,len(ls_total_1)]==i]],axis=0) for i in range(4)]).transpose(2,0,1)[:,:,:4]
+
+                    total_ref=np.eye(5)[total_rlist.astype(int)]
+                    total_ref[:,4]=0
+                    total_ref=total_ref[np.newaxis,:]
+                    mat=np.dstack([mat,np.zeros([4,mat.shape[1]])+np.eye(4)[base_to_num_map[ref_dict[v_pos]]][:,np.newaxis]])
+                    data=np.vstack([total_ref,np.multiply(mat,1-2*total_ref)])
+                    data=np.hstack([np.zeros([5,nbr_size-len(ls_total_1),5]),data,np.zeros([5, nbr_size-len(ls_total_2), 5])]).astype(np.int8)
+
+                    if instance=='pos':
+                        pileup_list['pos'].append((v_pos,tr_pos[v_pos][0],tr_pos[v_pos][1],tr_pos[v_pos][2], base_to_num_map[ref_dict[v_pos]],data))
+
+                    else:
+                        pileup_list['neg'].append((v_pos,0,base_to_num_map[ref_dict[v_pos]], base_to_num_map[ref_dict[v_pos]], base_to_num_map[ref_dict[v_pos]],data))
+               
+
+    return pileup_list, dct['type']
+
+def generate(params,mode='training'):
+    cores=params['cpu']
+
+    chrom_list=params['chrom_list']
+    chrom_lengths=params['chrom_lengths']
+    
+    threshold=params['threshold']
+    
+    
+    pool = mp.Pool(processes=args.cpu)
+    
+    in_dict_list=[]
+    
+    r_count=0
+    for chrom in chrom_list:
+        start, end=1, chrom_lengths[chrom]
+        for mbase in range(start, end,1000000):
+            r_count+=1
+            d = copy.deepcopy(params)
+            d['chrom']=chrom
+            d['start']=mbase
+            d['end']=min(end,mbase+1000000)
+            d['type']='test' if r_count%40==0 else 'train'
+            in_dict_list.append(d)
+        
+    results_dict=pool.imap_unordered(get_snp_training_pileups, in_dict_list)
+    
+    count=0
+    file_count=0
+    
+    pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count) 
+    pileup_file=open(pileup_file_name , "w")
+    
+    test_pileup_file_name=os.path.join(params['out_path'],'test.pileups') 
+    test_pileup_file=open(test_pileup_file_name , "w")
+    
+    total=len(in_dict_list)
+        
+    for result_pair in results_dict:
+        result, data_type=result_pair
+        
+        count+=1
+        if count%100==0:
+            pileup_file.close()
+            file_count+=1
+
+            pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count)
+            pileup_file=open(pileup_file_name , "w")
+            
+        for i in ['pos','neg']:
+            pileups=result[i]
+            if pileups:
+                if data_type=='train':
+                    for data in pileups:
+                            pos,gt,allele1,allele2,ref,mat=data
+                            mat=mat.reshape(-1)
+                            s='%s%d%d%d%d%d%s' %((11-len(str(pos)))*'0',pos,gt,allele1,allele2,ref,''.join([(6-len(x))*' '+x for x in mat.astype('<U6')]))
+                            pileup_file.write(s)
+                else:
+                    for data in pileups:
+                            pos,gt,allele1,allele2,ref,mat=data
+                            mat=mat.reshape(-1)
+                            s='%s%d%d%d%d%d%s' %((11-len(str(pos)))*'0',pos,gt,allele1,allele2,ref,''.join([(6-len(x))*' '+x for x in mat.astype('<U6')]))
+                            test_pileup_file.write(s)
+        print('%s: %d/%d regions done\n'%(str(datetime.datetime.now()),count,total),flush=True)
+
+    results_dict=None
+    elapsed=time.time()-t
+    
+    pileup_file.close()
+    test_pileup_file.close()
+    
+if __name__ == '__main__':
+    
+    t=time.time()
+
+    print('%s: Starting pileup generation.' %str(datetime.datetime.now()))
+    parser = argparse.ArgumentParser()
+
+    #-r chromosome region   -m mode   -bam bam file   -ref reference file   -vcf ground truth variants   -o output path
+    
+    parser.add_argument("-chrom", "--chrom", help="Chromosome region")
+    parser.add_argument("-seq",  "--sequencing",  help="Sequencing type, options are 'ont' and 'pacbio'", type=str, default='ont')
+    
+    parser.add_argument("-bam", "--bam", help="Bam file")
+    parser.add_argument("-ref", "--ref", help="Size")
+    parser.add_argument("-vcf", "--vcf", help="Ground truth variants")
+    
+    parser.add_argument("-o", "--output", help="Output path")
+    
+    parser.add_argument("-cpu", "--cpu", help="CPUs",type=int)
+    parser.add_argument("-include_bed", "--include_bed", help="Include BED file")
+    parser.add_argument("-exclude_bed", "--exclude_bed", help="Exclude BED file")
+
+    parser.add_argument("-start", "--start", help="start",type=int)
+    parser.add_argument("-end", "--end", help="end",type=int)
+    
+    parser.add_argument("-nbr_t",  "--neighbor_threshold",  help="SNP neighboring site thresholds with lower and upper bounds seperated by comma, for Nanopore reads '0.4,0.6' is recommended, for PacBio CCS anc CLR reads '0.3,0.7' and '0.3,0.6' are recommended respectively", type=str, default='0.4,0.6')
+    
+    parser.add_argument("-mincov",  "--mincov",  help="min coverage", type=int, default=8)
+    parser.add_argument("-maxcov",  "--maxcov",  help="max coverage", type=int, default=160)
+    
+    parser.add_argument('-wgs_contigs_type','--wgs_contigs_type', \
+                        help="""Options are "with_chr", "without_chr" and "all",\
+                        "with_chr" option will assume \
+                        human genome and run NanoCaller on chr1-22, "without_chr" will \
+                        run on chromosomes 1-22 if the BAM and reference genome files \
+                        use chromosome names without "chr". "all" option will run \
+                        NanoCaller on each contig present in reference genome FASTA file.""", \
+                        type=str, default='with_chr')
+    
+    args = parser.parse_args()
+    
+    os.makedirs(args.output, exist_ok=True)
+    
+    if args.chrom:
+        chrom_list= args.chrom.split()
+    
+    else:
+        if args.wgs_contigs_type=='with_chr':
+            chrom_list=['chr%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'without_chr':
+            chrom_list=['%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'all':
+            chrom_list=[]
+
+            try:
+                with open(args.ref+'.fai','r') as file:
+                    for line in file:
+                        chrom_list.append(line.split('\t')[0])
+
+            except FileNotFoundError:
+                print('%s: index file .fai required for reference genome file.\n' %str(datetime.datetime.now()), flush=True)
+                sys.exit(2)
+    
+    chrom_lengths={}
+    with open(args.ref+'.fai','r') as file:
+        for line in file:
+            chrom_lengths[line.split('\t')[0]]=int(line.split('\t')[1])
+            
+            
+    threshold=[float(args.neighbor_threshold.split(',')[0]), float(args.neighbor_threshold.split(',')[1])]
+    
+    
+    in_dict={'chrom_list':chrom_list,'chrom_lengths':chrom_lengths,\
+         'sam_path':args.bam, 'fasta_path':args.ref, 'vcf_path':args.vcf,\
+             'out_path':args.output,'seq':args.sequencing, \
+             'threshold':threshold, 'cpu':args.cpu, 'include_bed':args.include_bed,'exclude_bed':args.exclude_bed,\
+            'mincov':args.mincov,'maxcov':args.maxcov}    
+    
+    generate(in_dict)
+    elapsed=time.time()-t
+    print ('%s: Total Time Elapsed: %.2f seconds' %(str(datetime.datetime.now()), elapsed))
\ No newline at end of file
diff --git a/misc/training/generate_indel_pileups.py b/misc/training/generate_indel_pileups.py
new file mode 100644
index 0000000..cd1786d
--- /dev/null
+++ b/misc/training/generate_indel_pileups.py
@@ -0,0 +1,400 @@
+import sys,pysam, time,os,copy,argparse,subprocess, random, re, datetime
+import numpy as np
+import multiprocessing as mp
+from pysam import VariantFile
+from subprocess import Popen, PIPE, STDOUT
+from Bio import pairwise2
+from intervaltree import Interval, IntervalTree
+
+gt_map={(0,0):0, (1,1):1, (0,1):2, (1,0):2,(1,2):3, (2,1):3}
+
+mapping={'A':0,'G':1,'T':2,'C':3,'-':4}
+rev_base_map={0:'A',1:'G',2:'T',3:'C',4:'-'}
+
+allele_map={('I','I'):0, ('D','D'):1, ('N','I'):2, ('I','N'):3, ('N','D'):4, ('D','N'):5, \
+           ('I','D'):6, ('D','I'):7, ('N','N'):8}
+
+def pairwise(x,y):
+    alignments = pairwise2.align.globalms(x, y, 2, -1.0, -0.9, -0.1)
+
+    return alignments
+
+def msa(seq_list, ref, v_pos, mincov, maxcov):
+    np.random.seed(812)
+    sample=list(seq_list.keys())
+            
+    if len(sample) > maxcov:
+        sample=random.sample(sample,min(len(sample),maxcov))
+
+    sample=sorted(sample)
+    
+    fa_tmp_file=''.join(['>%s_SEQ\n%s\n'%(read_name,seq_list[read_name]) for read_name in sample])
+
+
+    fa_tmp_file+='>ref_SEQ\n%s' %ref
+    
+    gap_penalty=1.0
+    msa_process =Popen(['muscle', '-quiet','-gapopen','%.1f' %gap_penalty,'-maxiters', '1' ,'-diags1'], stdout=PIPE, stdin=PIPE, stderr=PIPE)
+    hap_file=msa_process.communicate(input=fa_tmp_file.encode('utf-8'))
+
+    if len(hap_file)==0:
+        print('hapfile length 0')
+
+
+    tmp=hap_file[0].decode('utf-8')[1:].replace('\n','').split('>')
+
+    zz_0=[]
+    for seq in tmp:
+        p1,p2=seq.split('_SEQ')
+        if p1!='ref':
+            zz_0.append(p2[:128])
+        else:
+            ref_real_0=p2
+
+    if len(zz_0)<mincov:
+        return (0,0,None,None,None)
+        
+    try:
+        ref_real_0_mat=np.eye(5)[[mapping[x] for x in ref_real_0[:128]]].transpose()
+    except UnboundLocalError:
+        return (0,0,None)
+    
+    mat=np.array([[mapping[c] for c in x] for x in zz_0])
+    h0_mat=np.sum(np.eye(5)[mat],axis=0).transpose()
+
+    return (1,1,np.dstack([h0_mat, ref_real_0_mat]))
+
+def get_indel_training_pileups(dct):
+    window_before,window_after=0,160
+    
+    if dct['seq']=='pacbio':
+        window_after=260
+    
+    chrom=dct['chrom']
+    start=dct['start']
+    end=dct['end']
+    sam_path=dct['sam_path']
+    fasta_path=dct['fasta_path']
+    samfile = pysam.Samfile(sam_path, "rb")
+    fastafile=pysam.FastaFile(fasta_path)
+    mincov,maxcov=dct['mincov'],dct['maxcov']
+    
+    include_intervals, exclude_intervals=None, None
+    
+    
+    if dct['include_bed']:
+        tbx = pysam.TabixFile(dct['include_bed'])
+        include_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+        
+        def in_bed(tree,pos):            
+            return tree.overlaps(pos)
+        
+        include_intervals=IntervalTree(include_intervals.overlap(start,end))
+        
+        if not include_intervals:
+            return None
+    
+        else:
+            start=max(start, min(x[0] for x in include_intervals))
+            end=min(end, max(x[1] for x in include_intervals))
+        
+    else:
+        def in_bed(tree, pos):
+            return True
+    
+    if dct['exclude_bed']:
+        tbx = pysam.TabixFile(dct['exclude_bed'])
+        try:
+            exclude_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+        
+            def ex_bed(tree, pos):
+                return tree.overlaps(pos)
+        
+        except ValueError:
+            def ex_bed(tree, pos):
+                return False 
+    else:
+        def ex_bed(tree, pos):
+            return False 
+    
+
+        
+    ref_dict={j:s.upper() if s in 'AGTC' else '' for j,s in zip(range(max(1,start-200),end+400+1),fastafile.fetch(chrom,max(1,start-200)-1,end+400)) }
+    
+    hap_dict={1:[],2:[]}
+    
+    for pread in samfile.fetch(chrom, max(0,dct['start']-100), dct['end']+100):
+        if pread.has_tag('HP'):
+            hap_dict[pread.get_tag('HP')].append(pread.qname)
+
+    hap_reads_0=set(hap_dict[1])
+    hap_reads_1=set(hap_dict[2])
+    
+    bcf_in = VariantFile(dct['vcf_path'])
+    tr_var={}
+    for rec in bcf_in.fetch(chrom,start,end+1):
+            gt=rec.samples.items()[0][1].get('GT')
+            tr_var[rec.pos-40]=gt_map[gt]
+            tr_var[rec.pos-30]=gt_map[gt]
+            tr_var[rec.pos-20]=gt_map[gt]
+            tr_var[rec.pos-10]=gt_map[gt]
+            tr_var[rec.pos]=gt_map[gt]
+            
+    tr_pos=set(tr_var.keys())        
+
+    output={'pos':[],'high':[],'low':[]}
+
+    for pcol in samfile.pileup(chrom,max(0,start-100),end,min_base_quality=0,\
+                                           flag_filter=0x4|0x100|0x200|0x400|0x800,truncate=True):
+        v_pos=pcol.pos+1
+        if v_pos in tr_pos:
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+            len_seq_0=len(read_names_0)
+            len_seq_1=len(read_names_1)
+            
+            if len(read_names_0)>=dct['mincov'] and len(read_names_1)>=dct['mincov']:
+                output['pos'].append(pcol.pos+1)
+                        
+        if in_bed(include_intervals, v_pos) and not ex_bed(exclude_intervals, v_pos):                
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+            len_seq_0=len(read_names_0)
+            len_seq_1=len(read_names_1)
+            
+            if len_seq_0>=dct['mincov']  and len_seq_1>=dct['mincov']:
+                seq=[x for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False, add_indels=True)]
+
+                tmp_seq_0=''.join([s for n,s in zip(read_names,seq) if n in read_names_0])
+                tmp_seq_1=''.join([s for n,s in zip(read_names,seq) if n in read_names_1])
+
+                del_freq_0=(tmp_seq_0.count('-'))/len_seq_0 if len_seq_0>0 else 0
+                ins_freq_0=tmp_seq_0.count('+')/len_seq_0 if len_seq_0>0 else 0
+
+                del_freq_1=(tmp_seq_1.count('-'))/len_seq_1 if len_seq_1>0 else 0
+                ins_freq_1=tmp_seq_1.count('+')/len_seq_1 if len_seq_1>0 else 0
+
+                if 0.3<=del_freq_0 or 0.3<=del_freq_1 or 0.3<=ins_freq_0 or 0.3<=ins_freq_1:
+                    output['high'].append(pcol.pos+1)
+
+                elif del_freq_0<=0.2 and del_freq_1<=0.2 and ins_freq_0<=0.1 and ins_freq_1<=0.1 and np.random.randint(100)==0:
+                    output['low'].append(pcol.pos+1)
+    
+    
+    if output['pos']:
+        tr_len=len(output['pos'])    
+    else:
+        tr_len=20
+        
+    sizes={'high':tr_len, 'low':tr_len}
+
+    output['pos']=set(output['pos'])
+    
+    
+    for i in ['high','low']:
+        if sizes[i]<len(output[i]):
+            perm=np.random.permutation(sizes[i])
+            var_list=np.take(output[i],perm,axis=0)
+            output[i]=set(var_list)
+            
+            
+    
+    total_variants=set.union(output['pos'],output['high'],output['low'])
+    
+    data=[]
+
+    for pcol in samfile.pileup(chrom,max(0,start-100),end,min_base_quality=0, flag_filter=0x4|0x100|0x200|0x400|0x800,truncate=True):
+        
+        v_pos=pcol.pos+1
+        
+        if v_pos in total_variants:
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+           
+
+            d={'hap0':{},'hap1':{}}
+            d_tot={}
+                        
+            ref=''.join([ref_dict[p] for p in range(v_pos-window_before,v_pos+window_after+1)])
+            
+            if len(ref)<window_after:
+                continue
+                    
+            for pread in pcol.pileups:
+                dt=pread.alignment.query_sequence[max(0,pread.query_position_or_next-window_before):pread.query_position_or_next+window_after]                    
+                d_tot[pread.alignment.qname]=dt
+                    
+                if pread.alignment.qname in read_names_0:
+                    d['hap0'][pread.alignment.qname]=dt
+
+                elif pread.alignment.qname in read_names_1:
+                    d['hap1'][pread.alignment.qname]=dt    
+            
+            seq_list=d['hap0']
+            flag0,indel_flag0,data_0=msa(seq_list,ref,v_pos,2,dct['maxcov'])
+
+            seq_list=d['hap1']
+            flag1,indel_flag1,data_1=msa(seq_list,ref,v_pos,2,dct['maxcov'])
+            
+            seq_list = d_tot
+            flag_total,indel_flag_total,data_total=msa(seq_list,ref,v_pos,dct['mincov'],dct['maxcov'])
+
+            if flag0 and flag1 and flag_total:
+                gt= tr_var[v_pos] if v_pos in tr_pos else 0
+                data.append((v_pos,gt,(data_0,data_1,data_total)))
+    
+    return data
+
+def generate(params,mode='training'):
+    cores=params['cpu']
+
+    chrom_list=params['chrom_list']
+    chrom_lengths=params['chrom_lengths']    
+    
+    pool = mp.Pool(processes=args.cpu)
+    
+    in_dict_list=[]
+    
+    r_count=0
+
+    for chrom in chrom_list:
+        start, end=1, chrom_lengths[chrom]
+        for mbase in range(start, end,100000):
+            d = copy.deepcopy(params)
+            d['chrom']=chrom
+            d['start']=mbase
+            d['end']=min(end,mbase+100000)
+            in_dict_list.append(d)
+            
+    results_dict=pool.imap_unordered(get_indel_training_pileups, in_dict_list)
+    
+    count=0
+    file_count=0
+    
+    pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count) 
+    pileup_file=open(pileup_file_name , "w")
+    
+    test_pileup_file_name=os.path.join(params['out_path'],'test.pileups') 
+    test_pileup_file=open(test_pileup_file_name , "w")
+    
+    total=len(in_dict_list)
+    for result in results_dict:        
+        count+=1
+        if count%1000==0:
+            pileup_file.close()
+            file_count+=1
+
+            pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count)
+            pileup_file=open(pileup_file_name , "w")
+            
+        if result:
+            if count%40!=0:
+                for data in result:
+                    pos,gt,phased_data=data
+                    mat_0=phased_data[0].reshape(-1).astype(np.int16)
+                    mat_1=phased_data[1].reshape(-1).astype(np.int16)
+                    mat_2=phased_data[2].reshape(-1).astype(np.int16)
+                    s='%s%d%d%s%s%s' %((11-len(str(pos)))*'0',pos,gt, ''.join([(4-len(x))*' '+x for x in mat_0.astype('<U4')]), ''.join([(4-len(x))*' '+x for x in mat_1.astype('<U4')]),''.join([(4-len(x))*' '+x for x in mat_2.astype('<U4')]))
+                    pileup_file.write(s)
+            else:
+                for data in result:
+                    pos,gt,phased_data=data
+                    mat_0=phased_data[0].reshape(-1).astype(np.int16)
+                    mat_1=phased_data[1].reshape(-1).astype(np.int16)
+                    mat_2=phased_data[2].reshape(-1).astype(np.int16)
+                    s='%s%d%d%s%s%s' %((11-len(str(pos)))*'0',pos,gt, ''.join([(4-len(x))*' '+x for x in mat_0.astype('<U4')]), ''.join([(4-len(x))*' '+x for x in mat_1.astype('<U4')]),''.join([(4-len(x))*' '+x for x in mat_2.astype('<U4')]))
+                    test_pileup_file.write(s)
+        print('%s: %d/%d regions done\n'%(str(datetime.datetime.now()),count,total),flush=True)
+
+    results_dict=None
+    elapsed=time.time()-t
+    
+    pileup_file.close()
+    test_pileup_file.close()
+    
+if __name__ == '__main__':
+    
+    t=time.time()
+
+    print('%s: Starting pileup generation.' %str(datetime.datetime.now()))
+    parser = argparse.ArgumentParser()
+
+    #-r chromosome region   -m mode   -bam bam file   -ref reference file   -vcf ground truth variants   -o output path
+    
+    parser.add_argument("-chrom", "--chrom", help="Chromosome region")
+    parser.add_argument("-seq",  "--sequencing",  help="Sequencing type, options are 'ont' and 'pacbio'", type=str, default='ont')
+    
+    parser.add_argument("-bam", "--bam", help="Bam file")
+    parser.add_argument("-ref", "--ref", help="Size")
+    parser.add_argument("-vcf", "--vcf", help="Ground truth variants")
+    
+    parser.add_argument("-o", "--output", help="Output path")
+    
+    parser.add_argument("-cpu", "--cpu", help="CPUs",type=int)
+    parser.add_argument("-include_bed", "--include_bed", help="Include BED file")
+    parser.add_argument("-exclude_bed", "--exclude_bed", help="Exclude BED file")
+
+    parser.add_argument("-start", "--start", help="start",type=int)
+    parser.add_argument("-end", "--end", help="end",type=int)
+    
+    parser.add_argument("-nbr_t",  "--neighbor_threshold",  help="SNP neighboring site thresholds with lower and upper bounds seperated by comma, for Nanopore reads '0.4,0.6' is recommended, for PacBio CCS anc CLR reads '0.3,0.7' and '0.3,0.6' are recommended respectively", type=str, default='0.4,0.6')
+    
+    parser.add_argument("-mincov",  "--mincov",  help="min coverage", type=int, default=8)
+    parser.add_argument("-maxcov",  "--maxcov",  help="max coverage", type=int, default=160)
+    
+    parser.add_argument('-wgs_contigs_type','--wgs_contigs_type', \
+                        help="""Options are "with_chr", "without_chr" and "all",\
+                        "with_chr" option will assume \
+                        human genome and run NanoCaller on chr1-22, "without_chr" will \
+                        run on chromosomes 1-22 if the BAM and reference genome files \
+                        use chromosome names without "chr". "all" option will run \
+                        NanoCaller on each contig present in reference genome FASTA file.""", \
+                        type=str, default='with_chr')
+    
+    args = parser.parse_args()
+    
+    os.makedirs(args.output, exist_ok=True)
+    
+    if args.chrom:
+        chrom_list= args.chrom.split()
+    
+    else:
+        if args.wgs_contigs_type=='with_chr':
+            chrom_list=['chr%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'without_chr':
+            chrom_list=['%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'all':
+            chrom_list=[]
+
+            try:
+                with open(args.ref+'.fai','r') as file:
+                    for line in file:
+                        chrom_list.append(line.split('\t')[0])
+
+            except FileNotFoundError:
+                print('%s: index file .fai required for reference genome file.\n' %str(datetime.datetime.now()), flush=True)
+                sys.exit(2)
+    
+    chrom_lengths={}
+    with open(args.ref+'.fai','r') as file:
+        for line in file:
+            chrom_lengths[line.split('\t')[0]]=int(line.split('\t')[1])
+            
+            
+    threshold=[float(args.neighbor_threshold.split(',')[0]), float(args.neighbor_threshold.split(',')[1])]
+    
+    
+    in_dict={'chrom_list':chrom_list,'chrom_lengths':chrom_lengths,\
+         'sam_path':args.bam, 'fasta_path':args.ref, 'vcf_path':args.vcf,\
+             'out_path':args.output,'seq':args.sequencing, \
+             'threshold':threshold, 'cpu':args.cpu, 'include_bed':args.include_bed,'exclude_bed':args.exclude_bed,\
+            'mincov':args.mincov,'maxcov':args.maxcov}    
+
+    generate(in_dict)
+    elapsed=time.time()-t
+    print ('%s: Total Time Elapsed: %.2f seconds' %(str(datetime.datetime.now()), elapsed))
\ No newline at end of file
diff --git a/misc/training/generate_indel_pileups_hifi.py b/misc/training/generate_indel_pileups_hifi.py
new file mode 100644
index 0000000..79b656a
--- /dev/null
+++ b/misc/training/generate_indel_pileups_hifi.py
@@ -0,0 +1,397 @@
+import sys,pysam, time,os,copy,argparse,subprocess, random, re, datetime
+import numpy as np
+import multiprocessing as mp
+from pysam import VariantFile
+from subprocess import Popen, PIPE, STDOUT
+from Bio import pairwise2
+from intervaltree import Interval, IntervalTree
+
+gt_map={(0,0):0, (1,1):1, (0,1):2, (1,0):2,(1,2):3, (2,1):3}
+
+mapping={'A':0,'G':1,'T':2,'C':3,'-':4}
+rev_base_map={0:'A',1:'G',2:'T',3:'C',4:'-'}
+
+allele_map={('I','I'):0, ('D','D'):1, ('N','I'):2, ('I','N'):3, ('N','D'):4, ('D','N'):5, \
+           ('I','D'):6, ('D','I'):7, ('N','N'):8}
+
+def pairwise(x,y):
+    alignments = pairwise2.align.globalms(x, y, 2, -1.0, -0.9, -0.1)
+
+    return alignments
+
+def msa(seq_list, ref, v_pos, mincov, maxcov):
+    np.random.seed(812)
+    sample=list(seq_list.keys())
+            
+    if len(sample) > maxcov:
+        sample=random.sample(sample,min(len(sample),maxcov))
+
+    sample=sorted(sample)
+    
+    fa_tmp_file=''.join(['>%s_SEQ\n%s\n'%(read_name,seq_list[read_name]) for read_name in sample])
+
+
+    fa_tmp_file+='>ref_SEQ\n%s' %ref
+    
+    gap_penalty=1.0
+    msa_process =Popen(['muscle', '-quiet','-gapopen','%.1f' %gap_penalty,'-maxiters', '1' ,'-diags1'], stdout=PIPE, stdin=PIPE, stderr=PIPE)
+    hap_file=msa_process.communicate(input=fa_tmp_file.encode('utf-8'))
+
+    if len(hap_file)==0:
+        print('hapfile length 0')
+
+
+    tmp=hap_file[0].decode('utf-8')[1:].replace('\n','').split('>')
+
+    zz_0=[]
+    for seq in tmp:
+        p1,p2=seq.split('_SEQ')
+        if p1!='ref':
+            zz_0.append(p2[:128])
+        else:
+            ref_real_0=p2
+
+    if len(zz_0)<mincov:
+        return (0,0,None,None,None)
+        
+    try:
+        ref_real_0_mat=np.eye(5)[[mapping[x] for x in ref_real_0[:128]]].transpose()
+    except UnboundLocalError:
+        return (0,0,None)
+    
+    mat=np.array([[mapping[c] for c in x] for x in zz_0])
+    h0_mat=np.sum(np.eye(5)[mat],axis=0).transpose()
+
+    return (1,1,np.dstack([h0_mat, ref_real_0_mat]))
+
+def get_indel_training_pileups(dct):
+    window_before,window_after=0,160
+    
+    if dct['seq']=='pacbio':
+        window_after=260
+    
+    chrom=dct['chrom']
+    start=dct['start']
+    end=dct['end']
+    sam_path=dct['sam_path']
+    fasta_path=dct['fasta_path']
+    samfile = pysam.Samfile(sam_path, "rb")
+    fastafile=pysam.FastaFile(fasta_path)
+    mincov,maxcov=dct['mincov'],dct['maxcov']
+    
+    include_intervals, exclude_intervals=None, None
+    
+    
+    if dct['include_bed']:
+        tbx = pysam.TabixFile(dct['include_bed'])
+        include_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+        
+        def in_bed(tree,pos):            
+            return tree.overlaps(pos)
+        
+        include_intervals=IntervalTree(include_intervals.overlap(start,end))
+        
+        if not include_intervals:
+            return None
+    
+        else:
+            start=max(start, min(x[0] for x in include_intervals))
+            end=min(end, max(x[1] for x in include_intervals))
+        
+    else:
+        def in_bed(tree, pos):
+            return True
+    
+    if dct['exclude_bed']:
+        tbx = pysam.TabixFile(dct['exclude_bed'])
+        try:
+            exclude_intervals=IntervalTree(Interval(int(row[1]), int(row[2]), "%s" % (row[1])) for row in tbx.fetch(chrom, parser=pysam.asBed()))
+        
+            def ex_bed(tree, pos):
+                return tree.overlaps(pos)
+        
+        except ValueError:
+            def ex_bed(tree, pos):
+                return False 
+    else:
+        def ex_bed(tree, pos):
+            return False 
+    
+
+        
+    ref_dict={j:s.upper() if s in 'AGTC' else '' for j,s in zip(range(max(1,start-200),end+400+1),fastafile.fetch(chrom,max(1,start-200)-1,end+400)) }
+    
+    hap_dict={1:[],2:[]}
+    
+    for pread in samfile.fetch(chrom, max(0,dct['start']-100), dct['end']+100):
+        if pread.has_tag('HP'):
+            hap_dict[pread.get_tag('HP')].append(pread.qname)
+
+    hap_reads_0=set(hap_dict[1])
+    hap_reads_1=set(hap_dict[2])
+    
+    bcf_in = VariantFile(dct['vcf_path'])
+    tr_var={}
+    for rec in bcf_in.fetch(chrom,start,end+1):
+            gt=rec.samples.items()[0][1].get('GT')
+            tr_var[rec.pos-10]=gt_map[gt]
+            tr_var[rec.pos]=gt_map[gt]
+            
+    tr_pos=set(tr_var.keys())        
+
+    output={'pos':[],'high':[],'low':[]}
+
+    for pcol in samfile.pileup(chrom,max(0,start-100),end,min_base_quality=0,\
+                                           flag_filter=0x4|0x100|0x200|0x400|0x800,truncate=True):
+        v_pos=pcol.pos+1
+        if v_pos in tr_pos:
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+            len_seq_0=len(read_names_0)
+            len_seq_1=len(read_names_1)
+            
+            if len(read_names_0)>=dct['mincov'] and len(read_names_1)>=dct['mincov']:
+                output['pos'].append(pcol.pos+1)
+                        
+        if in_bed(include_intervals, v_pos) and not ex_bed(exclude_intervals, v_pos):                
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+            len_seq_0=len(read_names_0)
+            len_seq_1=len(read_names_1)
+            
+            if len_seq_0>=dct['mincov']  and len_seq_1>=dct['mincov']:
+                seq=[x for x in pcol.get_query_sequences( mark_matches=False, mark_ends=False, add_indels=True)]
+
+                tmp_seq_0=''.join([s for n,s in zip(read_names,seq) if n in read_names_0])
+                tmp_seq_1=''.join([s for n,s in zip(read_names,seq) if n in read_names_1])
+
+                del_freq_0=(tmp_seq_0.count('-'))/len_seq_0 if len_seq_0>0 else 0
+                ins_freq_0=tmp_seq_0.count('+')/len_seq_0 if len_seq_0>0 else 0
+
+                del_freq_1=(tmp_seq_1.count('-'))/len_seq_1 if len_seq_1>0 else 0
+                ins_freq_1=tmp_seq_1.count('+')/len_seq_1 if len_seq_1>0 else 0
+
+                if 0.3<=del_freq_0 or 0.3<=del_freq_1 or 0.3<=ins_freq_0 or 0.3<=ins_freq_1:
+                    output['high'].append(pcol.pos+1)
+
+                elif del_freq_0<=0.2 and del_freq_1<=0.2 and ins_freq_0<=0.1 and ins_freq_1<=0.1 and np.random.randint(100)==0:
+                    output['low'].append(pcol.pos+1)
+    
+    
+    if output['pos']:
+        tr_len=len(output['pos'])    
+    else:
+        tr_len=20
+        
+    sizes={'high':tr_len, 'low':tr_len}
+
+    output['pos']=set(output['pos'])
+    
+    
+    for i in ['high','low']:
+        if sizes[i]<len(output[i]):
+            perm=np.random.permutation(sizes[i])
+            var_list=np.take(output[i],perm,axis=0)
+            output[i]=set(var_list)
+            
+            
+    
+    total_variants=set.union(output['pos'],output['high'],output['low'])
+    
+    data=[]
+
+    for pcol in samfile.pileup(chrom,max(0,start-100),end,min_base_quality=0, flag_filter=0x4|0x100|0x200|0x400|0x800,truncate=True):
+        
+        v_pos=pcol.pos+1
+        
+        if v_pos in total_variants:
+            read_names=pcol.get_query_names()
+            read_names_0=set(read_names) & hap_reads_0
+            read_names_1=set(read_names) & hap_reads_1
+           
+
+            d={'hap0':{},'hap1':{}}
+            d_tot={}
+                        
+            ref=''.join([ref_dict[p] for p in range(v_pos-window_before,v_pos+window_after+1)])
+            
+            if len(ref)<window_after:
+                continue
+                    
+            for pread in pcol.pileups:
+                dt=pread.alignment.query_sequence[max(0,pread.query_position_or_next-window_before):pread.query_position_or_next+window_after]                    
+                d_tot[pread.alignment.qname]=dt
+                    
+                if pread.alignment.qname in read_names_0:
+                    d['hap0'][pread.alignment.qname]=dt
+
+                elif pread.alignment.qname in read_names_1:
+                    d['hap1'][pread.alignment.qname]=dt    
+            
+            seq_list=d['hap0']
+            flag0,indel_flag0,data_0=msa(seq_list,ref,v_pos,2,dct['maxcov'])
+
+            seq_list=d['hap1']
+            flag1,indel_flag1,data_1=msa(seq_list,ref,v_pos,2,dct['maxcov'])
+            
+            seq_list = d_tot
+            flag_total,indel_flag_total,data_total=msa(seq_list,ref,v_pos,dct['mincov'],dct['maxcov'])
+
+            if flag0 and flag1 and flag_total:
+                gt= tr_var[v_pos] if v_pos in tr_pos else 0
+                data.append((v_pos,gt,(data_0,data_1,data_total)))
+    
+    return data
+
+def generate(params,mode='training'):
+    cores=params['cpu']
+
+    chrom_list=params['chrom_list']
+    chrom_lengths=params['chrom_lengths']    
+    
+    pool = mp.Pool(processes=args.cpu)
+    
+    in_dict_list=[]
+    
+    r_count=0
+
+    for chrom in chrom_list:
+        start, end=1, chrom_lengths[chrom]
+        for mbase in range(start, end,100000):
+            d = copy.deepcopy(params)
+            d['chrom']=chrom
+            d['start']=mbase
+            d['end']=min(end,mbase+100000)
+            in_dict_list.append(d)
+            
+    results_dict=pool.imap_unordered(get_indel_training_pileups, in_dict_list)
+    
+    count=0
+    file_count=0
+    
+    pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count) 
+    pileup_file=open(pileup_file_name , "w")
+    
+    test_pileup_file_name=os.path.join(params['out_path'],'test.pileups') 
+    test_pileup_file=open(test_pileup_file_name , "w")
+    
+    total=len(in_dict_list)
+    for result in results_dict:        
+        count+=1
+        if count%1000==0:
+            pileup_file.close()
+            file_count+=1
+
+            pileup_file_name=os.path.join(params['out_path'],'train.pileups.%d' %file_count)
+            pileup_file=open(pileup_file_name , "w")
+            
+        if result:
+            if count%40!=0:
+                for data in result:
+                    pos,gt,phased_data=data
+                    mat_0=phased_data[0].reshape(-1).astype(np.int16)
+                    mat_1=phased_data[1].reshape(-1).astype(np.int16)
+                    mat_2=phased_data[2].reshape(-1).astype(np.int16)
+                    s='%s%d%d%s%s%s' %((11-len(str(pos)))*'0',pos,gt, ''.join([(4-len(x))*' '+x for x in mat_0.astype('<U4')]), ''.join([(4-len(x))*' '+x for x in mat_1.astype('<U4')]),''.join([(4-len(x))*' '+x for x in mat_2.astype('<U4')]))
+                    pileup_file.write(s)
+            else:
+                for data in result:
+                    pos,gt,phased_data=data
+                    mat_0=phased_data[0].reshape(-1).astype(np.int16)
+                    mat_1=phased_data[1].reshape(-1).astype(np.int16)
+                    mat_2=phased_data[2].reshape(-1).astype(np.int16)
+                    s='%s%d%d%s%s%s' %((11-len(str(pos)))*'0',pos,gt, ''.join([(4-len(x))*' '+x for x in mat_0.astype('<U4')]), ''.join([(4-len(x))*' '+x for x in mat_1.astype('<U4')]),''.join([(4-len(x))*' '+x for x in mat_2.astype('<U4')]))
+                    test_pileup_file.write(s)
+        print('%s: %d/%d regions done\n'%(str(datetime.datetime.now()),count,total),flush=True)
+
+    results_dict=None
+    elapsed=time.time()-t
+    
+    pileup_file.close()
+    test_pileup_file.close()
+    
+if __name__ == '__main__':
+    
+    t=time.time()
+
+    print('%s: Starting pileup generation.' %str(datetime.datetime.now()))
+    parser = argparse.ArgumentParser()
+
+    #-r chromosome region   -m mode   -bam bam file   -ref reference file   -vcf ground truth variants   -o output path
+    
+    parser.add_argument("-chrom", "--chrom", help="Chromosome region")
+    parser.add_argument("-seq",  "--sequencing",  help="Sequencing type, options are 'ont' and 'pacbio'", type=str, default='ont')
+    
+    parser.add_argument("-bam", "--bam", help="Bam file")
+    parser.add_argument("-ref", "--ref", help="Size")
+    parser.add_argument("-vcf", "--vcf", help="Ground truth variants")
+    
+    parser.add_argument("-o", "--output", help="Output path")
+    
+    parser.add_argument("-cpu", "--cpu", help="CPUs",type=int)
+    parser.add_argument("-include_bed", "--include_bed", help="Include BED file")
+    parser.add_argument("-exclude_bed", "--exclude_bed", help="Exclude BED file")
+
+    parser.add_argument("-start", "--start", help="start",type=int)
+    parser.add_argument("-end", "--end", help="end",type=int)
+    
+    parser.add_argument("-nbr_t",  "--neighbor_threshold",  help="SNP neighboring site thresholds with lower and upper bounds seperated by comma, for Nanopore reads '0.4,0.6' is recommended, for PacBio CCS anc CLR reads '0.3,0.7' and '0.3,0.6' are recommended respectively", type=str, default='0.4,0.6')
+    
+    parser.add_argument("-mincov",  "--mincov",  help="min coverage", type=int, default=8)
+    parser.add_argument("-maxcov",  "--maxcov",  help="max coverage", type=int, default=160)
+    
+    parser.add_argument('-wgs_contigs_type','--wgs_contigs_type', \
+                        help="""Options are "with_chr", "without_chr" and "all",\
+                        "with_chr" option will assume \
+                        human genome and run NanoCaller on chr1-22, "without_chr" will \
+                        run on chromosomes 1-22 if the BAM and reference genome files \
+                        use chromosome names without "chr". "all" option will run \
+                        NanoCaller on each contig present in reference genome FASTA file.""", \
+                        type=str, default='with_chr')
+    
+    args = parser.parse_args()
+    
+    os.makedirs(args.output, exist_ok=True)
+    
+    if args.chrom:
+        chrom_list= args.chrom.split()
+    
+    else:
+        if args.wgs_contigs_type=='with_chr':
+            chrom_list=['chr%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'without_chr':
+            chrom_list=['%d' %d for d in range(1,23)]
+
+        elif args.wgs_contigs_type == 'all':
+            chrom_list=[]
+
+            try:
+                with open(args.ref+'.fai','r') as file:
+                    for line in file:
+                        chrom_list.append(line.split('\t')[0])
+
+            except FileNotFoundError:
+                print('%s: index file .fai required for reference genome file.\n' %str(datetime.datetime.now()), flush=True)
+                sys.exit(2)
+    
+    chrom_lengths={}
+    with open(args.ref+'.fai','r') as file:
+        for line in file:
+            chrom_lengths[line.split('\t')[0]]=int(line.split('\t')[1])
+            
+            
+    threshold=[float(args.neighbor_threshold.split(',')[0]), float(args.neighbor_threshold.split(',')[1])]
+    
+    
+    in_dict={'chrom_list':chrom_list,'chrom_lengths':chrom_lengths,\
+         'sam_path':args.bam, 'fasta_path':args.ref, 'vcf_path':args.vcf,\
+             'out_path':args.output,'seq':args.sequencing, \
+             'threshold':threshold, 'cpu':args.cpu, 'include_bed':args.include_bed,'exclude_bed':args.exclude_bed,\
+            'mincov':args.mincov,'maxcov':args.maxcov}    
+
+    generate(in_dict)
+    elapsed=time.time()-t
+    print ('%s: Total Time Elapsed: %.2f seconds' %(str(datetime.datetime.now()), elapsed))
diff --git a/misc/training/model_architect.py b/misc/training/model_architect.py
new file mode 100644
index 0000000..e944ac7
--- /dev/null
+++ b/misc/training/model_architect.py
@@ -0,0 +1,182 @@
+import tensorflow as tf
+import numpy as np
+        
+def conv2d(x, W, b, strides=1,pad_type='VALID'):
+    x = tf.nn.conv2d(x, W, strides=strides, padding=pad_type)
+    x = tf.nn.bias_add(x, b)
+    return tf.nn.selu(x)
+
+
+
+
+def conv_net(x,A_ref,G_ref,T_ref,C_ref,weights, biases,keep):
+    conv1_1 = conv2d(x, weights['wc1_1'], biases['bc1_1'],strides=[1,1,1,1],pad_type='SAME')
+    conv1_2 = conv2d(x, weights['wc1_2'], biases['bc1_2'],strides=[1,1,1,1],pad_type='SAME')
+    conv1_3 = conv2d(x, weights['wc1_3'], biases['bc1_3'],strides=[1,1,1,1],pad_type='SAME')
+    
+    merge_conv1=tf.concat([conv1_1, conv1_2,conv1_3],3)
+    conv2 = conv2d(merge_conv1, weights['wc2'], biases['bc2'],strides=[1,1,2,1])
+    
+    #pooled=tf.nn.max_pool(merge_conv1,ksize=[1,5,1,1],strides=[1,1,1,1],padding='SAME')
+    #conv2 = conv2d(pooled, weights['wc2'], biases['bc2'],strides=[1,1,2,1])
+    
+    conv3 = conv2d(conv2, weights['wc3'], biases['bc3'],strides=[1,1,2,1])
+
+    flat_nn = tf.reshape(conv3, [-1, weights['wd1'].get_shape().as_list()[0]])
+    
+    drop_out_flat =flat_nn# tf.nn.dropout(flat_nn, keep)
+    
+    fc1 = tf.add(tf.matmul(drop_out_flat, weights['wd1']), biases['bd1'])
+    fc1 = tf.nn.selu(fc1)
+    
+    drop_out_fc1=tf.nn.dropout(fc1, keep)
+    
+    '''fc2 = tf.add(tf.matmul(drop_out_fc1, weights['wd2']), biases['bd2'])
+    fc2 = tf.nn.selu(fc2)
+    
+    drop_out_fc2=tf.nn.dropout(fc2, keep)'''
+    
+    fa = tf.add(tf.matmul(drop_out_fc1, weights['w-tr']), biases['w-tr'])
+    fa = tf.nn.selu(fa)
+    
+    out_A = tf.add(tf.matmul(tf.concat([fa,A_ref],1), weights['out-A']), biases['out-A'])
+    out_G = tf.add(tf.matmul(tf.concat([fa,G_ref],1), weights['out-G']), biases['out-G'])
+    out_T = tf.add(tf.matmul(tf.concat([fa,T_ref],1), weights['out-T']), biases['out-T'])
+    out_C = tf.add(tf.matmul(tf.concat([fa,C_ref],1), weights['out-C']), biases['out-C'])
+    
+    #fc2 = tf.add(tf.matmul(tf.concat([drop_out_fc1, out_A, out_T, out_G, out_C],1), weights['wd-gtp_0']), biases['bd-gtp_0'])
+    
+    fc2 = tf.add(tf.matmul(drop_out_fc1, weights['wd-gtp_0']), biases['bd-gtp_0'])
+    gt_fc1 = tf.nn.selu(fc2)
+    
+    
+    gt_fc2=tf.nn.selu(tf.add(tf.matmul(tf.concat([gt_fc1,out_A,out_G,out_T,out_C],1), weights['wd-gtp']), biases['bd-gtp']))
+    out_gt= tf.add(tf.matmul(gt_fc2, weights['out-gtp']), biases['out-gtp'])
+    
+    
+    
+    
+    return out_gt,out_A,out_G,out_T,out_C
+
+def get_tensors(n_input,learning_rate=0):
+    h_in,w_in,depth=n_input
+    
+    #h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    #w_in=int(np.ceil(float(w_in-4+1)/float(1)))
+    #h_in=28
+       
+    h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    w_in=int(np.ceil(float(w_in-3+1)/float(2)))
+
+    h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    w_in=int(np.ceil(float(w_in-3+1)/float(2)))
+    
+    a,b,c=16,16,16
+    d,e,f=32,64,48
+    q=32
+    g,h=16,16
+    j=8
+    i=0
+
+    weights = {
+        'wc1_1': tf.get_variable('W0_1', shape=(1,5,depth,a), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wc1_2': tf.get_variable('W0_2', shape=(5,1,depth,b), initializer=tf.contrib.layers.xavier_initializer()),
+        'wc1_3': tf.get_variable('W0_3', shape=(5,5,depth,c), initializer=tf.contrib.layers.xavier_initializer()),
+        'wc2': tf.get_variable('W1', shape=(2,3,a+b+c,d), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wc3': tf.get_variable('W2', shape=(2,3,d,e), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wd1': tf.get_variable('W3', shape=(h_in*w_in*e,f), initializer=tf.contrib.layers.xavier_initializer()),
+        #'wd2': tf.get_variable('W3_2', shape=(f,q), initializer=tf.contrib.layers.xavier_initializer()),
+        'w-tr': tf.get_variable('W4', shape=(f,g), initializer=tf.contrib.layers.xavier_initializer()),
+        'wd-gtp_0': tf.get_variable('W5', shape=(f,h), initializer=tf.contrib.layers.xavier_initializer()),
+        'wd-gtp': tf.get_variable('W6', shape=(8+h,j), initializer=tf.contrib.layers.xavier_initializer()), 
+        'out-gtp': tf.get_variable('W7', shape=(j,2), initializer=tf.contrib.layers.xavier_initializer()),
+        
+        'out-A': tf.get_variable('W8', shape=(g+1,2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-G': tf.get_variable('W9', shape=(g+1,2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-T': tf.get_variable('W10', shape=(g+1,2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-C': tf.get_variable('W11', shape=(g+1,2), initializer=tf.contrib.layers.xavier_initializer())
+    }
+    
+    biases = {
+        'bc1_1': tf.get_variable('B0_1', shape=(a), initializer=tf.contrib.layers.xavier_initializer()), 
+        'bc1_2': tf.get_variable('B0_2', shape=(b), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc1_3': tf.get_variable('B0_3', shape=(c), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc2': tf.get_variable('B1', shape=(d), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc3': tf.get_variable('B2', shape=(e), initializer=tf.contrib.layers.xavier_initializer()),
+        'bd1': tf.get_variable('B3', shape=(f), initializer=tf.contrib.layers.xavier_initializer()),
+        #'bd2': tf.get_variable('B3_2', shape=(q), initializer=tf.contrib.layers.xavier_initializer()),
+        'w-tr': tf.get_variable('B4', shape=(g), initializer=tf.contrib.layers.xavier_initializer()),
+        'bd-gtp_0': tf.get_variable('B5', shape=(h), initializer=tf.contrib.layers.xavier_initializer()),
+        'bd-gtp': tf.get_variable('B6', shape=(j), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-gtp': tf.get_variable('B7', shape=(2), initializer=tf.contrib.layers.xavier_initializer()),
+        
+        'out-A': tf.get_variable('B8', shape=(2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-G': tf.get_variable('B9', shape=(2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-T': tf.get_variable('B10', shape=(2), initializer=tf.contrib.layers.xavier_initializer()),
+        'out-C': tf.get_variable('B11', shape=(2), initializer=tf.contrib.layers.xavier_initializer()),
+    }
+    
+    
+    x,GT_label = tf.placeholder("float", [None]+n_input), tf.placeholder("float", [None, 2])
+    
+    A_label,G_label,T_label,C_label=tf.placeholder("float", [None, 2]), tf.placeholder("float", [None, 2]), tf.placeholder("float", [None, 2]), tf.placeholder("float", [None, 2])
+    
+    A_ref,G_ref,T_ref,C_ref=tf.placeholder("float", [None, 1]), tf.placeholder("float", [None, 1]), tf.placeholder("float", [None, 1]), tf.placeholder("float", [None, 1])
+    
+    keep = tf.placeholder(tf.float32)
+    
+    GT_score,A_score,G_score,T_score,C_score = conv_net(x,A_ref,G_ref,T_ref,C_ref, weights, biases, keep)
+
+    gamma=0.001
+    reg_loss=tf.add_n([tf.nn.l2_loss(t) for t in weights.values()])
+    
+    prob_GT=tf.nn.softmax(logits=GT_score)
+    prob_A=tf.nn.softmax(logits=A_score)
+    prob_G=tf.nn.softmax(logits=G_score)
+    prob_T=tf.nn.softmax(logits=T_score)
+    prob_C=tf.nn.softmax(logits=C_score)
+    
+    cost_GT=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=GT_score, labels=GT_label))
+    cost_A=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=A_score, labels=A_label))
+    cost_G=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=G_score, labels=G_label))
+    cost_T=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=T_score, labels=T_label))
+    cost_C=tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=C_score, labels=C_label))
+    
+    cost = cost_A+cost_G+cost_T+cost_C+tf.constant(1.0)*cost_GT+tf.constant(gamma)*tf.reduce_mean(reg_loss)
+    
+    prediction_GT=tf.argmax(prob_GT, 1)
+    prediction_A=tf.argmax(prob_A, 1)
+    prediction_G=tf.argmax(prob_G, 1)
+    prediction_T=tf.argmax(prob_T, 1)
+    prediction_C=tf.argmax(prob_C, 1)
+    
+    
+    optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(cost)
+    
+
+    
+    prediction_accuracy_A= tf.equal(prediction_A, tf.argmax(A_label, 1))
+    prediction_accuracy_G= tf.equal(prediction_G, tf.argmax(G_label, 1))
+    prediction_accuracy_T= tf.equal(prediction_T, tf.argmax(T_label, 1))
+    prediction_accuracy_C= tf.equal(prediction_C, tf.argmax(C_label, 1))
+    prediction_accuracy_GT=tf.equal(prediction_GT, tf.argmax(GT_label, 1))
+    
+    accuracy_GT = tf.reduce_sum(tf.cast(prediction_accuracy_GT, tf.float32))
+    accuracy_A = tf.reduce_sum(tf.cast(prediction_accuracy_A, tf.float32))
+    accuracy_G = tf.reduce_sum(tf.cast(prediction_accuracy_G, tf.float32))
+    accuracy_T = tf.reduce_sum(tf.cast(prediction_accuracy_T, tf.float32))
+    accuracy_C = tf.reduce_sum(tf.cast(prediction_accuracy_C, tf.float32))
+    
+    
+    
+    tmp1=tf.logical_and(prediction_accuracy_GT, prediction_accuracy_A)
+    tmp2=tf.logical_and(prediction_accuracy_G, prediction_accuracy_T)
+    tmp3=tf.logical_and(tmp1,tmp2)
+    accuracy=tf.reduce_sum(tf.cast(tf.logical_and(prediction_accuracy_C,tmp3), tf.float32))
+    
+
+    tensors=(x,GT_label,A_label, G_label, T_label, C_label,GT_score, A_score, G_score, T_score, C_score, accuracy_GT, accuracy_A,  accuracy_G,  accuracy_T,  accuracy_C, prediction_accuracy_GT, prediction_accuracy_A,  prediction_accuracy_G,  prediction_accuracy_T,  prediction_accuracy_C, prediction_GT, prediction_A,  prediction_G,  prediction_T,  prediction_C, accuracy, cost, optimizer,  cost_GT, cost_A, cost_G, cost_T, cost_C,A_ref,G_ref,T_ref,C_ref,prob_GT,prob_A,prob_G,prob_T,prob_C,keep)
+
+    
+    
+    return weights,biases,tensors
diff --git a/misc/training/model_architect_indel.py b/misc/training/model_architect_indel.py
new file mode 100644
index 0000000..09be110
--- /dev/null
+++ b/misc/training/model_architect_indel.py
@@ -0,0 +1,180 @@
+import tensorflow as tf
+import numpy as np
+        
+def conv2d(x, W, b, strides=1,pad_type='VALID'):
+    x = tf.nn.conv2d(x, W, strides=strides, padding=pad_type)
+    x = tf.nn.bias_add(x, b)
+    return tf.nn.selu(x)
+
+
+
+
+def conv_net(x,weights, biases,rate):
+    conv1_1 = conv2d(x, weights['wc1_1'], biases['bc1_1'],strides=[1,1,1,1],pad_type='SAME')
+    conv1_2 = conv2d(x, weights['wc1_2'], biases['bc1_2'],strides=[1,1,1,1],pad_type='SAME')
+    conv1_3 = conv2d(x, weights['wc1_3'], biases['bc1_3'],strides=[1,1,1,1],pad_type='SAME')
+    
+    merge_conv1=tf.concat([conv1_1, conv1_2,conv1_3],3)
+    conv2 = conv2d(merge_conv1, weights['wc2'], biases['bc2'],strides=[1,1,2,1])
+    
+    #pooled=tf.nn.max_pool(merge_conv1,ksize=[1,5,1,1],strides=[1,1,1,1],padding='SAME')
+    #conv2 = conv2d(pooled, weights['wc2'], biases['bc2'],strides=[1,1,2,1])
+    
+    conv3 = conv2d(conv2, weights['wc3'], biases['bc3'],strides=[1,1,2,1])
+
+    flat_nn = tf.reshape(conv3, [-1, weights['wd1'].get_shape().as_list()[0]])
+
+    fc1 = tf.add(tf.matmul(flat_nn, weights['wd1']), biases['bd1'])
+    fc1 = tf.nn.selu(fc1)
+    drop_out_fc1=tf.nn.dropout(fc1, rate=rate)
+    
+    '''fc2 = tf.add(tf.matmul(drop_out_fc1, weights['wd2']), biases['bd2'])
+    fc2 = tf.nn.selu(fc2)
+    drop_out_fc2=tf.nn.dropout(fc2, rate)'''
+    
+    fc3 = tf.add(tf.matmul(drop_out_fc1, weights['wd3']), biases['bd3'])
+    fc3 = tf.nn.selu(fc3)
+    
+    out = tf.add(tf.matmul(fc3, weights['out']), biases['out'])
+    
+    return out
+
+def get_tensors(n_input,learning_rate=0):
+    n_input=[n_input[0]*3,n_input[1],n_input[2]]
+    h_in,w_in,depth=n_input
+    #h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    #w_in=int(np.ceil(float(w_in-4+1)/float(1)))
+    #h_in=28
+       
+    h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    w_in=int(np.ceil(float(w_in-3+1)/float(2)))
+
+    h_in=int(np.ceil(float(h_in-2+1)/float(1)))
+    w_in=int(np.ceil(float(w_in-3+1)/float(2)))
+    
+    a,b,c=8,8,8
+    d,e,f=32,48,32
+    q=32
+    g,h=24,16
+    j=8
+    i=0
+    
+    pi=0.01
+    bias=-np.log((1-pi)/pi)
+    
+    weights = {
+        'wc1_1': tf.compat.v1.get_variable('W0_1', shape=(1,5,depth,a), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wc1_2': tf.compat.v1.get_variable('W0_2', shape=(5,1,depth,b), initializer=tf.contrib.layers.xavier_initializer()),
+        'wc1_3': tf.compat.v1.get_variable('W0_3', shape=(5,5,depth,c), initializer=tf.contrib.layers.xavier_initializer()),
+        'wc2': tf.compat.v1.get_variable('W1', shape=(2,3,a+b+c,d), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wc3': tf.compat.v1.get_variable('W2', shape=(2,3,d,e), initializer=tf.contrib.layers.xavier_initializer()), 
+        'wd1': tf.compat.v1.get_variable('W3', shape=(h_in*w_in*e,f), initializer=tf.contrib.layers.xavier_initializer()),
+        #'wd2': tf.get_variable('W3_2', shape=(f,q), initializer=tf.contrib.layers.xavier_initializer()),
+        'wd3': tf.compat.v1.get_variable('W3_3', shape=(f,g), initializer=tf.contrib.layers.xavier_initializer()),
+      
+        'out': tf.compat.v1.get_variable('W13', shape=(g,4), initializer=tf.contrib.layers.xavier_initializer())
+    }
+    
+    biases = {
+        'bc1_1': tf.compat.v1.get_variable('B0_1', shape=(a), initializer=tf.contrib.layers.xavier_initializer()), 
+        'bc1_2': tf.compat.v1.get_variable('B0_2', shape=(b), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc1_3': tf.compat.v1.get_variable('B0_3', shape=(c), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc2': tf.compat.v1.get_variable('B1', shape=(d), initializer=tf.contrib.layers.xavier_initializer()),
+        'bc3': tf.compat.v1.get_variable('B2', shape=(e), initializer=tf.contrib.layers.xavier_initializer()),
+        'bd1': tf.compat.v1.get_variable('B3', shape=(f), initializer=tf.contrib.layers.xavier_initializer()),
+        #'bd2': tf.get_variable('B3_2', shape=(q), initializer=tf.contrib.layers.xavier_initializer()),
+        'bd3': tf.compat.v1.get_variable('B3_3', shape=(g), initializer=tf.contrib.layers.xavier_initializer()),
+        'out': tf.compat.v1.get_variable('B13', shape=(4), initializer=tf.contrib.layers.xavier_initializer()),}
+    #tf.constant([1, bias,bias,bias])}
+    
+    
+    
+    x0,x1,x2= tf.compat.v1.placeholder("float", [None]+n_input), tf.compat.v1.placeholder("float", [None]+n_input), tf.compat.v1.placeholder("float", [None]+n_input),
+    
+    gt=tf.compat.v1.placeholder("float", [None, 4])
+    
+    
+    rate = tf.compat.v1.placeholder(tf.float32)
+
+    score = conv_net(x2, weights, biases, rate)
+
+    gamma=1e-5
+    alpha=tf.constant([0.1, 0.3, 0.2, 0.4])
+    
+    
+    
+    reg_loss=tf.add_n([tf.nn.l2_loss(t) for t in weights.values()])
+    #ref_loss=tf.add_n([tf.reduce_sum(tf.abs(t)) for t in weights.values()])
+    
+    prob=tf.nn.softmax(logits=score)
+    prediction=tf.argmax(prob, 1)
+    
+    epsilon = 1.e-9
+    decay=2.0
+    
+    y_true = tf.convert_to_tensor(gt, tf.float32)
+
+    model_out = tf.add(prob, epsilon)
+    ce = tf.multiply(y_true, -tf.log(model_out))
+    weight = tf.multiply(y_true, tf.pow(tf.subtract(1., model_out), decay))
+    fl = tf.multiply(alpha, tf.multiply(weight, ce))
+    reduced_fl = tf.reduce_max(fl, axis=1)
+
+    fl_cost=tf.reduce_mean(reduced_fl)
+    #cost = fl_cost+tf.constant(gamma)*reg_loss
+    
+    cost = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=score, labels=gt))+tf.constant(gamma)*reg_loss
+    
+    
+    optimizer = tf.compat.v1.train.AdamOptimizer(learning_rate=learning_rate).minimize(cost)
+    
+
+    
+    prediction_accuracy= tf.equal(prediction, tf.argmax(gt, 1))
+    
+
+   
+    accuracy=tf.reduce_sum(tf.cast(prediction_accuracy, tf.float32))
+    
+    
+
+    tensors=(x0, x1,x2,gt, accuracy, cost, optimizer, prob, rate)
+    
+    
+    '''score_0 = conv_net(x0, weights, biases, rate)
+    score_1 = conv_net(x1, weights, biases, rate)
+
+    gamma=0.0001
+    reg_loss=tf.add_n([tf.nn.l2_loss(t) for t in weights.values()])
+    
+    prob_0=tf.nn.softmax(logits=score_0)
+    prob_1=tf.nn.softmax(logits=score_1)
+
+    prediction_0=tf.argmax(prob_0, 1)
+    prediction_1=tf.argmax(prob_1, 1)
+
+    
+    cost = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=score_0, labels=label_1))+tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=score_1, labels=label_2))+ tf.constant(gamma)*tf.reduce_mean(reg_loss)
+    
+    
+    optimizer = tf.compat.v1.train.AdamOptimizer(learning_rate=learning_rate).minimize(cost)
+    
+
+    
+    prediction_0_accuracy= tf.equal(prediction_0, tf.argmax(label_1, 1))
+    prediction_1_accuracy= tf.equal(prediction_1, tf.argmax(label_2, 1))
+    
+  
+    
+    accuracy_0 = tf.reduce_sum(tf.cast(prediction_0_accuracy, tf.float32))
+    accuracy_1 = tf.reduce_sum(tf.cast(prediction_1_accuracy, tf.float32))
+   
+    accuracy=tf.reduce_sum(tf.cast(tf.logical_and(prediction_0_accuracy,prediction_1_accuracy), tf.float32))
+    
+    
+
+    tensors=(x0, x1, label_1,label_2, accuracy, accuracy_0, accuracy_1, cost, optimizer, prob_0,prob_1, rate)'''
+
+    
+    
+    return weights,biases,tensors
diff --git a/misc/training/model_run.py b/misc/training/model_run.py
new file mode 100644
index 0000000..e191bfe
--- /dev/null
+++ b/misc/training/model_run.py
@@ -0,0 +1,290 @@
+from warnings import simplefilter 
+simplefilter(action='ignore', category=FutureWarning)
+
+import time, os, copy, argparse, subprocess, glob, re, datetime
+os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' 
+
+import numpy as np
+import multiprocessing as mp
+import tensorflow as tf
+from model_architect import *
+from utils import *
+from scipy.sparse import csr_matrix
+
+if type(tf.contrib) != type(tf): tf.contrib._warning = None
+
+
+config = tf.ConfigProto(device_count={"CPU": 32})
+#config = tf.ConfigProto()
+config.gpu_options.allow_growth = True
+
+num_to_base_map={0:'A',1:'G',2:'T',3:'C'}
+
+class Genome:
+    def __init__(self, name, path):
+        self.name = name
+        self.path = path
+        self.training_data = {}
+        self.testing_data = None
+        self.chunk_list=[]
+
+    def add_data(self, pool, sparse, val):
+        print('%s: Starting reading genome %s.' %(str(datetime.datetime.now()), self.name), flush=True)
+        train_list=glob.glob(os.path.join(self.path,'train*'))
+        chunk_list=[re.findall('pileups.(\d+)', x)[0] for x in train_list]
+        
+        print('\nProgress\n\nTotal: %s\nDone : ' %('.'*len(chunk_list)),end='', flush=True)
+        
+        for chunk, path in zip(chunk_list, train_list):
+            self.chunk_list.append(chunk)
+            x_train,y_train,train_allele,train_ref= get_data(path,pool)
+
+            if sparse:
+                self.training_data[chunk]=(csr_matrix(x_train.reshape(-1)),y_train,train_allele,train_ref)
+            else:
+                self.training_data[chunk]=(x_train,y_train,train_allele,train_ref)
+                
+            print('.',end='',flush=True)
+        
+        if val:
+            self.testing_data = get_data(os.path.join(self.path,'test.pileups'),pool)
+            
+        print('\n\n%s: Finished reading genome %s.\n' %(str(datetime.datetime.now()), self.name), flush=True)
+    
+def train_SNP_model(params):
+    
+    sparse=params['sparse']
+    
+    tf.reset_default_graph()
+    cpu=params['cpu']
+    params['val']=True
+    dims=[5,41,5]
+    n_input=dims
+    
+    pool = mp.Pool(processes=cpu)
+    
+    print('%s: Starting reading pileups.' %str(datetime.datetime.now()),flush=True)
+
+    genomes_list=[]
+    with open(params['train_path'],'r') as file:
+        for line in file:
+            x=line.rstrip('\n').split(',')
+            
+            current_genome=Genome(*x)
+            current_genome.add_data(pool,sparse, params['val'])
+            
+            genomes_list.append(current_genome)
+            
+    pool.close()  
+    pool.join()
+
+    print('\n%s: Finished reading pileups.' %str(datetime.datetime.now()),flush=True)
+    
+    training_iters, learning_rate, batch_size= params['iters'],\
+    params['rate'], params['size']
+
+    weights,biases,tensors=get_tensors(n_input,learning_rate)
+    (x,GT_label,A_label, G_label, T_label, C_label,GT_score, A_score, G_score, T_score, C_score, accuracy_GT, accuracy_A,  accuracy_G,  accuracy_T,  accuracy_C, prediction_accuracy_GT, prediction_accuracy_A,  prediction_accuracy_G,  prediction_accuracy_T,  prediction_accuracy_C, prediction_GT, prediction_A,  prediction_G,  prediction_T,  prediction_C, accuracy, cost, optimizer,  cost_GT, cost_A, cost_G, cost_T, cost_C,A_ref,G_ref,T_ref,C_ref,prob_GT,prob_A,prob_G,prob_T,prob_C,keep)=tensors
+
+    init = tf.global_variables_initializer()
+    saver = tf.train.Saver(max_to_keep=10000)
+    
+    
+    with tf.Session(config=config)  as sess:
+        sess.run(init)
+        sess.run(tf.local_variables_initializer())
+        if params['rt_path']:
+            print('%s: Retraining model:> %s' %(str(datetime.datetime.now()),params['rt_path']))
+            saver.restore(sess, params['rt_path'])        
+
+        stats,v_stats=[],[]
+        
+        count=0
+        
+        save_num=1        
+        
+        for k in range(training_iters):
+            t=time.time()
+            print('\n'+100*'-'+'\n',flush=True)
+            print('%s: Starting epoch #: %d\n' %(str(datetime.datetime.now()),k),flush=True)
+            
+            for genome in genomes_list:
+                print('%s: Training on genome %s \n' %(str(datetime.datetime.now()), genome.name), flush=True)
+                
+                chunk_list=genome.chunk_list
+                print('Progress\n\nTotal: %s\nDone : ' %('.'*len(chunk_list)),end='', flush=True)
+                
+                for chunk in chunk_list:
+                    training_loss, total_train_data, train_acc = 0, 0, 0
+                    
+                    x_train_sparse,y_train,train_allele,train_ref=genome.training_data[chunk]
+
+                    if sparse:
+                        x_train_flat=np.array(x_train_sparse.todense())
+                        x_train=x_train_flat.reshape(x_train_flat.shape[1]//1025,5,41,5)
+                    else:
+                        x_train=x_train_sparse
+
+                    for batch in range(len(x_train)//batch_size):
+                        batch_x = x_train[batch*batch_size:min((batch+1)*batch_size, len(x_train))]
+
+                        batch_y = y_train[batch*batch_size:min((batch+1)*batch_size, len(y_train))]
+
+                        batch_ref = train_ref[batch*batch_size :min((batch+1)*batch_size, len(train_ref))]
+
+                        batch_allele = train_allele[batch*batch_size :min((batch+1)*batch_size, len(train_allele))]
+
+                        opt,loss,batch_acc = sess.run([optimizer,cost,accuracy], feed_dict={x: batch_x, GT_label:batch_y, A_label:np.eye(2)[batch_allele[:,0]], G_label:np.eye(2)[batch_allele[:,1]], T_label:np.eye(2)[batch_allele[:,2]], C_label:np.eye(2)[batch_allele[:,3]] , A_ref:batch_ref[:,0][:,np.newaxis], G_ref:batch_ref[:,1][:,np.newaxis], T_ref:batch_ref[:,2][:,np.newaxis], C_ref:batch_ref[:,3][:,np.newaxis], keep:0.5})
+
+                        training_loss+=loss*len(batch_x)
+                        total_train_data+=len(batch_x)
+                        train_acc+=batch_acc
+
+                    
+
+                    print('.',end='',flush=True)
+
+                    training_loss=training_loss/total_train_data
+                    train_acc=train_acc/total_train_data
+                    
+                print('\n\nGenome: %s      Training Loss: %.4f      Training Accuracy: %.4f\n' %(genome.name, training_loss, train_acc), flush=True)
+                
+                if params['val']:
+                    print(50*'*'+'\n', flush=True)
+                    print('%s: Performing Validation\n' %str(datetime.datetime.now()), flush=True)
+                    
+                    for val_genome in genomes_list:
+                        vx_test, vy_test, vtest_allele, vtest_ref = val_genome.testing_data
+                        
+                        v_loss,v_acc,A_acc,G_acc,T_acc,C_acc,GT_acc,v_loss,v_acc=0,0,0,0,0,0,0,0,0
+
+                        for batch in range(int(np.ceil(len(vx_test)/(batch_size)))):
+                            vbatch_x = vx_test[batch*batch_size:min((batch+1)*batch_size,len(vx_test))]
+                            vbatch_y = vy_test[batch*batch_size:min((batch+1)*batch_size,len(vx_test))] 
+                            vbatch_allele = vtest_allele[batch*batch_size:min((batch+1)*batch_size,len(vx_test))]
+                            vbatch_ref = vtest_ref[batch*batch_size:min((batch+1)*batch_size,len(vx_test))]
+
+                            batch_loss,batch_acc, batch_GT_acc, batch_A_acc, batch_G_acc, batch_T_acc, batch_C_acc,batch_prediction_GT, batch_prediction_A,  batch_prediction_G,  batch_prediction_T,  batch_prediction_C,bGT_score,bA_score,bG_score,bT_score,bC_score = sess.run([cost, accuracy, accuracy_GT, accuracy_A,  accuracy_G,  accuracy_T,  accuracy_C,prediction_GT, prediction_A,  prediction_G,  prediction_T,  prediction_C,GT_score,A_score,G_score,T_score,C_score], feed_dict={x: vbatch_x,GT_label:vbatch_y, A_label:np.eye(2)[vbatch_allele[:,0]], G_label:np.eye(2)[vbatch_allele[:,1]], T_label:np.eye(2)[vbatch_allele[:,2]], C_label:np.eye(2)[vbatch_allele[:,3]],  A_ref:vbatch_ref[:,0][:,np.newaxis], G_ref:vbatch_ref[:,1][:,np.newaxis], T_ref:vbatch_ref[:,2][:,np.newaxis], C_ref:vbatch_ref[:,3][:,np.newaxis], keep:1.0})
+
+                            v_loss+=batch_loss*len(vbatch_x)
+                            v_acc+=batch_acc
+                            A_acc+=batch_A_acc
+                            G_acc+=batch_G_acc
+                            T_acc+=batch_T_acc
+                            C_acc+=batch_C_acc
+                            GT_acc+=batch_GT_acc
+
+                        print('Genome: %s    Validation Loss: %.4f    Validation Accuracy: %.4f' %(val_genome.name, v_loss/len(vx_test), v_acc/len(vx_test)), flush=True)
+
+                        print('Genome: %s    GT_acc=%.4f    A_acc=%.4f    G_acc=%.4f    T_acc=%.4f    C_acc=%.4f\n' %(val_genome.name, GT_acc/len(vx_test), A_acc/len(vx_test), G_acc/len(vx_test), T_acc/len(vx_test), C_acc/len(vx_test)), flush=True)
+
+                print(50*'*'+'\n', flush=True)
+            saver.save(sess, save_path=os.path.join(params['model'],'model'),global_step=save_num)
+            elapsed=time.time()-t
+            save_num+=1
+            
+            print ('%s: Time Taken for Iteration %d: %.2f seconds\n' %(str(datetime.datetime.now()),k,elapsed), flush=True)
+            
+            
+
+            
+def get_data(fname,pool,a=None, b=None,dims=(5,41,5)):
+    t=time.time()
+    l=os.stat(fname).st_size
+    
+    rec_size=15+dims[0]*dims[1]*dims[2]*6
+    if a!=None and b!=None:
+        my_array=[(fname,x,dims) for x in range(a,b,1000*rec_size)]
+    else:
+        my_array=[(fname,x,dims) for x in range(0,l,1000*rec_size)]
+
+    results = pool.map(read_pileups_from_file, my_array)
+    
+    pos=np.vstack([res[0][:,np.newaxis] for res in results])
+    mat=np.vstack([res[1] for res in results])
+    
+
+    ref=np.vstack([res[2] for res in results])
+    allele, gt=None, None
+
+    allele=np.vstack([res[3] for res in results])
+    
+    gt=np.vstack([res[4] for res in results])
+   
+
+    return mat,gt,allele,ref
+    
+def read_pileups_from_file(options):
+    
+    fname,n,dims=options
+    file= open(fname,'r')
+    file.seek(n)
+    mat=[]
+    pos=[]
+    ref=[]
+    allele1,allele2=[],[]
+    gt=[]
+    
+    dp=[]
+    freq=[]
+    
+    i=0
+    while i<1000:
+        i+=1
+        c=file.read(15)
+        if not c:
+            break
+        pos.append(int(c[:11]))
+        gt.append(int(c[11]))
+        allele1.append(int(c[12]))
+        allele2.append(int(c[13]))
+        ref.append(int(c[14]))
+
+
+
+        m=file.read(dims[0]*dims[1]*dims[2]*6)
+        p_mat=np.array([int(m[6*i:6*i+6]) for i in range(dims[0]*dims[1]*dims[2])]).reshape((dims[0],dims[1],dims[2]))
+
+        mat.append(p_mat)
+
+    mat=np.array(mat)    
+    pos=np.array(pos)
+    ref=np.eye(4)[np.array(ref)].astype(np.int8)
+    allele1=np.eye(4)[np.array(allele1)].astype(np.int8)
+    allele2=np.eye(4)[np.array(allele2)].astype(np.int8)
+    
+    allele=allele1+allele2
+    
+    allele[allele==2]=1
+    
+    gt=np.eye(2)[np.array(gt)].astype(np.int8)
+
+    return (pos,mat.astype(np.int16),ref,allele,gt)
+
+if __name__ == '__main__':
+    t=time.time()
+    
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-r", "--rate", help="Learning rate",type=float)
+    parser.add_argument("-i", "--iterations", help="Training iterations",type=int)
+    parser.add_argument("-s", "--size", help="Batch size",type=int)
+    parser.add_argument("-train", "--train", help="Train path")
+    parser.add_argument("-model", "--model", help="Model output path")
+    parser.add_argument("-cpu", "--cpu", help="CPUs",type=int)
+    parser.add_argument("-rt_path", "--retrain_path", help="Retrain saved model",type=str)
+    parser.add_argument('-sparse','--sparse', help='Stores features as sparse matrices', default=False, action='store_true')
+    parser.add_argument('-val','--validation', help='Perform validation', default=True, action='store_false')
+    
+    args = parser.parse_args()
+
+    
+    os.makedirs(args.model, exist_ok=True)
+    
+    in_dict={'cpu':args.cpu,'rate':args.rate, 'iters':args.iterations, 'size':args.size,\
+             'train_path':args.train, 'model':args.model, 'val':args.validation,'rt_path':args.retrain_path,\
+            'sparse':args.sparse}
+    
+    train_SNP_model(in_dict)
+
+    elapsed=time.time()-t
+    print ('Total Time Elapsed: %.2f seconds' %elapsed)
\ No newline at end of file
diff --git a/misc/training/model_run_indels.py b/misc/training/model_run_indels.py
new file mode 100644
index 0000000..a7f2233
--- /dev/null
+++ b/misc/training/model_run_indels.py
@@ -0,0 +1,260 @@
+from warnings import simplefilter 
+simplefilter(action='ignore', category=FutureWarning)
+
+import time, os, copy, argparse, subprocess, glob, re, datetime
+os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' 
+
+import numpy as np
+import multiprocessing as mp
+import tensorflow as tf
+from model_architect_indel import *
+
+if type(tf.contrib) != type(tf): tf.contrib._warning = None
+
+
+config = tf.ConfigProto(device_count={"CPU": 32})
+#config = tf.ConfigProto()
+config.gpu_options.allow_growth = True
+
+
+rev_gt_map={0:'hom-ref', 1:'hom-alt', 2:'het-ref', 3:'het-alt'}
+rev_base_map={0:'A',1:'G',2:'T',3:'C',4:'-'}
+
+class Genome:
+    def __init__(self, name, path):
+        self.name = name
+        self.path = path
+        self.training_data = {}
+        self.testing_data = None
+        self.chunk_list=[]
+
+    def add_data(self, pool, val):
+        print('%s: Starting reading genome %s.' %(str(datetime.datetime.now()), self.name), flush=True)
+        train_list=glob.glob(os.path.join(self.path,'train*'))
+        chunk_list=[re.findall('pileups.(\d+)', x)[0] for x in train_list]
+        
+        print('\nProgress\n\nTotal: %s\nDone : ' %('.'*len(chunk_list)),end='', flush=True)
+        
+        for chunk, path in zip(chunk_list, train_list):
+            self.chunk_list.append(chunk)
+            x_train, train_gt= get_data(path,pool)
+
+            self.training_data[chunk]=(x_train, train_gt)
+                
+            print('.',end='',flush=True)
+        
+        if val:
+            self.testing_data = get_data(os.path.join(self.path,'test.pileups'),pool)
+            
+        print('\n\n%s: Finished reading genome %s.\n' %(str(datetime.datetime.now()), self.name), flush=True)
+    
+def train_indel_model(params):
+    
+    
+    tf.reset_default_graph()
+    cpu=params['cpu']
+    params['val']=True
+    dims=[5,41,5]
+    n_input=dims
+    
+    pool = mp.Pool(processes=cpu)
+    
+    print('%s: Starting reading pileups.' %str(datetime.datetime.now()),flush=True)
+
+    genomes_list=[]
+    with open(params['train_path'],'r') as file:
+        for line in file:
+            x=line.rstrip('\n').split(',')
+            
+            current_genome=Genome(*x)
+            current_genome.add_data(pool, params['val'])
+            
+            genomes_list.append(current_genome)
+            
+    pool.close()  
+    pool.join()
+
+    print('\n%s: Finished reading pileups.' %str(datetime.datetime.now()),flush=True)
+    
+    training_iters, learning_rate, batch_size= params['iters'],\
+    params['rate'], params['size']
+
+    training_iters, learning_rate, batch_size= params['iters'],\
+    params['rate'], params['size']
+
+    weights,biases,tensors=get_tensors([5,128,2],learning_rate)
+    (x0, x1,x2,gt, accuracy, cost, optimizer, prob, rate)=tensors
+
+    
+    init = tf.global_variables_initializer()
+    saver = tf.train.Saver(max_to_keep=10000)
+    
+    
+    with tf.Session(config=config)  as sess:
+        sess.run(init)
+        sess.run(tf.local_variables_initializer())
+        if params['rt_path']:
+            print('%s: Retraining model:> %s' %(str(datetime.datetime.now()),params['rt_path']))
+            saver.restore(sess, params['rt_path'])        
+
+        stats,v_stats=[],[]
+        
+        count=0
+        
+        save_num=1        
+        
+        for k in range(training_iters):
+            t=time.time()
+            print('\n'+100*'-'+'\n',flush=True)
+            print('%s: Starting epoch #: %d\n' %(str(datetime.datetime.now()),k),flush=True)
+            
+            for genome in genomes_list:
+                print('%s: Training on genome %s \n' %(str(datetime.datetime.now()), genome.name), flush=True)
+                
+                chunk_list=genome.chunk_list
+                print('Progress\n\nTotal: %s\nDone : ' %('.'*len(chunk_list)),end='', flush=True)
+                
+                for chunk in chunk_list:
+                    training_loss, total_train_data, train_acc = 0, 0, 0
+                    
+                    x_train, train_gt=genome.training_data[chunk]
+
+                    for batch in range(len(x_train)//batch_size):
+                        batch_x = x_train[batch*batch_size:min((batch+1)*batch_size, len(x_train))]
+                        batch_gt = train_gt[batch*batch_size:min((batch+1)*batch_size, len(train_gt))]
+
+                        opt, loss, batch_acc = sess.run([optimizer,cost,accuracy], feed_dict={x2: batch_x, gt:batch_gt, rate:0.5})
+
+                        training_loss+=loss*len(batch_x)
+                        total_train_data+=len(batch_x)
+                        train_acc+=batch_acc
+
+                    
+
+                    print('.',end='',flush=True)
+
+                    training_loss=training_loss/total_train_data
+                    train_acc=train_acc/total_train_data
+                    
+                print('\n\nGenome: %s      Training Loss: %.4f      Training Accuracy: %.4f\n' %(genome.name, training_loss, train_acc), flush=True)
+                
+                if params['val']:
+                    print(50*'*'+'\n', flush=True)
+                    print('%s: Performing Validation\n' %str(datetime.datetime.now()), flush=True)
+                    
+                    for val_genome in genomes_list:
+                        vx_test, vtest_gt= val_genome.testing_data
+                        
+                        v_loss,v_acc,total_test_data=0,0,0
+
+                        for batch in range(int(np.ceil(len(vx_test)/(batch_size)))):
+                            vbatch_x = vx_test[batch*batch_size:min((batch+1)*batch_size,len(vx_test))]
+                            vbatch_gt = vtest_gt[batch*batch_size:min((batch+1)*batch_size,len(vtest_gt))] 
+
+                            batch_loss,batch_acc= sess.run([cost, accuracy], feed_dict={x2: vbatch_x, gt:vbatch_gt, rate:0.0})
+                            v_loss+=batch_loss*len(vbatch_x)
+                            v_acc+=batch_acc
+                            total_test_data+=len(vbatch_x)
+
+                        print('Genome: %s    Validation Loss: %.4f    Validation Accuracy: %.4f' %(val_genome.name, v_loss/total_test_data, v_acc/total_test_data), flush=True)
+
+                print(50*'*'+'\n', flush=True)
+            saver.save(sess, save_path=os.path.join(params['model'],'model'),global_step=save_num)
+            elapsed=time.time()-t
+            save_num+=1
+            
+            print ('%s: Time Taken for Iteration %d: %.2f seconds\n' %(str(datetime.datetime.now()),k,elapsed), flush=True)
+            
+            
+
+            
+def get_data(fname,pool):
+    t=time.time()
+    l=os.stat(fname).st_size
+    dims=(5,128,2)
+    
+    rec_size=12+dims[0]*dims[1]*dims[2]*4*3
+    my_array=[(fname,x,'train',dims) for x in range(0,l,1000*rec_size)]
+    
+    results = pool.map(read_pileups_from_file, my_array)
+
+    mat=np.vstack([res[0] for res in results])
+    gt=np.vstack([res[1] for res in results])
+    
+    return mat,gt
+
+def read_pileups_from_file(options):    
+    fname,n,mode,dims=options
+    file= open(fname,'r')
+    file.seek(n)
+    
+    mat_0,mat_1,mat_2=[],[],[]
+    pos=[]
+    gt=[]
+    
+    i=0
+
+    while i<1000:
+        i+=1
+        c=file.read(12)
+        if not c:
+            break             
+
+        gt.append(int(c[11]))
+
+        m=file.read(dims[0]*dims[1]*dims[2]*4)
+        p_mat_0=np.array([int(float(m[4*i:4*i+4])) for i in range(dims[0]*dims[1]*dims[2])]).reshape((dims[0],dims[1],dims[2]))
+        mat_0.append(p_mat_0)
+
+        m=file.read(dims[0]*dims[1]*dims[2]*4)
+        p_mat_1=np.array([int(float(m[4*i:4*i+4])) for i in range(dims[0]*dims[1]*dims[2])]).reshape((dims[0],dims[1],dims[2]))
+        mat_1.append(p_mat_1)
+
+        m=file.read(dims[0]*dims[1]*dims[2]*4)
+        p_mat_2=np.array([int(float(m[4*i:4*i+4])) for i in range(dims[0]*dims[1]*dims[2])]).reshape((dims[0],dims[1],dims[2]))
+        mat_2.append(p_mat_2)
+
+    mat_0=norm(np.array(mat_0))
+    mat_1=norm(np.array(mat_1))
+    mat_2=norm(np.array(mat_2))
+
+    mat=np.hstack([mat_0, mat_1, mat_2])
+
+    gt=np.array(gt)
+    gt=np.eye(4)[gt].astype(np.int8)
+
+    return (mat,gt)
+
+
+def norm(batch_x0):
+    batch_x0=batch_x0.astype(np.float32)
+    batch_x0[:,:,:,0]=batch_x0[:,:,:,0]/(np.sum(batch_x0[:,:,:,0],axis=1)[:,np.newaxis,:])-batch_x0[:,:,:,1]
+    return batch_x0
+
+if __name__ == '__main__':
+    t=time.time()
+    
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-r", "--rate", help="Learning rate",type=float)
+    parser.add_argument("-i", "--iterations", help="Training iterations",type=int)
+    parser.add_argument("-s", "--size", help="Batch size",type=int)
+    parser.add_argument("-train", "--train", help="Train path")
+    parser.add_argument("-model", "--model", help="Model output path")
+    parser.add_argument("-cpu", "--cpu", help="CPUs",type=int)
+    parser.add_argument("-rt_path", "--retrain_path", help="Retrain saved model",type=str)
+    parser.add_argument('-sparse','--sparse', help='Stores features as sparse matrices', default=False, action='store_true')
+    parser.add_argument('-val','--validation', help='Perform validation', default=True, action='store_false')
+    
+    args = parser.parse_args()
+
+    print('aasas',flush=True)
+    os.makedirs(args.model, exist_ok=True)
+    
+    in_dict={'cpu':args.cpu,'rate':args.rate, 'iters':args.iterations, 'size':args.size,\
+             'train_path':args.train, 'model':args.model, 'val':args.validation,'rt_path':args.retrain_path,\
+            'sparse':args.sparse}
+    
+    train_indel_model(in_dict)
+
+    elapsed=time.time()-t
+    print ('Total Time Elapsed: %.2f seconds' %elapsed)
\ No newline at end of file