1919from .parser import parse_ref , parse_bed , check_fasta
2020from .helper import logger , map_fusion_to_iso , fix_bed , generate_bed
2121from collections import defaultdict
22+ import tempfile
2223
2324__author__ = [
2425 'Xiao-Ou Zhang (zhangxiaoou@picb.ac.cn)' ,
3031
3132@logger
3233def annotate (options ):
34+ # create temporary annotated fusion file
35+ fusion_tmp = tempfile .TemporaryFile (mode = 'w+' )
3336 # annotate fusion junctions
34- annotate_fusion (options ['--ref' ], options ['--bed' ],
37+ annotate_fusion (options ['--ref' ], options ['--bed' ], fusion_tmp ,
3538 secondary_flag = options ['--low-confidence' ])
3639 # fix fusion juncrions
37- fix_fusion (options ['--ref' ], options ['--genome' ], options ['--output' ],
40+ fix_fusion (options ['--ref' ], options ['--genome' ], fusion_tmp , options ['--output' ],
3841 options ['--no-fix' ], secondary_flag = options ['--low-confidence' ])
3942
4043
41- def annotate_fusion (ref_f , junc_bed , secondary_flag = 0 , denovo_flag = 0 ):
44+ def annotate_fusion (ref_f , junc_bed , fusion_tmp , secondary_flag = 0 , denovo_flag = 0 ):
4245 """
4346 Align fusion juncrions to gene annotations
4447 """
@@ -48,96 +51,96 @@ def annotate_fusion(ref_f, junc_bed, secondary_flag=0, denovo_flag=0):
4851 fusion_bed = junc_bed
4952 fusions , fusion_index = parse_bed (fusion_bed ) # fusion junctions
5053 total = set ()
51- annotated_fusion_f = 'annotated_fusion.txt.tmp'
52- with open (annotated_fusion_f , 'w' ) as outf :
53- for chrom in chrom_info :
54- # overlap gene annotations with fusion juncrions
55- result = []
56- # overlap genes
57- if chrom in genes :
58- result += Interval .overlapwith (genes [chrom ].interval ,
59- fusions [chrom ])
60- # overlap novel genes in denovo mode
61- if denovo_flag and chrom in novel_genes :
62- result += Interval .overlapwith (novel_genes [chrom ].interval ,
63- fusions [chrom ])
64- for itl in result :
65- # extract gene annotations
66- iso = list ([x for x in itl [2 :] if x .startswith ('iso' )])
67- # for each overlapped fusion junction
68- for fus in itl [(2 + len (iso )):]:
69- reads = fus .split ()[1 ]
70- fus_start , fus_end = fusion_index [fus ]
71- fus_loc = '%s\t %d\t %d\t FUSIONJUNC/%s' % (chrom , fus_start ,
72- fus_end , reads )
73- edge_annotations = [] # first or last exon flag
74- secondary_exon = defaultdict (dict ) # secondary exons
75- annotate_flag = 0
76- for iso_id in iso :
77- g , i , c , s = iso_id .split ()[1 :]
78- start = gene_info [iso_id ][0 ][0 ]
79- end = gene_info [iso_id ][- 1 ][- 1 ]
80- # fusion junction excesses boundaries of gene
81- # annotation
82- if fus_start < start - 10 or fus_end > end + 10 :
83- if not secondary_flag :
84- continue
85- (fusion_info ,
86- index ,
87- edge ,
88- secondary ) = map_fusion_to_iso (fus_start ,
89- fus_end , s ,
90- gene_info [iso_id ])
91- if fusion_info :
92- annotate_flag += 1
93- bed_info = '\t ' .join ([fus_loc , '0' , s ,
94- str (fus_start ),
95- str (fus_start ), '0,0,0' ])
96- bed = '\t ' .join ([bed_info , fusion_info , g , i ,
97- index ])
98- if not edge : # not first or last exon
99- outf .write (bed + '\n ' )
100- total .add (fus )
101- else : # first or last exon
102- edge_annotations .append (bed )
103- elif secondary_flag and secondary is not None :
104- li , ri = secondary
105- gene = ':' .join ([g , s ])
106- if li is not None :
107- li = str (li )
108- secondary_exon ['left' ][gene ] = ':' .join ([i ,
109- li ])
110- if ri is not None :
111- ri = str (ri )
112- secondary_exon ['right' ][gene ] = ':' .join ([i ,
113- ri ])
114- if edge_annotations :
115- for bed in edge_annotations :
54+ outf = fusion_tmp
55+ for chrom in chrom_info :
56+ # overlap gene annotations with fusion juncrions
57+ result = []
58+ # overlap genes
59+ if chrom in genes :
60+ result += Interval .overlapwith (genes [chrom ].interval ,
61+ fusions [chrom ])
62+ # overlap novel genes in denovo mode
63+ if denovo_flag and chrom in novel_genes :
64+ result += Interval .overlapwith (novel_genes [chrom ].interval ,
65+ fusions [chrom ])
66+ for itl in result :
67+ # extract gene annotations
68+ iso = list ([x for x in itl [2 :] if x .startswith ('iso' )])
69+ # for each overlapped fusion junction
70+ for fus in itl [(2 + len (iso )):]:
71+ reads = fus .split ()[1 ]
72+ fus_start , fus_end = fusion_index [fus ]
73+ fus_loc = '%s\t %d\t %d\t FUSIONJUNC/%s' % (chrom , fus_start ,
74+ fus_end , reads )
75+ edge_annotations = [] # first or last exon flag
76+ secondary_exon = defaultdict (dict ) # secondary exons
77+ annotate_flag = 0
78+ for iso_id in iso :
79+ g , i , c , s = iso_id .split ()[1 :]
80+ start = gene_info [iso_id ][0 ][0 ]
81+ end = gene_info [iso_id ][- 1 ][- 1 ]
82+ # fusion junction excesses boundaries of gene
83+ # annotation
84+ if fus_start < start - 10 or fus_end > end + 10 :
85+ if not secondary_flag :
86+ continue
87+ (fusion_info ,
88+ index ,
89+ edge ,
90+ secondary ) = map_fusion_to_iso (fus_start ,
91+ fus_end , s ,
92+ gene_info [iso_id ])
93+ if fusion_info :
94+ annotate_flag += 1
95+ bed_info = '\t ' .join ([fus_loc , '0' , s ,
96+ str (fus_start ),
97+ str (fus_start ), '0,0,0' ])
98+ bed = '\t ' .join ([bed_info , fusion_info , g , i ,
99+ index ])
100+ if not edge : # not first or last exon
116101 outf .write (bed + '\n ' )
117- total .add (fus )
118- if secondary_flag and not annotate_flag :
119- for gene in secondary_exon ['left' ]:
120- if gene in secondary_exon ['right' ]:
121- left = secondary_exon ['left' ][gene ]
122- right = secondary_exon ['right' ][gene ]
123- g , s = gene .split (':' )
124- # for avoid dup, use fus_loc_new
125- fus_loc_new = fus_loc + '\t 0\t %s' % s
126- outf .write ('%s\t %s:%s\t %s:%s\n ' % (fus_loc_new ,
127- g , left , g ,
128- right ))
102+ total .add (fus )
103+ else : # first or last exon
104+ edge_annotations .append (bed )
105+ elif secondary_flag and secondary is not None :
106+ li , ri = secondary
107+ gene = ':' .join ([g , s ])
108+ if li is not None :
109+ li = str (li )
110+ secondary_exon ['left' ][gene ] = ':' .join ([i ,
111+ li ])
112+ if ri is not None :
113+ ri = str (ri )
114+ secondary_exon ['right' ][gene ] = ':' .join ([i ,
115+ ri ])
116+ if edge_annotations :
117+ for bed in edge_annotations :
118+ outf .write (bed + '\n ' )
119+ total .add (fus )
120+ if secondary_flag and not annotate_flag :
121+ for gene in secondary_exon ['left' ]:
122+ if gene in secondary_exon ['right' ]:
123+ left = secondary_exon ['left' ][gene ]
124+ right = secondary_exon ['right' ][gene ]
125+ g , s = gene .split (':' )
126+ # for avoid dup, use fus_loc_new
127+ fus_loc_new = fus_loc + '\t 0\t %s' % s
128+ outf .write ('%s\t %s:%s\t %s:%s\n ' % (fus_loc_new ,
129+ g , left , g ,
130+ right ))
129131 print ('Annotated %d fusion junctions!' % len (total ))
130132
131133
132- def fix_fusion (ref_f , genome_fa , out_file , no_fix , secondary_flag = 0 ,
134+ def fix_fusion (ref_f , genome_fa , fusion_tmp , out_file , no_fix , secondary_flag = 0 ,
133135 denovo_flag = 0 ):
134136 """
135137 Realign fusion juncrions
136138 """
137139 print ('Start to fix fusion junctions...' )
138140 fa = check_fasta (genome_fa )
139141 ref = parse_ref (ref_f , 2 )
140- annotated_fusion_f = 'annotated_fusion.txt.tmp'
142+ annotated_fusion_f = fusion_tmp
143+ annotated_fusion_f .seek (0 ) ## from the start
141144 fusions , fusion_names , fixed_flag = fix_bed (annotated_fusion_f , ref , fa ,
142145 no_fix , denovo_flag )
143146 total = 0
@@ -216,6 +219,6 @@ def fix_fusion(ref_f, genome_fa, out_file, no_fix, secondary_flag=0,
216219 if secondary_flag :
217220 secondary_f .close ()
218221
219- os . remove ( 'annotated_fusion.txt.tmp' )
222+ fusion_tmp . close ( )
220223
221224 print ('Fixed %d fusion junctions!' % total )
0 commit comments