/tools/solid_tools/solid_qual_stats.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 140 lines · 114 code · 24 blank · 2 comment · 41 complexity · 5560d37474ff2b3d8ebf369c696fb229 MD5 · raw file

  1. #!/usr/bin/env python
  2. #Guruprasad Ananda
  3. import sys, os, zipfile, tempfile
  4. QUAL_UPPER_BOUND = 41
  5. QUAL_LOWER_BOUND = 1
  6. def stop_err( msg ):
  7. sys.stderr.write( "%s\n" % msg )
  8. sys.exit()
  9. def unzip( filename ):
  10. zip_file = zipfile.ZipFile( filename, 'r' )
  11. tmpfilename = tempfile.NamedTemporaryFile().name
  12. for name in zip_file.namelist():
  13. file( tmpfilename, 'a' ).write( zip_file.read( name ) )
  14. zip_file.close()
  15. return tmpfilename
  16. def __main__():
  17. infile_score_name = sys.argv[1].strip()
  18. fout = open(sys.argv[2].strip(),'r+w')
  19. infile_is_zipped = False
  20. if zipfile.is_zipfile( infile_score_name ):
  21. infile_is_zipped = True
  22. infile_name = unzip( infile_score_name )
  23. else:
  24. infile_name = infile_score_name
  25. readlen = None
  26. invalid_lines = 0
  27. j = 0
  28. for line in file( infile_name ):
  29. line = line.strip()
  30. if not(line) or line.startswith("#") or line.startswith(">"):
  31. continue
  32. elems = line.split()
  33. try:
  34. for item in elems:
  35. int(item)
  36. if not readlen:
  37. readlen = len(elems)
  38. if len(elems) != readlen:
  39. print "Note: Reads in the input dataset are of variable lengths."
  40. j += 1
  41. except ValueError:
  42. invalid_lines += 1
  43. if j > 10:
  44. break
  45. position_dict = {}
  46. print >>fout, "column\tcount\tmin\tmax\tsum\tmean\tQ1\tmed\tQ3\tIQR\tlW\trW"
  47. for k,line in enumerate(file( infile_name )):
  48. line = line.strip()
  49. if not(line) or line.startswith("#") or line.startswith(">"):
  50. continue
  51. elems = line.split()
  52. if position_dict == {}:
  53. for pos in range(readlen):
  54. position_dict[pos] = [0]*QUAL_UPPER_BOUND
  55. if len(elems) != readlen:
  56. invalid_lines += 1
  57. continue
  58. for ind,item in enumerate(elems):
  59. try:
  60. item = int(item)
  61. position_dict[ind][item]+=1
  62. except:
  63. pass
  64. invalid_positions = 0
  65. for pos in position_dict:
  66. carr = position_dict[pos] #count array for position pos
  67. total = sum(carr) #number of bases found in this column.
  68. med_elem = int(round(total/2.0))
  69. lowest = None #Lowest quality score value found in this column.
  70. highest = None #Highest quality score value found in this column.
  71. median = None #Median quality score value found in this column.
  72. qsum = 0.0 #Sum of quality score values for this column.
  73. q1 = None #1st quartile quality score.
  74. q3 = None #3rd quartile quality score.
  75. q1_elem = int(round((total+1)/4.0))
  76. q3_elem = int(round((total+1)*3/4.0))
  77. try:
  78. for ind,cnt in enumerate(carr):
  79. qsum += ind*cnt
  80. if cnt!=0:
  81. highest = ind
  82. if lowest==None and cnt!=0: #first non-zero count
  83. lowest = ind
  84. if q1==None:
  85. if sum(carr[:ind+1]) >= q1_elem:
  86. q1 = ind
  87. if median==None:
  88. if sum(carr[:ind+1]) < med_elem:
  89. continue
  90. median = ind
  91. if total%2 == 0: #even number of elements
  92. median2 = median
  93. if sum(carr[:ind+1]) < med_elem+1:
  94. for ind2,elem in enumerate(carr[ind+1:]):
  95. if elem != 0:
  96. median2 = ind+ind2+1
  97. break
  98. median = (median + median2)/2.0
  99. if q3==None:
  100. if sum(carr[:ind+1]) >= q3_elem:
  101. q3 = ind
  102. mean = qsum/total #Mean quality score value for this column.
  103. iqr = q3-q1
  104. left_whisker = max(q1 - 1.5*iqr,lowest)
  105. right_whisker = min(q3 + 1.5*iqr,highest)
  106. print >>fout,"%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" %(pos+1,total,lowest,highest,qsum,mean,q1,median,q3,iqr,left_whisker,right_whisker)
  107. except:
  108. invalid_positions += 1
  109. nullvals = ['NA']*11
  110. print >>fout,"%s\t%s" %(pos+1,'\t'.join(nullvals))
  111. if invalid_lines:
  112. print "Skipped %d reads as invalid." %invalid_lines
  113. if invalid_positions:
  114. print "Skipped stats computation for %d read positions." %invalid_positions
  115. if __name__=="__main__":
  116. __main__()