comparison pyPRADA_1.2/prada-fusion @ 0:acc2ca1a3ba4

Uploaded

0:acc2ca1a3ba4

#!/usr/bin/env python

#PRADA: Pipeline for RnAseq Data Analysis
#Fusion detection module, algorithm published by Michael Berger et al (Genome Res, 2010) at Broad Institute. 
#Implemented by Siyuan Zheng, Wandaliz Torres-Garcia and Rahul Vegesna. 
#Copy Right: The Univ of Texas MD Anderson Cancer Center, Department of Bioinformatics and Computational Biology
#Contact: Roel Verhaak ([e-mail])
#Citation: to be added
#Tested with Python v2.6 and v2.7 
#The program requires NM tag and high quality mapping reads with mapping score more than -minmapq. 
#Last modified: 04/11/2013

######################################################################################
import sys
import time
import os
import os.path
import subprocess
import operator
import pysam
import bioclass
import gfclass 
import ioprada
import privutils
from Bio import SeqIO,Seq

######################################################################################
args=sys.argv

help_menu='''\nPipeline for RNAseq Data Analaysis - fusion detection (PRADA).
    **Command**:
    prada-fusion -bam XX.bam -conf xx.txt -tag XX -mm 1 -junL XX -minmapq 30 -outdir ./ 
    **Parameters**:
    -h		print help message
    -bam	input BAM file, must has a .bam suffix. BAM is the output from PRADA preprocess module.
    -conf	config file for references and parameters. Use conf.txt in py-PRADA installation folder if none specified.
    -tag	a tag to describe the sample, used to name output files. Default is ''.
    -mm		number of mismatches allowed in discordant pairs and fusion spanning reads.Default is 1.  
    -junL	length of sequences taken from EACH side of exons when making hypothetical junctions. No default.
    -minmapq	minimum read mapping quality to be considered as fusion evidences. Default is 30.
    -outdir	output directory.
    -v		print version
'''

if '-h' in args or '-help' in args or len(args)==1:
    print help_menu
    sys.exit(0)

if '-v' in args:
    import version
    print version.version
    sys.exit(0)

if '-bam' not in args:
    sys.exit('ERROR: BAM file needed')
i=args.index('-bam')
bampath=args[i+1]
bampath=os.path.abspath(bampath)
bam=os.path.basename(bampath)
if bam[-4:] != '.bam':
    sys.exit('ERROR: BAM file must have suffix .bam')

#Mismatch allowed. This filter is applied at the very end of the pipeline.
#I strongly suggest 1. We also record how many junction spanning reads (JSRs) are perfect matched.
if '-mm' not in args:
    mm=1
else:
    i=args.index('-mm')
    mm=int(args[i+1])

#junL should be less than the read length in the dataset. I suggest it to be read_length*0.8
if '-junL' not in args:
    sys.exit('ERROR: junL must be specified')
i=args.index('-junL')
overlap=int(args[i+1])

#minimum mapping quality for reads as fusion evidences
if '-minmapq' not in args:
    minmapq=30
else:
    i=args.index('-minmapq')
    minmapq=int(args[i+1])

if '-outdir' not in args:
    outpath=os.path.abspath('./')
else:
    i=args.index('-outdir')
    outpath=os.path.abspath(args[i+1])
if os.path.exists(outpath):
    print 'WARNING: outdir %s exists'%outpath
else:
    os.mkdir(outpath)

if '-tag' not in args:
    docstring='prada'
else:
    i=args.index('-tag')
    docstring=args[i+1]

#by default, search conf.txt in the prada folder
prada_path=os.path.dirname(os.path.abspath(__file__))   #path to find libraries and executives.
ref_search_path=[prada_path,os.getcwd()]                #search path for conf file if not specified in command line

if '-conf' in args:
    i=args.index('-ref')
    reffile=args[i+1]
    if os.path.exists(reffile):
        pass
    else:
        for pth in ref_search_path:
            new_reffile='%s/%s'%(pth, os.path.basename(reffile))  
            if os.path.exists(new_reffile):
                reffile=new_reffile
                break
        else:
            sys.exit('ERROR: conf file %s not found'%reffile)
