comparison defuse_trinity_analysis.py @ 39:90127ee1eae5

Fix defuse_trinity_analysis.py

38:9c815a3721b3

      fusion_type = 'inter' if fusion['gene_chromosome1'] != fusion['gene_chromosome2'] else 'intra' if abs(fusion['genomic_break_pos1'] - fusion['genomic_break_pos2']) > options.ticdist else 'TIC'
      fusion['fusion_type'] = fusion_type
      fusion['transcripts'] = []
      fusion['Transcript'] = 'No'
      fusion['Protein'] = 'No'
      print >> sys.stdout, "%4d\t%6s\t%s\t%s\t%s\t%s\t%s" % (i,fusion['cluster_id'],fwd_seq,rev_seq,fusion_type,fusion['gene_name1'],fusion['gene_name2'])
  inputFile.close()

  ## Process Trinity data and compare to deFuse
  matched_transcripts = dict()
  matched_orfs = dict()

          if fusion['fwd_seq'] in seq or fusion['rev_seq'] in seq:
            fusions_with_transcripts.add(i)
            matched_transcripts[name] = seq
            fusion['transcripts'].append(name)
            fusion['Transcript'] = 'Yes'
    print >> sys.stdout, "fusions_with_transcripts: %d  %s\n matched_transcripts: %d" % (len(fusions_with_transcripts),fusions_with_transcripts,len(matched_transcripts))
    for i,fusion in enumerate(fusions):
      print >> sys.stdout, "%4d\t%6s\t%s\t%s\t%s\t%s\t%s\t%s" % (i,fusion['cluster_id'],fusion['fwd_seq'],fusion['rev_seq'],fusion['fusion_type'],fusion['gene_name1'],fusion['gene_name2'], fusion['transcripts'])
    ## Process ORFs and compare to matched deFuse and Trinity data.
    ## Proteins must be at least 100 aa long, starting at the first "M" and must end with an "*".
    if options.peptides: 
      with open(options.peptides) as fp:
        for name, seq in read_fasta(fp):

                  if pep_len - start < options.min_pep_len:
                    continue
                  fusions_with_orfs.add(i)
                  matched_orfs[name] = seq
                  fusion['Protein'] = 'Yes'
                  # fwd or reverse
                  tx_seq = matched_transcripts(tx_id)
                  pos = tx_seq.find(fusion['fwd_seq'])
                  if pos < 0:
                    pos = tx_seq.find(fusion['rev_seq'])
                  # locate fusion in transcript
                  # locate fusion in ORF
                  fusion['prior_pep_seq'] = ''
                  fusion['novel_pep_seq'] = ''
      print >> sys.stdout, "fusions_with_orfs: %d  %s\n matched_orfs: %d" % (len(fusions_with_orfs),fusions_with_orfs,len(matched_orfs))
  ## Write reports
  report_fields = ['gene_name1','gene_name2','span_count','probability','gene_chromosome1','gene_location1','gene_chromosome2','gene_location2','fusion_type','Transcript','Protein']
  report_colnames = {'gene_name1':'Gene 1','gene_name2':'Gene 2','span_count':'Span cnt','probability':'Probability','gene_chromosome1':'From Chr','gene_location1':'Fusion point','gene_chromosome2':'To Chr','gene_location2':'Fusion point','fusion_type':'Type','Transcript':'Transcript?','Protein':'Protein?' }
  print >> outputFile,"%s\t%s" % ('#','\t'.join([report_colnames[x] for x in report_fields]))
  for i,fusion in enumerate(fusions): 

39:90127ee1eae5

      fusion_type = 'inter' if fusion['gene_chromosome1'] != fusion['gene_chromosome2'] else 'intra' if abs(fusion['genomic_break_pos1'] - fusion['genomic_break_pos2']) > options.ticdist else 'TIC'
      fusion['fusion_type'] = fusion_type
      fusion['transcripts'] = []
      fusion['Transcript'] = 'No'
      fusion['Protein'] = 'No'
      #print >> sys.stdout, "%4d\t%6s\t%s\t%s\t%s\t%s\t%s" % (i,fusion['cluster_id'],fwd_seq,rev_seq,fusion_type,fusion['gene_name1'],fusion['gene_name2'])
  inputFile.close()

  ## Process Trinity data and compare to deFuse
  matched_transcripts = dict()
  matched_orfs = dict()

          if fusion['fwd_seq'] in seq or fusion['rev_seq'] in seq:
            fusions_with_transcripts.add(i)
            matched_transcripts[name] = seq
            fusion['transcripts'].append(name)
            fusion['Transcript'] = 'Yes'
    #print >> sys.stdout, "fusions_with_transcripts: %d  %s\n matched_transcripts: %d" % (len(fusions_with_transcripts),fusions_with_transcripts,len(matched_transcripts))
    print >> sys.stdout, "fusions_with_transcripts: %d unique_transcripts: %d" % (len(fusions_with_transcripts),len(matched_transcripts))
    #for i,fusion in enumerate(fusions):
    #  print >> sys.stdout, "%4d\t%6s\t%s\t%s\t%s\t%s\t%s\t%s" % (i,fusion['cluster_id'],fusion['fwd_seq'],fusion['rev_seq'],fusion['fusion_type'],fusion['gene_name1'],fusion['gene_name2'], fusion['transcripts'])
    ## Process ORFs and compare to matched deFuse and Trinity data.
    ## Proteins must be at least 100 aa long, starting at the first "M" and must end with an "*".
    if options.peptides: 
      with open(options.peptides) as fp:
        for name, seq in read_fasta(fp):

                  if pep_len - start < options.min_pep_len:
                    continue
                  fusions_with_orfs.add(i)
                  matched_orfs[name] = seq
                  fusion['Protein'] = 'Yes'
                  """
                  # fwd or reverse
                  tx_seq = matched_transcripts(tx_id)
                  pos = tx_seq.find(fusion['fwd_seq'])
                  if pos < 0:
                    pos = tx_seq.find(fusion['rev_seq'])
                  # locate fusion in transcript
                  # locate fusion in ORF
                  fusion['prior_pep_seq'] = ''
                  fusion['novel_pep_seq'] = ''
                  """
      #print >> sys.stdout, "fusions_with_orfs: %d  %s\n matched_orfs: %d" % (len(fusions_with_orfs),fusions_with_orfs,len(matched_orfs))
      print >> sys.stdout, "fusions_with_orfs: %d  unique_orfs: %d" % (len(fusions_with_orfs),len(matched_orfs))
  ## Write reports
  report_fields = ['gene_name1','gene_name2','span_count','probability','gene_chromosome1','gene_location1','gene_chromosome2','gene_location2','fusion_type','Transcript','Protein']
  report_colnames = {'gene_name1':'Gene 1','gene_name2':'Gene 2','span_count':'Span cnt','probability':'Probability','gene_chromosome1':'From Chr','gene_location1':'Fusion point','gene_chromosome2':'To Chr','gene_location2':'Fusion point','fusion_type':'Type','Transcript':'Transcript?','Protein':'Protein?' }
  print >> outputFile,"%s\t%s" % ('#','\t'.join([report_colnames[x] for x in report_fields]))
  for i,fusion in enumerate(fusions):