Computes UnrootedGraph.expr(), UnrootedGraph.delta(), UnrootedGraph.connectivity(), UnrootedGraph.connectivity_pvalue() and Benjamini-Holchberg adjusted q-values for all genes that are expressed in more than 'filtercells' cells and whose maximum expression value is a
(self, n=500, filtercells=0, filterexp=0.0, annotation={})
| 1048 | return pi |
| 1049 | |
| 1050 | def save(self, n=500, filtercells=0, filterexp=0.0, annotation={}): |
| 1051 | """ |
| 1052 | Computes UnrootedGraph.expr(), UnrootedGraph.delta(), UnrootedGraph.connectivity(), |
| 1053 | UnrootedGraph.connectivity_pvalue() and Benjamini-Holchberg adjusted q-values for all |
| 1054 | genes that are expressed in more than 'filtercells' cells and whose maximum expression |
| 1055 | value is above 'filterexp'. The optional argument 'annotation' allos to include a dictionary |
| 1056 | with lists of genes to be annotated in the table. The output is stored in a tab separated |
| 1057 | file called name.genes.txt. |
| 1058 | """ |
| 1059 | pol = [] |
| 1060 | with open(self.name + '.genes.tsv', 'w') as ggg: |
| 1061 | cul = 'Gene\tCells\tMean\tMin\tMax\tConnectivity\tp_value\tq-value (BH)\t' |
| 1062 | for m in sorted(annotation.keys()): |
| 1063 | cul += m + '\t' |
| 1064 | ggg.write(cul[:-1] + '\n') |
| 1065 | lp = sorted(self.dicgenes.keys()) |
| 1066 | for gi in lp: |
| 1067 | if self.expr(gi) > filtercells and self.delta(gi)[2] > filterexp: |
| 1068 | pol.append(self.connectivity_pvalue(gi, n=n)) |
| 1069 | por = benjamini_hochberg(pol) |
| 1070 | mj = 0 |
| 1071 | for gi in lp: |
| 1072 | po = self.expr(gi) |
| 1073 | m1, m2, m3 = self.delta(gi) |
| 1074 | if po > filtercells and m3 > filterexp: |
| 1075 | cul = gi + '\t' + str(po) + '\t' + str(m1) + '\t' + str(m2) + '\t' + str(m3) + '\t' + \ |
| 1076 | str(self.connectivity(gi)) + '\t' + str(pol[mj]) + '\t' + str(por[mj]) + '\t' |
| 1077 | for m in sorted(annotation.keys()): |
| 1078 | if gi in annotation[m]: |
| 1079 | cul += 'Y' + '\t' |
| 1080 | else: |
| 1081 | cul += 'N' + '\t' |
| 1082 | ggg.write(cul[:-1] + '\n') |
| 1083 | mj += 1 |
| 1084 | centr = [] |
| 1085 | disp = [] |
| 1086 | centr2 = [] |
| 1087 | disp2 = [] |
| 1088 | f = open(self.name + '.genes.tsv', 'r') |
| 1089 | for n, line in enumerate(f): |
| 1090 | if n > 0: |
| 1091 | sp = line[:-1].split('\t') |
| 1092 | if float(sp[7]) <= 0.05: |
| 1093 | centr.append(float(sp[1])) |
| 1094 | disp.append(float(sp[5])) |
| 1095 | else: |
| 1096 | centr2.append(float(sp[1])) |
| 1097 | disp2.append(float(sp[5])) |
| 1098 | f.close() |
| 1099 | pylab.scatter(centr2, disp2, alpha=0.2, s=9, c='b') |
| 1100 | pylab.scatter(centr, disp, alpha=0.3, s=9, c='r') |
| 1101 | pylab.xlabel('Cells') |
| 1102 | pylab.ylabel('Connectivity') |
| 1103 | pylab.yscale('log') |
| 1104 | pylab.ylim(0.01, 1) |
| 1105 | pylab.xlim(0, max(centr+centr2)) |
| 1106 | pylab.show() |
| 1107 |
nothing calls this directly
no test coverage detected