else:
    reffile='%s/conf.txt'%prada_path
    if not os.path.exists(reffile):
        sys.exit('ERROR: No default conf.txt found and none specified')

#Now print all input parameters
print 'CMD: %s'%('\t'.join(args))

#reference files pointing to the annotation files. 
refdict=ioprada.read_conf(reffile)
featurefile=refdict['--REF--']['feature_file']
txseqfile=refdict['--REF--']['tx_seq_file']
txcatfile=refdict['--REF--']['txcat_file']
cdsfile=refdict['--REF--']['cds_file']

#underlying utilities, automatically detected
#these are customized tools. update is needed. 
samtools='%s/tools/samtools-0.1.16/samtools'%prada_path
bwa='%s/tools/bwa-0.5.7-mh/bwa'%prada_path
blastn='%s/tools/blastn'%prada_path

######################################################################################
print 'step 0: loading gene annotations @ %s'%time.ctime()
#call functions in ioprada module // 
txdb,genedb=ioprada.read_feature(featurefile,verbose=True)
tx_primary=ioprada.read_tx_cat(txcatfile)
tx_cds=ioprada.read_cds(cdsfile)

######################################################################################
print 'step 1: finding discordant pairs @ %s'%time.ctime()

#We sift through all exons of protein coding genes and get the mapping reads. 
#Within, we exclude low mapping quality reads and PCR duplicates. For pairs that the two ends
#map to different genes, we all it a discordant pair. 
#This is a step for finding all possible candidate fusions. 

samfile=pysam.Samfile(bampath,'rb')

read1_ab={}
read2_ab={}
db1={}
db2={}

i,N=0,len(genedb)
for gene in genedb:
    i+=1
    if i%200==0:
        print '%d/%d genes processed for discordant pairs'%(i,N)
    g=genedb[gene]
    exons=g.get_exons()
    for e in exons.values():
        for rd in samfile.fetch(e.chr,e.start-1,e.end):
            if rd.mapq < minmapq:
                continue
            if rd.is_duplicate:
                continue
            if rd.mate_is_unmapped:   #at this point, only consider pairs
                continue
            if rd.rnext == rd.tid and rd.mpos <= g.end and rd.mpos >= g.start-1: #remove reads that fall into the same gene range
                continue
            if rd.is_read1:
                if read1_ab.has_key(rd.qname):
                    read1_ab[rd.qname].add(gene)
                else:
                    read1_ab[rd.qname]=set([gene])
                db1[rd.qname]=rd
            if rd.is_read2:
                if read2_ab.has_key(rd.qname):
                    read2_ab[rd.qname].add(gene)
                else:
                    read2_ab[rd.qname]=set([gene])
                db2[rd.qname]=rd

##output the discordant pairs and determine the orientation of candidate fusions
discordant={}   #catalogue all discordant pairs, using gene pairs as keys
outfile=open('%s/%s.discordant.txt'%(outpath,docstring),'w')
title=['read','gene1','gene1_chr','read1_pos','read1_mm','read1_strand','read1_orient','gene2',\
      'gene2_chr','read2_pos','read2_mm','read2_strand','read2_orient']
