# HG changeset patch # User cmmt # Date 1327314260 18000 # Node ID 519e45d9539fcedc2bde31196e474b77b6c04319 # Parent 116d726d2ed49cdb0b0d8d1b391aa6be37f2a768 Uploaded diff -r 116d726d2ed4 -r 519e45d9539f chromatral.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/chromatral.py Mon Jan 23 05:24:20 2012 -0500 @@ -0,0 +1,283 @@ +import os; +import re; +import csv; +import sys; +import math; +import numpy as np; +import matplotlib; +matplotlib.use('Agg'); +from pylab import *; + +def extractGffFeatures(gffString): + # GFF v3 + # 1. seqname - Must be a chromosome or scaffold. + # 2. source - The program that generated this feature. + # 3. feature - The name of this type of feature. Some examples of standard feature types are "CDS", "start_codon", "stop_codon", and "exon". + # 4. start - The starting position of the feature in the sequence. The first base is numbered 1. + # 5. end - The ending position of the feature (inclusive). + # 6. score - A score between 0 and 1000. If there is no score value, enter ".". + # 7. strand - Valid entries include '+', '-', or '.' (for don't know/care). + # 8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'. + # 9. group - All lines with the same group are linked together into a single item. + + res = []; + + if len(gffString) > 1: # check for empty line (i.e. row in GFF file only contains \n) + try: + gffSeqName, gffSource, gffFeature, gffStart, gffEnd, gffScr, gffStrand, gffFrame, gffGroup = gffString.split('\t'); + except: + return []; # exit, since we don't have 9 features + + # handle the gffSeqName + try: + chrNo = int(re.findall(r'[0-9]+', gffSeqName)[0]); + except: + return []; + + # handle gffStart + try: + probeStart = int(gffStart); + except ValueError: + return []; # exit, invalid gffStart + + # handle gffEnd + try: + probeEnd = int(gffEnd); + except ValueError: + return []; # exit, invalid gffEnd + + # handle gffScr + try: + probeScr = float(gffScr); + except ValueError: + return []; # exit, invalid gffScr + + # if the input string was well-formed, return the extracted features + res = [chrNo, probeStart, probeEnd, probeScr]; + return res; # return feature set (if any was found in this record) + + +# Define Constants +tickPlotAdj = 0.5; # adjust tick marks on x-axis such that they are exactly at the bin boundary (and not in the bin center) + +# Initialize Variables +feaNames = []; +feaChrNoList = []; +feaStrandList = []; +feaStartList = []; +feaEndList = []; + +gffChrNoList = []; +gffProbeCenterList = []; +gffMatScrList = []; + +# get command line arguments +gffFileName = sys.argv[1]; +feaFileName = sys.argv[2]; +outFileFormat = sys.argv[3]; +stepWidth = int(sys.argv[4]); +upStreamStretch = int(sys.argv[5]); # in bp +maxX = int(sys.argv[6]); # in bp, if -1 maxX <= max(feature length), see below +tickSpacing = int(sys.argv[7]); # in bp +plotTitle = sys.argv[8]; +outFileName = sys.argv[9]; + +# Read transcript information +feaFileReader = csv.reader(open(feaFileName, 'rb'), delimiter = '\t'); +for row in feaFileReader: + feaName, feaChrNo, feaStrand, feaStart, feaEnd = row; + feaNames.append(feaName); + feaChrNoList.append(int(feaChrNo)); + feaStrandList.append(int(feaStrand)); + feaStartList.append(int(feaStart)); + feaEndList.append(int(feaEnd)); + +# Read enrichment scores +gffFileHandle = open(gffFileName, "rU"); +for row in gffFileHandle: + cache = extractGffFeatures(row); + if len(cache) == 4: # if curr row wasn't well-formed, skip row + chrNo, probeStart, probeEnd, probeScr = cache; + gffChrNoList.append(chrNo); + gffProbeCenterList.append(int(float(probeEnd - probeStart)/2.0 + probeStart)); # calc the center pos of each probe + gffMatScrList.append(probeScr); + +# convert numerical lists into vectors +feaChrNos = np.array(feaChrNoList); +feaStrands = np.array(feaStrandList); +feaStarts = np.array(feaStartList); +feaEnds = np.array(feaEndList); + +gfChrNos = np.array(gffChrNoList); +gfProbeLocs = np.array(gffProbeCenterList); +gfMatScrs = np.array(gffMatScrList); + +# preallocate the results array, bins and 1 extra cols for feature length +avgFeatureMatScrs = np.empty((len(feaStarts),\ + (int(math.ceil(float(max(feaEnds - feaStarts))/float(stepWidth)) + upStreamStretch/stepWidth + 1))),\ + ); +avgFeatureMatScrs[:] = np.NAN; + +featureCounter = 0; + +for chr in range(min(feaChrNos),max(feaChrNos) + 1): + #get all rows from the feature file that contain information about the i'th chromosome + currChrIdxs = (feaChrNos == chr); + currStrands = feaStrands[currChrIdxs]; + currStarts = feaStarts[currChrIdxs]; + currEnds = feaEnds[currChrIdxs]; + + #get all rows from ChIP-* data sets + currGfChrIdxs = (gfChrNos == chr); + currGfProbeLocs = gfProbeLocs[currGfChrIdxs]; + currGfScrs = gfMatScrs[currGfChrIdxs]; + + # iterate over all features of the current chromosome + for j in range(0,len(currStarts)): + feaStrand = currStrands[j]; + feaStart = currStarts[j]; + feaEnd = currEnds[j]; + + #calc the corresponding probe positions for the feature start + if feaStrand == 1: # +strand + refFeaSt = 0; + while feaStart >= currGfProbeLocs[refFeaSt]: + refFeaSt += 1; + else: # -strand + refFeaSt = 0; + while feaEnd >= currGfProbeLocs[refFeaSt]: + refFeaSt += 1; + #now refFeaSt points to the correct probe according to strand orientation + + #calc the avg probe scores per bin + if feaStrand == 1: + feaChunkStart = feaStart; # bp level pointer + refChunkStart = refFeaSt; # probe level pointer + refChunkEnd = refFeaSt; # probe level pointer + + # determine values for the upstream bins + for k in range(1,upStreamStretch/stepWidth + 1): + feaChunkEnd = feaStart - k * stepWidth; # bp level pointer + + # find the correct ref for the current feaChunkEnd; catch potential IndexOutOfBoundEx when upstream region too close to chr boundary + try: + while feaChunkEnd <= currGfProbeLocs[refChunkEnd]: + refChunkEnd -= 1; + except: + pass; # don't handle the exception further + + #calc probe mean + avgFeatureMatScrs[featureCounter, upStreamStretch/stepWidth - k] = np.mean(currGfScrs[refChunkEnd:refChunkStart + 1]); + + #adjust pointers + refChunkStart = refChunkEnd; + + # determine values for bins in feature body + for k in range(1, int(math.ceil(float(feaEnd - feaStart)/float(stepWidth)) + 1)): + feaChunkEnd = feaStart + k*stepWidth; # bp level pointer + + try: + while feaChunkEnd >= currGfProbeLocs[refChunkEnd]: + refChunkEnd += 1; + except: + pass; # don't handle the exception further + + # calc probe mean + avgFeatureMatScrs[featureCounter, k - 1 + upStreamStretch/stepWidth] = np.mean(currGfScrs[refChunkStart:refChunkEnd + 1]); + + # adjust chunk pointers + refChunkStart = refChunkEnd; + + # write feature length in last col of the array + avgFeatureMatScrs[featureCounter, -1] = feaEnd - feaStart; + featureCounter += 1; + else: # -1 strand + feaChunkStart = feaEnd; # bp level pointer + refChunkStart = refFeaSt; # probe level pointer + refChunkEnd = refFeaSt; # probe level pointer + + # determine values for the upStream region bins + for k in range(1,upStreamStretch/stepWidth + 1): + feaChunkEnd = feaEnd + k*stepWidth; # bp level pointer + + # find the ref for the current feaChunkEnd, catch potential IndexOutOfBoundEx when upstream region too close to chr boundary + try: + while feaChunkEnd >= currGfProbeLocs[refChunkEnd]: + refChunkEnd += 1; + except: + pass; + + # calc probe mean + avgFeatureMatScrs[featureCounter, upStreamStretch/stepWidth - k] = np.mean(currGfScrs[refChunkStart:refChunkEnd+1]); + + # adjust chunk pointers + refChunkStart = refChunkEnd; + + # determine values for bins in feature body + for k in range(1, int(math.ceil(float(feaEnd - feaStart)/float(stepWidth)) + 1)): + feaChunkEnd = feaEnd - k*stepWidth; # bp level pointer + + try: + while feaChunkEnd <= currGfProbeLocs[refChunkEnd]: + refChunkEnd -= 1; + except: + pass; + + # calc probe mean + avgFeatureMatScrs[featureCounter, k-1 + upStreamStretch/stepWidth] = np.mean(currGfScrs[refChunkEnd:refChunkStart+1]); + + # adjust chunk pointers + refChunkStart = refChunkEnd; + + # write feature length in last col of the array + avgFeatureMatScrs[featureCounter, -1] = feaEnd - feaStart; + featureCounter += 1; + + +# assign names and data format to all cols of the results array and sort array by feature length (last column) +sortByLengthCol = 'col' + str(np.size(avgFeatureMatScrs,1)-1); +matDtype = {'names':['col%i'%i for i in range(np.size(avgFeatureMatScrs,1))],'formats':(np.size(avgFeatureMatScrs,1))*[np.float]}; +avgFeatureMatScrs.dtype = matDtype; +avgFeatureMatScrs.sort(axis = 0, order = sortByLengthCol); + +# generate colormap +BYcmapData = {'red': ((0.0, 0.0, 0.094), (0.5, 0.0, 0.0), (1.0, 0.992, 1.0)), 'green': ((0.0, 0.0,0.658), (0.5, 0.0, 0.0), (1.0, 0.996, 1.0)), 'blue': ((0.0, 1.0, 0.828), (0.5, 0.0, 0.0), (1.0, 0.0, 0.0))}; +BYcmap = matplotlib.colors.LinearSegmentedColormap('BYcmap', BYcmapData,9); + +lowerHinge = matplotlib.mlab.prctile(gfMatScrs, 25.0); +upperHinge = matplotlib.mlab.prctile(gfMatScrs, 75.0); +cmLim = (upperHinge - lowerHinge) * 1.5 + upperHinge; # upper inner fence + +# remove custom array dtype used for sorting and plot avgScrs (without the last col which stores the length info) +imgplot = imshow(avgFeatureMatScrs.view(np.float)[:,0:-1], aspect='auto', interpolation='nearest', origin='lower', vmin = -cmLim, vmax = cmLim); +axvline(x = int(float(upStreamStretch)/float(stepWidth)) - tickPlotAdj, color = 'w', linewidth = 1); + +# if maxX wasn't set by hand, set it to max feature length and set x-coord accordingly +if maxX == -1: + maxX = max(feaEnds - feaStarts); + imgplot.axes.set_xbound(-tickPlotAdj, int(float(maxX)/float(stepWidth)) + int(float(upStreamStretch)/float(stepWidth)) - tickPlotAdj); +else: + imgplot.axes.set_xbound(-tickPlotAdj, int(float(maxX)/float(stepWidth)) - tickPlotAdj); + +xTicks = [k * int(float(tickSpacing)/float(stepWidth)) + int(float(upStreamStretch)/float(stepWidth)) - tickPlotAdj for k in range(0, int(math.ceil(float(maxX)/float(tickSpacing)))+1)]; +xTicks.insert(-1,-tickPlotAdj); +xTicksLabel = [k * tickSpacing for k in range(0, int(math.ceil(float(maxX)/float(tickSpacing)))+1)]; +xTicksLabel.insert(-1, str(-upStreamStretch)); +imgplot.axes.set_xticks(xTicks); +imgplot.axes.set_xticklabels(xTicksLabel); + +imgplot.axes.set_title(plotTitle); +imgplot.axes.set_xlabel('Distance from feature start (bp)'); +imgplot.axes.set_ylabel(''); +imgplot.axes.set_yticks([]); + +imgplot.set_cmap(BYcmap); +colorbar(); +show(); + +savefig('chromatra_l_tmp_plot', format=outFileFormat); +data = file('chromatra_l_tmp_plot', 'rb').read(); +fp = open(outFileName, 'wb'); +fp.write(data); +fp.close(); +os.remove('chromatra_l_tmp_plot');