POSCO/plotGoAnalysis_final.py at main · hendrixlab/POSCO · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import re,sys,os
import copy
import math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import csv
from matplotlib.colors import LinearSegmentedColormap


#plt.rcParams.update({"text.usetex": True})
def addLine(ax,x1,x2,y1,y2,lw):
    line = plt.Line2D([x1,x2],[y1,y2], color='black',lw = lw)
    line.set_clip_on(False)
    ax.add_line(line)

def readCSV(csvFile):

    data = []
    with open(csvFile) as c:
        csvReader = csv.reader(c, delimiter = '\t')
        #for x in range(0,12): next(csvReader)
        for x in range(0,12):
            next(csvReader)
        for row in csvReader:
            #print(row)
            sig = -1 *  math.log(float(row[-1]),10)

            data.append((row[0],float(row[2]),float(row[3]),row[4],float(row[-1]), sig))

    # sort in ascending order, from least to most significant
    data.sort(key=lambda x:x[-1])
    #index 1 : observed (y)
    #index 2 : expected (x)
    #index 3 : significance = p-value
    #index 4: -log(p-value)
    for d in data: print(d)
    print(len(data))
    return data


def plotScatter(data):

   fig,ax = plt.subplots()

   labels,yVals,xVals,folds,pVals,clrs = zip(*data)
   #colors = ['midnightblue','lime']
   colors = ['darkred','gold']
   cm = LinearSegmentedColormap.from_list('dah',colors, N=100)

   maxX = max(xVals)
   maxY = int(max(yVals) + 10000)
   lineData = [x for x in range(-1,int(maxY)+1)]
   print(max(yVals))
   stop = math.log(max(yVals) - 600,10)

   # plot the line y=x
   plt.plot(lineData,lineData,color = "grey",alpha = 0.5)

   # plot the scatter of points
   plt.scatter(xVals,yVals,c = clrs,cmap = cm, edgecolor = 'none',alpha=0.8, linewidths = 5)
   plt.xlabel("Expected Count")
   plt.xscale('symlog')
   plt.yscale('symlog')
   ax.set_xlim(-1,maxY)
   ax.set_ylim(-1,maxY)
   plt.ylabel("Observed Count")
   plt.colorbar(label = r'$-\log_{10} FDR$',orientation="horizontal",pad=0.2)
   plt.title("GO Terms for Lowest POSCO Quintile", loc='left')
   #plt.title("GO Terms for Highest POSCO Quintile", loc='left')
   # plot the labels
   numToPlot = 16
   topData = data[-numToPlot:]
   topData.sort(key=lambda x: x[2])
   for d in topData: print(d)
   logLabels = np.logspace(0,4,num = numToPlot)
   xLogLabels = np.logspace(2,6,num = numToPlot)
   n = 0
   offset = 0.35
   count = 0
   for i in range(len(topData)):
       l,y,x,f,p,s = topData[i]
       ll = logLabels[count]
       xll = xLogLabels[count]
       label = str(l.split('(')[0])
       n += 1
       uc = ''
       print(f)
       if 'biological_process' not in label and 'Unclassified' not in label:
           count += 1
           if f == '+':
               plt.text(xll,ll,u'\u25b2' + " " +label)
           else:
               plt.text(xll,ll,u'\u25bc' + " " +label)

           addLine(ax,x,xll,y,ll,0.2)
   ax.set_xlim(-1,math.pow(10,7.5))
   plt.text(0.5,1.0,"y = x",rotation = 38)
   ax.spines['right'].set_visible(False)
   ax.spines['top'].set_visible(False)
   plt.savefig(f"goTermsScatter_{suffix}.pdf",bbox_inches = "tight")
   #print(n)

########
# main #
########

usage = "Usage: " + sys.argv[0] + " <GO Term CSV>"
if len(sys.argv) != 3:
    print(usage)
    sys.exit()

csvFile = sys.argv[1]
suffix = sys.argv[2]
sortedData = readCSV(csvFile)
#for s in sortedData : print(s)
plotScatter(sortedData)