outfile.write('%s\n'%('\t'.join(title)))
i=0
for rdid in read1_ab:
    i+=1
    if i%10000==0:
        print '%d discordant pairs processed'%i
    if not read2_ab.has_key(rdid):   #skip if not all ends are catalogued
        continue
    g1set=read1_ab[rdid]    #consider all combinations if a read maps to multiple genes
    g2set=read2_ab[rdid]    #consider all combinations if a read maps to multiple genes
    r1,r2=db1[rdid],db2[rdid]
    read1strd='-1' if r1.is_reverse else '1'
    read2strd='-1' if r2.is_reverse else '1'
    for g1 in g1set:
        for g2 in g2set:
            if g1==g2:    #for some uncasual cases
                continue
            g1obj,g2obj=genedb[g1],genedb[g2]
            read1orient='F' if read1strd == g1obj.strand else 'R' #read1 --> gene1
            read2orient='F' if read2strd == g2obj.strand else 'R' #read2 --> gene2
            fkey=''
            if read1orient=='F' and read2orient=='R':    ##scenario I, gene1-gene2
                fkey=g1+'_'+g2
            if read1orient=='R' and read2orient=='F':    ##scenario II, gene2-gene1
                fkey=g2+'_'+g1
            if fkey:
                if discordant.has_key(fkey):
                    discordant[fkey].update({rdid:'%s:%s'%(read1orient,read2orient)})
                else:
                    discordant[fkey]={rdid:'%s:%s'%(read1orient,read2orient)}
                ##output
                nm1=str([x[1] for x in r1.tags if x[0]=='NM'][0])   #output mismatch, but does not consider it at this point
                nm2=str([x[1] for x in r2.tags if x[0]=='NM'][0])
                uvec=[rdid,g1,g1obj.chr,str(r1.pos+1),nm1,read1strd,read1orient,g2,g2obj.chr,str(r2.pos+1),nm2,read2strd,read2orient]
                outfile.write('%s\n'%('\t'.join(uvec)))
outfile.close()
      
##########################################################################
print 'step 2: finding recurrent pairs (candidates) @ %s'%time.ctime()    

#step 2 finds all candidates that have at least 2 discordant pairs. Meanwhile, filter out potential
#read through events. read through is defined as reads with mapping position less than 1M, while meeting
#the strand expectation.
 
guess=[]
outfile=open('%s/%s.recurrent.txt'%(outpath,docstring),'w')
title=['geneA','geneA_chr','geneB','geneB_chr','num_pairs','IDs']
outfile.write('%s\n'%('\t'.join(title)))
for pp in discordant:
    if len(discordant[pp]) < 2:   #consider only "recurrent" (more than 1 pair support) cases
        continue
    gene1,gene2=pp.split('_')
    g1obj,g2obj=genedb[gene1],genedb[gene2]
    rdset=discordant[pp].keys()
    #filter read-through
    #readthrough is defined at read level, regardless of mapping genes
    for rd in rdset:
        r1,r2=db1[rd],db2[rd]
        read1strd='-1' if r1.is_reverse else '1'
        read2strd='-1' if r2.is_reverse else '1'
        readthrough=False
        if db1[rd].tid == db2[rd].tid and abs(db1[rd].pos - db2[rd].pos) <= 1000000:
            if discordant[pp][rd]=='F:R':
                if read1strd=='1' and read2strd=='-1' and db1[rd].pos < db2[rd].pos:
                    readthrough=True
                if read1strd=='-1' and read2strd=='1' and db1[rd].pos > db2[rd].pos:
                    readthrough=True
            if discordant[pp][rd]=='R:F':
                if read2strd=='1' and read1strd=='-1' and db2[rd].pos < db1[rd].pos:
                    readthrough=True
                if read2strd=='-1' and read1strd=='1' and db2[rd].pos > db1[rd].pos:
                    readthrough=True
            if readthrough:
                del discordant[pp][rd]    #in-place deletion!!!! Change the discordant variable in place!
    if len(discordant[pp]) < 2:    #skip all that have less than 2 supporting discordant read pairs after readthrough filter
        continue
    guess.append(pp)
    #output
    uvec2=[gene1,g1obj.chr,gene2,g2obj.chr,str(len(discordant[pp])),'|'.join(discordant[pp])]
    outfile.write('%s\n'%('\t'.join(uvec2)))
outfile.close()

##########################################################################
print 'step 3: finding potential junction spanning reads @ %s'%time.ctime()

#For all candidates, find potential junction spanning reads (JSRs). A JSR is defined as a unmapped read but with the mate mapping
#to either F or R partner, with high mapping quality. Since the JSR is unmapped, it is not practical to consider PCR duplicate 
#because they are not properly marked. 

Fpartners=set()
Rpartners=set()
for pp in guess:
    gs=pp.split('_')
    Fpartners.add(gs[0])
    Rpartners.add(gs[1])
