/tools/ngs_rna/filter_transcripts_via_tracking.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 70 lines · 38 code · 11 blank · 21 comment · 9 complexity · 34756a9ee27636884e723c75347a3633 MD5 · raw file

  1. #!/usr/bin/env python
  2. import os, sys, tempfile
  3. assert sys.version_info[:2] >= ( 2, 4 )
  4. def __main__():
  5. """
  6. Utility script for analyzing Cufflinks data: uses a tracking file (produced by cuffcompare) to filter a GTF file of transcripts (usually the transcripts
  7. produced by cufflinks). Filtering is done by extracting transcript IDs from tracking file and then filtering the GTF so that the output GTF contains only
  8. transcript found in the tracking file. Because a tracking file has multiple samples, a sample number is used to filter transcripts for
  9. a particular sample.
  10. """
  11. # Read parms.
  12. tracking_file_name = sys.argv[1]
  13. transcripts_file_name = sys.argv[2]
  14. output_file_name = sys.argv[3]
  15. sample_number = int ( sys.argv[4] )
  16. # Open files.
  17. transcripts_file = open( transcripts_file_name, 'r' )
  18. output_file = open( output_file_name, 'w' )
  19. # Read transcript IDs from tracking file.
  20. transcript_ids = {}
  21. for i, line in enumerate( file( tracking_file_name ) ) :
  22. # Split line into elements. Line format is
  23. # [Transfrag ID] [Locus ID] [Ref Gene ID] [Ref Transcript ID] [Class code] [qJ:<gene_id>|<transcript_id>|<FMI>|<FPKM>|<conf_lo>|<conf_hi>]
  24. line = line.rstrip( '\r\n' )
  25. elems = line.split( '\t' )
  26. # Get transcript info.
  27. if sample_number == 1:
  28. transcript_info = elems[4]
  29. elif sample_number == 2:
  30. transcript_info = elems[5]
  31. if not transcript_info.startswith('q'):
  32. # No transcript for this sample.
  33. continue
  34. # Get and store transcript id.
  35. transcript_id = transcript_info.split('|')[1]
  36. transcript_id = transcript_id.strip('"')
  37. transcript_ids[transcript_id] = ""
  38. # Filter transcripts file using transcript_ids
  39. for i, line in enumerate( file( transcripts_file_name ) ):
  40. # GTF format: chrom source, name, chromStart, chromEnd, score, strand, frame, attributes.
  41. elems = line.split( '\t' )
  42. # Get attributes.
  43. attributes_list = elems[8].split(";")
  44. attributes = {}
  45. for name_value_pair in attributes_list:
  46. pair = name_value_pair.strip().split(" ")
  47. name = pair[0].strip()
  48. if name == '':
  49. continue
  50. # Need to strip double quote from values
  51. value = pair[1].strip(" \"")
  52. attributes[name] = value
  53. # Get element's transcript id.
  54. transcript_id = attributes['transcript_id']
  55. if transcript_id in transcript_ids:
  56. output_file.write(line)
  57. # Clean up.
  58. output_file.close()
  59. if __name__ == "__main__": __main__()