QAseq-Fusion_v2.1/py_script/breakpoint_search_v6.py

#!/usr/bin/python
# -*- coding:utf-8 -*-

import argparse,time,os

BASE_DICT={'A':'T','G':'C','T':'A','C':'G','N':'N'}

def read_sam(insam):

    ID_dict={}
    with open(insam,'r') as inputSam:
        for line in inputSam:
            if not line.startswith("@"):
                info=line.strip("\n").split("\t")
                #ID_info = info[0].split("_")
                ID = info[0]
                if ID not in ID_dict.keys():
                    ID_dict[ID]=[info]
                else:
                    ID_dict[ID].append(info)
    return ID_dict

def check_clip(reads_info):

    if len(reads_info)>=3:
        return "hard_clip"
    elif len(reads_info)==2:
        for item in reads_info:
            if "S" in item[5]:
                return "soft_clip"
    return "no_clip"

def get_map_info(map_str):
    #99M51S
    num_list=[]
    map_dict={}
    count=0
    for item in map_str:
        if item.isalpha():
            num=len(num_list)
            index=0
            for i in range(num):
                index += int(num_list[i])*10**(num-i-1)
            map_dict[count]=item+":"+str(index)
            count +=1
            num_list=[]
        elif item.isdigit():
            num_list.append(item)
    return map_dict

def check_match(match1,match2):
    if match1<30 or match2<30:
        return False
    return True

def adj_info(seq,info1):

    flag = True
    for key in info1.keys():
        item = info1[key]
        if 'D' in item or 'I' in item:
            flag = False

    if 'H' in info1[0]:
        flag = False

    #不需要调整
    if flag:
        return seq,info1
    else:
        new_seq_list=[]
        new_info1={}
        start=0
        old_item=''
        for key in info1.keys():
            item = info1[key]
            temp_len=int(item.split(":")[1])
            end = start+temp_len
            if item.startswith('M') or item.startswith('S'):
                for index in range(start,end):
                    new_seq_list.append(seq[index])
                old_item=item

            elif item.startswith('D') and old_item.startswith('M'):
                end = end-temp_len
                for index in range(temp_len):
                    new_seq_list.append('-')

            elif item.startswith('H'):
                end = end-temp_len
            start = end

        #
        num=len(info1)
        index =0
        new_index= 0 
        while index+3 <= num:
            item1=info1[index]
            item2=info1[index+1]
            item3=info1[index+2]
            if item1.startswith("M") and item2.startswith("D") and item3.startswith("M"):
                pos1=int(item1.split(":")[1])
                pos2=int(item2.split(":")[1])
                pos3=int(item3.split(":")[1])
                new_info1[new_index] = 'M:'+str(pos1+pos2+pos3)
                new_index +=1
                index +=3

            elif item1.startswith("M") and item2.startswith("I") and item3.startswith("M"):
                pos1=int(item1.split(":")[1])
                pos3=int(item3.split(":")[1])
                new_info1[new_index] = 'M:'+str(pos1+pos3)
                new_index +=1
                index +=3

            elif item1.startswith("H") and item2.startswith("M"):
                new_info1[new_index] = item2
                new_index +=1
                index +=2
            else:
                new_info1[new_index] = item1
                new_index += 1
                index += 1
        
        while index < num:
            new_info1[new_index]=info1[index]
            new_index += 1
            index += 1
        return "".join(new_seq_list),new_info1