AllPartners=Fpartners|Rpartners

samfile.reset()
posjun={}    ##catalogue all JSRs, with track of the mate mapping genes and orientation.
i,N=0,len(AllPartners)
for gene in AllPartners:
    i+=1
    if i%200==0:
        print '%d/%d genes processed for potential junc reads'%(i,N)
    g=genedb[gene]
    exons=g.get_exons().values()
    for e in exons:
        for rd in samfile.fetch(e.chr,e.start-1,e.end):
            if rd.mapq < minmapq:
                continue
            if not rd.mate_is_unmapped:
                continue
            readstrd='-1' if rd.is_reverse else '1'
            orient='F' if readstrd == g.strand else 'R'  #mapping info of mate read. JSR per se is unmapped in BAM
            posjun[rd.qname]={'gene':gene,'orient':orient}

samfile.reset()
outfile=open('%s/%s.pos_junc_unmapped.fastq'%(outpath,docstring),'w')
i,N=0,len(posjun)
for rd in samfile:
    if rd.mate_is_unmapped:  #since the read is potential jun spanning read, all mate map to A or B
        continue
    if rd.is_unmapped:
        if posjun.has_key(rd.qname):
            i+=1
            if i%10000==0:
                print 'extracted %d/%d potential junc reads'%(i,N)
            rdname='%s_prada_%s_prada_%s'%(rd.qname,posjun[rd.qname]['gene'],posjun[rd.qname]['orient']) #_prada_ as split
            outfile.write('@%s\n'%rdname)
            outfile.write('%s\n'%rd.seq)
            outfile.write('+\n')
            outfile.write('%s\n'%rd.qual)
outfile.close()

######################################################################################
print 'step 4: building junction database @ %s'%time.ctime()

#Make hypothetical junctions between candidate fusion partners. To improve speed, we make a big junction database comprising 
#exons of all candidates, instead of by candidate individually. This also gives the possibility to assess the JSR mapping ambiguity
#across many junctions. It turned out very useful in filtering out false positives. 
#Note that in PRADA transcript sequence file, all sequences are + strand sequences. For - strand transcript, one need to 
#reverse complement the sequence to get the real transcript sequences. 

seqdb={}
for record in SeqIO.parse(txseqfile,'fasta'):
    seqdb[record.name]=record

outfile=open('%s/%s.junction.fasta'%(outpath,docstring),'w')
i,N=0,len(guess)
for pp in guess:
    i+=1
    if i%100==0:
        print 'building junction for %d/%d pairs'%(i,N)
    gene1,gene2=pp.split('_')
    g1obj,g2obj=genedb[gene1],genedb[gene2]
    eset1=g1obj.get_exons()    #unique exons in gene 1
    eset2=g2obj.get_exons()    #unique exons in gene 2
    #collect unique junctions
    juncseqdict={}    #save junction sequences
    for e1 in eset1.values():
        for e2 in eset2.values():
            e1_jun_name='%s:%s:%s'%(gene1,e1.chr,e1.end) if e1.strand=='1' else '%s:%s:%s'%(gene1,e1.chr,e1.start)
            e2_jun_name='%s:%s:%s'%(gene2,e2.chr,e2.start) if e2.strand=='1' else '%s:%s:%s'%(gene2,e2.chr,e2.end)
            jun_name=e1_jun_name+'_'+e2_jun_name
            tx1,tx2=txdb[e1.transcript],txdb[e2.transcript]
            max_a=tx1.exon_relative_pos()[e1.name][1]
            min_a=0 if max_a - overlap < 0 else max_a - overlap
            min_b=tx2.exon_relative_pos()[e2.name][0]-1
            max_b=tx2.length if min_b + overlap > tx2.length else min_b + overlap
            #reverse complementary when necessary
            try:
                tx1seq=seqdb[tx1.name].seq.tostring() if tx1.strand=='1' else seqdb[tx1.name].reverse_complement().seq.tostring()
                tx2seq=seqdb[tx2.name].seq.tostring() if tx2.strand=='1' else seqdb[tx2.name].reverse_complement().seq.tostring()
            except KeyError:   #in case transcript not found in sequence file
                continue
            jun_seq=tx1seq[min_a:max_a]+tx2seq[min_b:max_b]
            juncseqdict[jun_name]=jun_seq
    for junc in juncseqdict:
        outfile.write('>%s\n'%junc)
        outfile.write('%s\n'%juncseqdict[junc])
