bestref.py

#!/usr/bin/env python
# Copyright (C) <2020> PMBL;South China Agricultural University. All rights reserved
__author__ = "Write by Fangping Li"
__version__ = '0.1.0'

import argparse
import sys
import os
import collections

def get_options():    
    example = "bestref.py -i pangenome.fa -hap  goin.hapc  -g guide.goc -l location.lg -o output " 
    description = "Create a optimal reference genome based on short reads from pan-genome: " + example
    parser = argparse.ArgumentParser(description = description,prog = 'bestref.py')
    parser.add_argument('-i', '--inpan', action='store',type=str,help='input your reference .fasta')
    parser.add_argument('-o', '--output', action='store',type=str,
                         help='name of the  pan-genome coverage output',default="output")
    parser.add_argument('-g', '--gocguide', action='store',type=str,help='input your .goc file generated by mappingtools.py')
    parser.add_argument('-hap', '--hapcguide', action='store',type=str,help='input your .hapc file generated by mappingtools.py')
    return parser.parse_args()

def seqfileread(file):#readfa
    file = open(file,"r")
    lines = list(file.readlines())
    dic={}
    temp = ""
    for i in lines:
        if i.startswith(">"):
            temp = i.strip()[1:]
            dic[temp] = ""
            #print(temp)
        else:
            i = i.strip()
            dic[temp] = dic[temp] + i
    file.close()
    return dic

gosome = get_options()
seq = gosome.inpan
goc = gosome.gocguide
hap = gosome.hapcguide

seqdic = seqfileread(seq)

gocf = open(goc,"r")
gocl = list(gocf.readlines())
gocf.close()

dicgoc = {}
for i in gocl:
    if i.find("Start") == -1:
        i = i.split()
        dicgoc[i[0].strip()] = []
        if "more" in i:
            lc = i.index("more")
            locstart = i[2]
            locend = i[lc-2]
        else:
            locstart = i[2]
            locend = i[5]
            lc = 7    

        maplc = i[2:lc-1]
        #print(maplc)
        mapc = 2
        c = 0 
        for k in range(len(maplc)): #build the [[xxx,xxx],[xxx,xxx]] .goc dic
            if mapc%2 == 0:
                    #print(mapc)
                loc1 = maplc[mapc-2]
                    #print(loc1)
                loc2 = maplc[mapc-1]
                #print(loc1,loc2)
                
                dicgoc[i[0]].append([loc1,loc2])
                
            mapc += 1
                
    
                
print(dicgoc)
hapf = open(hap,"r")
hapl = list(hapf.readlines())
hapf.close()
cumlength = 0
temp = ""

for i in hapl[1:]:
    i = i.split()
    #dicgoc[i[0]] = []


    mappc = i[3:-1]
    #print(mappc)
    name = i[0].strip()
    seqname = i[0].split("-")[2]
    if seqname != temp:
        cumlength = 0
    for j in range(len(mappc)):
        if mappc[j] == "0":
            print(dicgoc[name])
            #print(name.strip())
            #print(dicgoc["1-1-chr10"])
            
            a = int(dicgoc[name][j][0])-cumlength
            b = int(dicgoc[name][j][1])-cumlength
            seqdic[seqname] = seqdic[seqname][:a] + seqdic[seqname][b:]
            cumlength += b - a

    temp = seqname
            
outputfile = open(gosome.output+".fasta","w")  

for i in seqdic.keys():
    print(">"+i,file = outputfile)
    print(seqdic[i], file = outputfile)    
outputfile.close()