def get_points(info1,info2,OVERLAP=10):
    
    #OVERLAP默认值是10，默认不用这个参数值
    if len(info1)>=5 and len(info2)>=5:
        map_info1=get_map_info(info1[5])
        map_info2=get_map_info(info2[5])
        seq1=info1[9]
        seq2=info2[9]
        #print(map_info1)
        #print(seq1)
        info1_temp=adj_info(seq1,map_info1)
        seq=info1_temp[0]
        map_info1=info1_temp[1]
        #print(map_info1)
        #print(seq)


        #print(map_info2)
        info2_temp=adj_info(seq2,map_info2)
        seq2=info2_temp[0]
        map_info2=info2_temp[1]
        #print(map_info2)

        chr1=info1[2]
        pos1=int(info1[3])

        chr2=info2[2]
        pos2=int(info2[3])

        flag = int(info1[1])
        
        overlap=0
        breakpoint1=0
        breakpoint1_1=0
        breakpoint2=0

        breakpoint1_e=0
        breakpoint2_e=0


        fusion_seq=""

        if len(map_info1)==2 and len(map_info2)==2:
            #soft-clip end,hard-clip front
            #{0: 'M:131', 1: 'S:19'}
            #{0: 'H:118', 1: 'M:32'}
            if 'S' in map_info1[1] and 'H' in map_info2[0]:
                match1=int(map_info1[0].split(":")[1])
                skip1=int(map_info1[1].split(":")[1])
                match2=int(map_info2[1].split(":")[1])

                if not check_match(match1,match2):
                    return False,0
                
                overlap=abs(match2-skip1)

                breakpoint1=pos1+match1-1-overlap
                breakpoint2=pos2
                breakpoint1_1=pos1+match1-1

                if skip1 >30:
                    len2=30
                else:
                    len2=skip1
                fusion_seq=seq[match1-30:match1+len2]

                breakpoint1_e=breakpoint1_1-30
                breakpoint2_e=breakpoint2+30
                
            #soft-clip end,hard-clip end
            #{0: 'M:99', 1: 'S:51'}
            #{0: 'M:53', 1: 'H:97'}
            elif 'S' in map_info1[1] and 'H' in map_info2[1]:
                match1=int(map_info1[0].split(":")[1])
                skip1=int(map_info1[1].split(":")[1])
                match2=int(map_info2[0].split(":")[1])
                hide2=int(map_info2[1].split(":")[1])

                if not check_match(match1,match2):
                    return False,0
                
                overlap=abs(abs(match2-skip1)-abs(match1+skip1-match2-hide2))

                breakpoint1=pos1+match1-1-overlap
                breakpoint2=pos2+match2-1
                breakpoint1_1=pos1+match1-1

                if skip1 >30:
                    len2=30
                else:
                    len2=skip1
                #print(match1)
                fusion_seq=seq[match1-30:match1+len2]
                #print(fusion_seq)

                breakpoint1_e=breakpoint1_1-30
                breakpoint2_e=breakpoint2-30

            #soft-clip front,hard-clip end
            #{0: 'S:51', 1: 'M:99'}
            #{0: 'M:53', 1: 'H:97'}
            elif 'S' in map_info1[0] and 'H' in map_info2[1]:
                match1=int(map_info1[1].split(":")[1])
                skip1=int(map_info1[0].split(":")[1])
                match2=int(map_info2[0].split(":")[1])
                hide2=int(map_info2[1].split(":")[1])

                if not check_match(match1,match2):
                    return False,0
                
                overlap=abs(abs(match2-skip1)-abs(match1+skip1-match2-hide2))

                breakpoint1=pos1+overlap
                breakpoint2=pos2+match2-1
                breakpoint1_1=pos1

                fusion_seq=seq[skip1-30:skip1+30] if skip1 > 30 else seq[:skip1+30]
                
                breakpoint1_e=breakpoint1_1+30
                breakpoint2_e=breakpoint2-30

            #soft-clip front,hard-clip front
            #{0: 'S:51', 1: 'M:99'}
            #{0: 'H:97', 1: 'M:53'}
            elif 'S' in map_info1[0] and 'H' in map_info2[0]:

                match1=int(map_info1[1].split(":")[1])
                skip1=int(map_info1[0].split(":")[1])
                match2=int(map_info2[1].split(":")[1])

                if not check_match(match1,match2):
                    return False,0

                overlap=abs(match2-skip1)

                breakpoint1=pos1+overlap
                breakpoint2=pos2
                breakpoint1_1=pos1

                fusion_seq=seq[skip1-30:skip1+30] if skip1 > 30 else seq[:skip1+30]

                breakpoint1_e=breakpoint1_1+30
                breakpoint2_e=breakpoint2+30

        elif len(map_info1)==3 and len(map_info2)==2:

            flag_bin=bin(flag)
            if flag_bin[-5] == '1':
                #soft-clip front
                if 'H' in map_info2[0]:
                    #hard-clip front
                    #{0: 'S:51', 1: 'M:99', 2:'S:5'}
                    #{0: 'H:97', 1: 'M:53'}
                    match1=int(map_info1[1].split(":")[1])
                    skip1=int(map_info1[0].split(":")[1])
                    match2=int(map_info2[1].split(":")[1])

                    if not check_match(match1,match2):
                        return False,0

                    overlap=abs(match2-skip1)

                    breakpoint1=pos1+overlap
                    breakpoint2=pos2
                    breakpoint1_1=pos1

                    fusion_seq=seq[skip1-30:skip1+30] if skip1 > 30 else seq[:skip1+30]

                    breakpoint1_e=breakpoint1_1+30
                    breakpoint2_e=breakpoint2+30

                elif 'H' in map_info2[1]:
                    #hard-clip end
                    #{0: 'S:51', 1: 'M:99', 2:'S:5'}
                    #{0: 'M:53', 1: 'H:97'}
                    match1=int(map_info1[1].split(":")[1])
                    skip1=int(map_info1[0].split(":")[1])
                    match2=int(map_info2[0].split(":")[1])

                    if not check_match(match1,match2):
                        return False,0
                    
                    overlap=abs(match2-skip1)

                    breakpoint1=pos1+overlap
                    breakpoint2=pos2+match2-1
                    breakpoint1_1=pos1

                    fusion_seq=seq[skip1-30:skip1+30] if skip1 > 30 else seq[:skip1+30]

                    breakpoint1_e=breakpoint1_1+30
                    breakpoint2_e=breakpoint2-30

            elif flag_bin[-5] == '0':
                #soft-clip end
                if 'H' in map_info2[0]:
                    #hard-clip front
                    #{0: 'S:5', 1: 'M:99', 2:'S:51'}
                    #{0: 'H:97', 1: 'M:53'}
                    former_skip=int(map_info1[0].split(":")[1])
                    match1=int(map_info1[1].split(":")[1])
                    skip1=int(map_info1[2].split(":")[1])
                    match2=int(map_info2[1].split(":")[1])

                    if not check_match(match1,match2):
                        return False,0

                    overlap=abs(match2-skip1)

                    breakpoint1=pos1+match1-1-overlap
                    breakpoint2=pos2
                    breakpoint1_1=pos1+match1-1

                    if skip1 >30:
                        len2=30
                    else:
                        len2=skip1

                    fusion_seq=seq[former_skip+match1-30:former_skip+match1+len2]

                    breakpoint1_e=breakpoint1_1-30
                    breakpoint2_e=breakpoint2+30

                elif 'H' in map_info2[1]:
                    #hard-clip end
                    #{0: 'S:5', 1: 'M:99', 2:'S:51'}
                    #{0: 'M:53', 1: 'H:97'}
                    former_skip=int(map_info1[0].split(":")[1])
                    match1=int(map_info1[1].split(":")[1])
                    skip1=int(map_info1[2].split(":")[1])
                    match2=int(map_info2[0].split(":")[1])

                    if not check_match(match1,match2):
                        return False,0

                    overlap=abs(match2-skip1)

                    breakpoint1=pos1+match1-1-overlap
                    breakpoint2=pos2+match2-1
                    breakpoint1_1=pos1+match1-1

                    if skip1 >30:
                        len2=30
                    else:
                        len2=skip1
                    fusion_seq=seq[former_skip+match1-30:+former_skip+match1+len2]

                    breakpoint1_e=breakpoint1_1-30
                    breakpoint2_e=breakpoint2-30

        if overlap <= OVERLAP and breakpoint1 !=0:
            return True,chr1+":"+str(breakpoint1),chr2+":"+str(breakpoint2),chr1+":"+str(breakpoint1_1),fusion_seq,str(overlap),chr1+":"+str(breakpoint1_e),chr2+":"+str(breakpoint2_e)

        return False,0
    return False,0