outfile.close()
samfile.close()

#for memory efficiecy, del seqdb
del seqdb

########################################################################################
print 'step 5: aligning potential junction reads to junction database @ %s'%time.ctime()

#Mapping potential JSRs to hypothetical junction database. 
#Allow 4 mismatches at the beginning.

idx_cmd='%s index %s/%s.junction.fasta'%(bwa,outpath,docstring)
os.system(idx_cmd)
aln_cmd='%s aln -n 4 -R 100 %s/%s.junction.fasta %s/%s.pos_junc_unmapped.fastq > %s/%s.juncmap.sai'%(bwa,outpath,docstring,outpath,docstring,outpath,docstring)
os.system(aln_cmd)
samse_cmd='%s samse -n 1000 -s %s/%s.junction.fasta %s/%s.juncmap.sai %s/%s.pos_junc_unmapped.fastq > %s/%s.juncmap.sam'%(bwa,outpath,docstring,outpath,docstring,outpath,docstring,outpath,docstring)
os.system(samse_cmd)
sam2bam_cmd='%s view -bS %s/%s.juncmap.sam -o %s/%s.juncmap.bam'%(samtools,outpath,docstring,outpath,docstring)
os.system(sam2bam_cmd)

jsam=pysam.Samfile('%s/%s.juncmap.bam'%(outpath,docstring),'rb')
#get the junction name directory
junctions=jsam.references
junname=dict(zip(range(0,len(junctions)),junctions))   #this is essential for quick speed.
junc_align={}

#go through the BAM file for meaningful (meeting fusion orientation etc) reads
strd_right_reads={}
rdb={}    #collect all junction spanning reads
i=0
for rd in jsam:
    i+=1
    if i%100000==0:
        print '%d junction alignments parsed'%i
    if rd.is_unmapped:
        continue
    read,mate_gene,mate_orient=rd.qname.split('_prada_')
    junc=junname[rd.tid]
    tmp=junc.split('_')
    gene1,gene2=[x.split(':')[0] for x in tmp]
    if gene1==mate_gene:
        if mate_orient=='F':
            if rd.is_reverse:
                if strd_right_reads.has_key(rd.qname):
                    strd_right_reads[rd.qname]+=1    #count how many times the read maps
                else:  
                    strd_right_reads[rd.qname]=1
                rdb[rd.qname]={'read':rd,'gene1':gene1,'gene2':gene2,'junc':junc} #will overwrite, but it is OK since we only look at unique ones
    elif gene2==mate_gene:
        if mate_orient == 'R':
            if not rd.is_reverse:
                if strd_right_reads.has_key(rd.qname):
                    strd_right_reads[rd.qname]+=1
                else:
                    strd_right_reads[rd.qname]=1
                rdb[rd.qname]={'read':rd,'gene1':gene1,'gene2':gene2,'junc':junc} #will overwrite, but it is OK since we only look at unique ones

#find uniquely mapped reads and their gene pairs
junc_map={}  #a dictionary from junction to mapping reads
for rdname in strd_right_reads:
    if strd_right_reads[rdname] > 1: #remove non-unique junction spanning reads
        continue
    infodict=rdb[rdname]
    pp=infodict['gene1']+'_'+infodict['gene2']
    if junc_map.has_key(pp):
        junc_map[pp].add(rdname)
    else:
        junc_map[pp]=set([rdname])
    
########################################################################################
print 'step 6: summarizing fusion evidences @ %s'%time.ctime()

#Now, time to apply mismatch filter and summarize the results
#Candidate fusions --> guess
#Discordant pairs --> discordant, db1, db2
#Junction reads --> junc_map, rdb
#Gene info --> genedb

outfile_s=open('%s/%s.fus.candidates.txt'%(outpath,docstring),'w')
outfile_d=open('%s/%s.fus.evidences.txt'%(outpath,docstring),'w')

title=['Gene_A','Gene_B','A_chr','B_chr','A_strand','B_strand','Discordant_n','JSR_n','perfectJSR_n','Junc_n','Position_Consist','Junction']
outfile_s.write('%s\n'%('\t'.join(title)))

for pp in junc_map: #all pairs with junc reads
    gene1,gene2=pp.split('_')
    g1obj,g2obj=genedb[gene1],genedb[gene2]
    fus_disc=[]									#collecting discordant pairs
    for rdname in discordant[pp]:
        #arrange read1/read2 into F/R so it will be easier for GeneFusion obj to handle
        r1,r2=db1[rdname],db2[rdname]
        orient=discordant[pp][rdname]
        if orient=='F:R':
            fus_disc.append((r1,r2))
        elif orient=='R:F':
            fus_disc.append((r2,r1))
    fus_jsr=[]    
    if junc_map.has_key(pp):
        for rdname in junc_map[pp]:
            r=rdb[rdname]['read']
            junc=rdb[rdname]['junc']
            jsr=gfclass.JSR(r,junc)
            fus_jsr.append(jsr)
    gf=gfclass.GeneFusion(gene1,gene2,fus_disc,fus_jsr)
    gf_new=gf.update(mm=mm)     ##apply the mismatch parameter, default is 1
    #output the results
    disc_n=str(len(gf_new.discordantpairs))
    junctions=sorted(gf_new.get_junction_freq(),key=operator.itemgetter(1),reverse=True) #sort junc by # of JSRs
    junc_n=str(len(junctions))
    junc_str='|'.join([','.join([x[0],str(x[1])]) for x in junctions])
    jsr_n=str(len(gf_new.fusionreads))
    pjsr_n=str(len(gf_new.get_perfect_JSR()))
    pos_consist=gf_new.positioncheck()
    svec=[gene1,gene2,g1obj.chr,g2obj.chr,g1obj.strand,g2obj.strand,disc_n,jsr_n,pjsr_n,junc_n,pos_consist,junc_str]
    outfile_s.write('%s\n'%('\t'.join(svec)))
    outfile_d.write('@@\t%s,%s,%s\t%s,%s,%s\n'%(gene1,g1obj.chr,g1obj.strand,gene2,g2obj.chr,g2obj.strand))
    outfile_d.write('\n')
    outfile_d.write('>discordant\n')
    for rp in gf_new.discordantpairs:
        rf,rr=rp
        nm1=[x[1] for x in rf.tags if x[0]=='NM'][0]
        nm2=[x[1] for x in rr.tags if x[0]=='NM'][0]    
        outfile_d.write('%s\tF\t%s.%s.mm%d\n'%(rf.qname,gene1,rf.pos+1,nm1))    ##0-based coordinates
        outfile_d.write('%s\tR\t%s.%s.mm%d\n'%(rr.qname,gene2,rr.pos+1,nm2))    ##0-based coordinates
    outfile_d.write('\n')
    outfile_d.write('>spanning\n') 
    for jsr in gf_new.fusionreads:
        r=jsr.read
        nm=[x[1] for x in r.tags if x[0]=='NM'][0]
        outfile_d.write('%s\t%s.mm%d\n'%(r.qname,jsr.junction,nm))
    outfile_d.write('\n')
    outfile_d.write('>junction\n')
    for junc_info in junctions:
        outfile_d.write('%s\t%d\n'%(junc_info[0],junc_info[1]))
    outfile_d.write('\n')
    outfile_d.write('>summary\n')
    outfile_d.write('Number of Discordant Pairs = %s\n'%disc_n)
    outfile_d.write('Number of Fusion Reads = %s\n'%jsr_n)
    outfile_d.write('Number of Perfect Fusion Reads = %s\n'%pjsr_n)
    outfile_d.write('Number of Distinct Junctions = %s\n'%junc_n)
    outfile_d.write('Position Consistency = %s\n'%pos_consist)
    outfile_d.write('\n')