def hard_clip_fusion(reads_info):

    #R1:hard-clip
    R1_soft_info=reads_info[0]
    R1_hard_info=reads_info[1]
    if int(R1_hard_info[1])>256 and int(R1_soft_info[4])>=30:
        results=get_points(R1_soft_info,R1_hard_info)
        if results[0]:
            out_list=[results[1]+"-"+results[2],results[3],results[2],results[4],results[5],results[6],results[7]]
            return out_list
    #R2:hard-clip
    R2_soft_info=reads_info[-2]
    R2_hard_info=reads_info[-1]
    if int(R2_hard_info[1])>256 and int(R2_soft_info[4])>=30:
        results=get_points(R2_soft_info,R2_hard_info)
        if results[0]:
            out_list=[results[1]+"-"+results[2],results[3],results[2],results[4],results[5],results[6],results[7]]
            return out_list
    return []

def primer_loc(primer,seq):
    primer_num=len(primer)
    seq_primer=seq[:primer_num]
    seq_rev_comp="".join([BASE_DICT[base] for base in seq][::-1])
    seq_rev_comp_primer=seq_rev_comp[:primer_num]
    if primer == seq_primer:
        return "Front"
    elif primer == seq_rev_comp_primer:
        return "End"
    return "NO"

def get_pos_end(map_info,pos):
    
    chr=pos.split(":")[0]
    pos_temp=int(pos.split(":")[1])
    new_pos=0
    correct_pos=pos_temp
    if map_info[0][0]=='M':
        match_len=int(map_info[0].split(":")[1])
        pos_temp=pos_temp+match_len-1
        correct_pos=pos_temp
        new_pos=pos_temp-30
    elif map_info[0][0]=='S':
        new_pos=pos_temp+30
    else:
        return pos

    return chr+":"+str(correct_pos),chr+":"+str(new_pos)

def soft_clip_fusion(reads_info):

    R1_read_info=reads_info[0]
    R2_read_info=reads_info[1]
    ID_info = R1_read_info[0].split("_")
    ID = R1_read_info[0]
    R1_primer_loc=""
    R1_seq=R1_read_info[9]
    R1_mismatch_seq=""

    #比对质量值至少一个MAQ>=60
    if int(R1_read_info[4])<60 and int(R2_read_info[4])<60:
        return False,"",""

    R1_loc=""
    R1_map_info={}
    R2_map_info={}
    #R1 primer存在
    if len(ID_info)>=4:
        R1_primer=ID_info[3]
        R1_map_info=get_map_info(R1_read_info[5])
        R1_loc=primer_loc(R1_primer,R1_read_info[9])
        if len(R1_map_info)==3 and R1_map_info[0][0] == "S" and R1_map_info[2][0]=="S":
            if R1_loc == "End":
                R1_primer_loc="Not the Same"
                R1_seq="".join([BASE_DICT[base] for base in R1_seq][::-1])
                mismatch_len=int(R1_map_info[0].split(":")[1])
                R1_mismatch_seq=R1_seq[-1*mismatch_len:]

            if R1_loc == "Front":
                R1_primer_loc="Not the Same"
                mismatch_len=int(R1_map_info[2].split(":")[1])
                R1_mismatch_seq=R1_seq[-1*mismatch_len:]

        elif len(R1_map_info)==2:
            if R1_loc == "End" and R1_map_info[0][0]=="S":
                R1_primer_loc="Not the Same"
                R1_seq="".join([BASE_DICT[base] for base in R1_seq][::-1])
                mismatch_len=int(R1_map_info[0].split(":")[1])
                R1_mismatch_seq=R1_seq[-1*mismatch_len:]

            elif R1_loc == "End" and R1_map_info[1][0]=="S":
                R1_primer_loc="Same"
                R1_seq="".join([BASE_DICT[base] for base in R1_seq][::-1])
                mismatch_len=int(R1_map_info[1].split(":")[1])

                if len(R1_primer)+5 >=mismatch_len:
                    return False,"",""
                R1_mismatch_seq=R1_seq[:mismatch_len]

            elif R1_loc == "Front" and R1_map_info[0][0]=="S":
                R1_primer_loc="Same"
                mismatch_len=int(R1_map_info[0].split(":")[1])
                if len(R1_primer)+5 >=mismatch_len:
                    return False,"",""
                R1_mismatch_seq=R1_seq[:mismatch_len]

            elif R1_loc == "Front" and R1_map_info[1][0] =="S":
                R1_primer_loc="Not the Same"
                mismatch_len=int(R1_map_info[1].split(":")[1])
                R1_mismatch_seq=R1_seq[-1*mismatch_len:]

    R2_primer_loc=""
    R2_seq=R2_read_info[9]
    R2_mismatch_seq=""
    R2_loc=""
    R2_map_info={}
    #R1,R2都存在primer
    if len(ID_info)==5:
        R2_primer=ID_info[4]
        R2_map_info=get_map_info(R2_read_info[5])
        R2_loc=primer_loc(R2_primer,R2_read_info[9])
        if len(R2_map_info)==3 and R2_map_info[0][0] == "S" and R2_map_info[2][0]=="S":
            if R2_loc == "End":
                R2_primer_loc="Not the Same"
                mismatch_len=int(R2_map_info[0].split(":")[1])
                R2_mismatch_seq=R2_seq[:mismatch_len]

            if R2_loc == "Front":
                R2_primer_loc="Not the Same"
                R2_seq="".join([BASE_DICT[base] for base in R2_seq][::-1])
                mismatch_len=int(R2_map_info[2].split(":")[1])
                R2_mismatch_seq=R2_seq[:mismatch_len]

        elif len(R2_map_info)==2:
            if R2_loc == "End" and R2_map_info[0][0]=="S":
                R2_primer_loc="Not the Same"
                mismatch_len=int(R2_map_info[0].split(":")[1])
                R2_mismatch_seq=R2_seq[:mismatch_len]

            elif R2_loc == "End" and R2_map_info[1][0]=="S":
                R2_primer_loc="Same"
                mismatch_len=int(R2_map_info[1].split(":")[1])
                if len(R2_primer)+5 >=mismatch_len:
                    return False,"",""
                R2_mismatch_seq=R2_seq[-1*mismatch_len:]

            elif R2_loc == "Front" and R2_map_info[0][0]=="S":
                R2_primer_loc="Same"
                R2_seq="".join([BASE_DICT[base] for base in R2_seq][::-1])
                mismatch_len=int(R2_map_info[0].split(":")[1])
                if len(R2_primer)+5 >= mismatch_len:
                    return False,"",""
                R2_mismatch_seq=R2_seq[-1*mismatch_len:]

            elif R2_loc == "Front" and R2_map_info[1][0] =="S":
                R2_primer_loc="Not the Same"
                R2_seq="".join([BASE_DICT[base] for base in R2_seq][::-1])
                mismatch_len=int(R2_map_info[1].split(":")[1])
                R2_mismatch_seq=R2_seq[:mismatch_len]


    if R1_primer_loc=="Not the Same" and R2_primer_loc == "Not the Same":
        R1_mismatch_len=len(R1_mismatch_seq)
        R2_mismatch_len=len(R2_mismatch_seq)
        if R1_mismatch_len>=5 and R2_mismatch_len >=5:
            for index in range(5):
                R2_match_seq = R2_seq[index+R2_mismatch_len:index+R2_mismatch_len+R1_mismatch_len]
                if R2_match_seq == R1_mismatch_seq:
                    R2_match_seq_tmp = R2_seq[:index+R2_mismatch_len+R1_mismatch_len]
                    R1_match_seq_tmp = R1_seq[-1*len(R2_match_seq_tmp):]
                    if R2_match_seq_tmp == R1_match_seq_tmp:
                        FusionSeq=R1_seq[-1*R1_mismatch_len-30:-1*R1_mismatch_len]+R2_seq[R2_mismatch_len+index:R2_mismatch_len+index+30]
                        R1_point=R1_read_info[2]+":"+R1_read_info[3]
                        if R2_loc == "End":
                            R2_point=R2_read_info[2]+":"+str(int(R2_read_info[3])+index)
                        elif R2_loc == "Front":
                            R2_point=R2_read_info[2]+":"+str(int(R2_read_info[3])-index)
                        else:
                            return True,"Double",outstr
                        R1_point_list=get_pos_end(R1_map_info,R1_point)
                        R2_point_list=get_pos_end(R2_map_info,R2_point)
                        outstr=R1_point+"-"+R2_point+"\t"+R1_read_info[2]+":"+R1_read_info[3]+"\t"+R2_read_info[2]+":"+R2_read_info[3]+"\t"+FusionSeq+"\t"+str(index)+"\t"+R1_point_list[1]+"\t"+R2_point_list[1]+"\t"+ID+"\n"
                        return True,"Double",outstr

    elif R1_primer_loc=="Not the Same" and R2_primer_loc == "Same":
        R1_mismatch_len=len(R1_mismatch_seq)
        R2_mismatch_len=len(R2_mismatch_seq)
        common_len=min(R1_mismatch_len,R2_mismatch_len)
        R1_mismatch_tmp=R1_mismatch_seq[-1*R1_mismatch_len:-1*common_len]
        R2_mismatch_tmp=R2_mismatch_seq[-1*R2_mismatch_len:-1*common_len]
        R2_match_seq=R2_seq[:-1*R2_mismatch_len]
        if R1_mismatch_tmp==R2_mismatch_tmp:
            R2_point=R2_read_info[2]+":"+R2_read_info[3]
            R2_point_list=get_pos_end(R2_map_info,R2_point)
            otustr=R2_point_list[0]+"\tDown\t"+R2_mismatch_seq+"\t"+R2_match_seq+"\t"+R2_point_list[1]+"\t"+ID+"\n"
            return True,"Single",otustr

    elif R1_primer_loc=="Same" and R2_primer_loc == "Same":
        return False,"",""
    
    elif R1_primer_loc=="Same" and R2_primer_loc == "Not the Same":
        R1_mismatch_len=len(R1_mismatch_seq)
        R2_mismatch_len=len(R2_mismatch_seq)
        common_len=min(R1_mismatch_len,R2_mismatch_len)
        R1_mismatch_tmp=R1_mismatch_seq[R1_mismatch_len-1*common_len:R1_mismatch_len]
        R2_mismatch_tmp=R2_mismatch_seq[R2_mismatch_len-1*common_len:R2_mismatch_len]
        R1_match_seq=R1_seq[R1_mismatch_len:]
        if R1_mismatch_tmp == R2_mismatch_tmp:
            R1_point=R1_read_info[2]+":"+R1_read_info[3]
            R1_point_list=get_pos_end(R1_map_info,R1_point)
            otustr=R1_point_list[0]+"\tUp\t"+R1_mismatch_seq+"\t"+R1_match_seq+"\t"+R1_point_list[1]+"\t"+ID+"\n"
            return True,"Single",otustr
    else:
        return False,"",""

    return False,"",""