outfile_s.close()
outfile_d.close() 

########################################################################################
print 'step 7: generating fusion lists @ %s'%time.ctime()

#For convenience, filter the lists to candidates with 
# 1) at least 2 discordant pairs
# 2) at least 1 perfect JSR
#meanwhile, calculate sequence similarity for each pair
#user may need to manually filter the lists per this measure.

#The following code is a copy of prada-homology
outfile_o=open('%s/%s.fus.summary.txt'%(outpath,docstring),'w')
ifname='%s/%s.fus.candidates.txt'%(outpath,docstring) 
if not os.path.exists(ifname):
    sys.exit('ERROR: %s was not found'%ifname)

blastseq_tmp_dir='%s/blast_tmp/'%outpath
if not os.path.exists(blastseq_tmp_dir):
    os.mkdir(blastseq_tmp_dir)

flists=[]
infile=open(ifname)
iN=0
for line in open(ifname):
    info=line.strip().split('\t')
    if iN==0:
        iN+=1		#skip title
        flists.append(info)
        continue
    else:
        if int(info[6])>=2 and int(info[8])>=1 and info[10] in ['PARTIALLY','YES']:
            flists.append(info) 
infile.close()

if len(flists)==1:  						#if no candidate passes the filters
    outfile_o.write('%s\n'%'\t'.join(flists[0]))
    outfile_o.close()
    print 'step done @ %s'%time.ctime()
    sys.exit(0)

candidates={}
for line in flists[1:]:
    geneA,geneB=line[0],line[1]
    key='%s_%s'%(geneA,geneB)
    candidates[key]=''

selecttranscript={}
for gene in genedb:
    txs=genedb[gene].transcript
    stx=txs[0]
    initlen=stx.length
    for tx in txs:
        if tx.length >= initlen:
            stx=tx
            initlen=stx.length
    selecttranscript[gene]=stx.name

allpartners=set()
for item in candidates:
    sset=set(item.split('_'))
    allpartners=allpartners.union(sset)

presenttxs=[]    #tx that is present in our annotation
absent=[]     #tx that is not in our annotation
for gene in allpartners:
    if selecttranscript.has_key(gene):
        presenttxs.append(selecttranscript[gene])
    else:
        absent.append(gene)

for seq_record in SeqIO.parse(txseqfile,'fasta'):
    sid=seq_record.id
    seq=seq_record.seq
    if sid in presenttxs:
        g=txdb[sid].gene
        fastafile=open('%s/%s.fasta'%(blastseq_tmp_dir,g),'w')
        SeqIO.write(seq_record,fastafile,'fasta')
        fastafile.close()

for gp in candidates:
    geneA,geneB=gp.split('_')
    if geneA in absent or geneB in absent:
        candidates[gp]=['NA']*4
    else:
        gaseq='%s/%s.fasta'%(blastseq_tmp_dir,geneA)
        gaobj=SeqIO.parse(gaseq,'fasta').next()
        gbseq='%s/%s.fasta'%(blastseq_tmp_dir,geneB)
        gbobj=SeqIO.parse(gbseq,'fasta').next()
        ga_len,gb_len=str(len(gaobj.seq)),str(len(gbobj.seq))
        a=privutils.seqblast(gaseq,gbseq,blastn)
        if a==None:
            candidates[gp]=['NA','NA','100','0']
        else:
            candidates[gp]=a

header=flists[0][:]
header.extend(['Identity','Align_Len','Evalue','BitScore'])
outfile_o.write('%s\n'%('\t'.join(header)))

for info in flists[1:]:
    geneA,geneB=info[0],info[1]
    key='%s_%s'%(geneA,geneB)
    vv=candidates[key]
    row=info[:]
    row.extend(vv)
    outfile_o.write('%s\n'%('\t'.join(row)))
outfile_o.close()

########################################################################################
print 'step done @ %s'%time.ctime()