def main(insam,outdir):

    t1=time.time()
    #将sam内容按照ID区分，放到ID_info中
    ID_dict=read_sam(insam)
    #存放soft_clip的reads ID
    soft_clip_ID=[]

    #输出到本地hard clip fusion
    hard_clip_file=outdir+"/hard_clip_fusion_breakpoints.txt"
    hard_clip_out=open(hard_clip_file,'w')
    for ID in ID_dict.keys():
        reads_info = ID_dict[ID]
        clip_stat=check_clip(reads_info)
        if clip_stat == "soft_clip":
            soft_clip_ID.append(ID)
        elif clip_stat == "hard_clip":
            out_list=hard_clip_fusion(reads_info)
            if len(out_list) !=0:
                out_list.append(ID)
                hard_clip_out.write("\t".join(out_list)+"\n")
    hard_clip_out.close()

    #输出到本地single
    single_points=outdir+"/soft_clip_single_breakpoints.txt"
    SingleInput=open(single_points,'w')
    #输出到本地double
    double_points=outdir+"/soft_clip_double_breakpoints.txt"
    DoubleInput=open(double_points,'w')
    
    #
    for soft_ID in soft_clip_ID:
        reads_info = ID_dict[soft_ID]
        soft_clip_results=soft_clip_fusion(reads_info)
        if soft_clip_results[0]:
            if soft_clip_results[1]=="Double":
                DoubleInput.write(soft_clip_results[2])
            elif soft_clip_results[1]=="Single":
                SingleInput.write(soft_clip_results[2])
    SingleInput.close()
    DoubleInput.close()

    double_points_file=outdir+"/double_breakpoints.txt"
    os.system("cat "+hard_clip_file+" > "+double_points_file+" && cat "+double_points+" >> "+double_points_file)

    single_points_file=outdir+"/single_breakpoints.txt"
    os.system("cp "+single_points+" "+single_points_file)

    t2=time.time()
    print("Times: "+str(int(t2-t1))+"s")
    print("BreakPoints Search Done!")


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='get reads mapping location')
    parser.add_argument('-i', required=True, type=str, help="insam")
    parser.add_argument('-o', required=True, type=str, help="outdir")
    args = parser.parse_args()
    main(args.i, args.o)