|
0
|
1 # ----------------------------------------------------------------------#
|
|
|
2 # Copyright (c) 2011, Richard Lupat & Jason Li.
|
|
|
3 #
|
|
|
4 # > Source License <
|
|
|
5 # This file is part of CONTRA.
|
|
|
6 #
|
|
|
7 # CONTRA is free software: you can redistribute it and/or modify
|
|
|
8 # it under the terms of the GNU General Public License as published by
|
|
|
9 # the Free Software Foundation, either version 3 of the License, or
|
|
|
10 # (at your option) any later version.
|
|
|
11 #
|
|
|
12 # CONTRA is distributed in the hope that it will be useful,
|
|
|
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
15 # GNU General Public License for more details.
|
|
|
16 #
|
|
|
17 # You should have received a copy of the GNU General Public License
|
|
|
18 # along with CONTRA. If not, see <http://www.gnu.org/licenses/>.
|
|
|
19 #
|
|
|
20 #
|
|
|
21 #-----------------------------------------------------------------------#
|
|
|
22 # Last Updated : 12 Oct 2011 11:00AM
|
|
|
23
|
|
|
24 def applyThreshold(outputName, bufTable, threshold, maxGap):
|
|
|
25 srcFile = outputName + ".txt"
|
|
|
26 outFile = bufTable + ".LargeVariations.txt"
|
|
|
27 bedFile = bufTable + ".BED"
|
|
|
28 fFile = outputName + ".DetailsFILTERED.txt"
|
|
|
29 ts = float(threshold)
|
|
|
30
|
|
|
31 # Read and open files
|
|
|
32 srcTable = file.readlines(open(srcFile))
|
|
|
33 outTable = open(outFile, "w")
|
|
|
34 bedOut = open(bedFile, "w")
|
|
|
35 filteredTable = open(fFile, "w")
|
|
|
36
|
|
|
37
|
|
|
38 #header
|
|
|
39 outTable.write("Chr \tStartCoordinate \tEndCoordinate \tGenes \tGain.Loss \n")
|
|
|
40 filteredTable.write(srcTable[0])
|
|
|
41
|
|
|
42 prevChr = ''
|
|
|
43 prevStatus = ''
|
|
|
44 prevEnd = -1
|
|
|
45 genes = []
|
|
|
46 chrList = []
|
|
|
47
|
|
|
48 for exons in srcTable:
|
|
|
49 exon = exons.split()
|
|
|
50 try:
|
|
|
51 adjPVal = float(exon[12])
|
|
|
52 except:
|
|
|
53 continue
|
|
|
54
|
|
|
55 if adjPVal <= ts:
|
|
|
56 chr = exon[3]
|
|
|
57 gene = exon[2]
|
|
|
58 status = exon[13]
|
|
|
59 start = exon[4]
|
|
|
60 end = exon[5]
|
|
|
61
|
|
|
62 # For first row
|
|
|
63 if prevEnd == -1:
|
|
|
64 gap = 0
|
|
|
65 else:
|
|
|
66 gap = int(prevEnd) - int(start)
|
|
|
67
|
|
|
68 # Write Filtered Table
|
|
|
69 filteredTable.write(exons)
|
|
|
70
|
|
|
71 # Write Bed File
|
|
|
72 bedOut.write(chr.strip("chr") +"\t" +start +"\t"+ end+"\t"+
|
|
|
73 chr.strip("chr")+":"+start+"-"+end+":"+str(adjPVal)+"\n")
|
|
|
74
|
|
|
75 if prevChr == '' and prevStatus == '':
|
|
|
76 if chr not in chrList:
|
|
|
77 print chr
|
|
|
78 chrList.append(chr)
|
|
|
79 elif (chr == prevChr) and (status == prevStatus) and (gap < maxGap):
|
|
|
80 start = prevStart
|
|
|
81 else:
|
|
|
82 outTable.write(prevChr +"\t" +prevStart +"\t" +prevEnd + "\t")
|
|
|
83 for gsym in genes:
|
|
|
84 outTable.write(gsym + ", ")
|
|
|
85 outTable.write("\t" + prevStatus + "\n")
|
|
|
86 genes=[]
|
|
|
87
|
|
|
88 if gene not in genes:
|
|
|
89 genes.append(gene)
|
|
|
90 prevChr = chr
|
|
|
91 prevStatus = status
|
|
|
92 prevStart = start
|
|
|
93 prevEnd = end
|
|
|
94 elif len(genes) > 0:
|
|
|
95 outTable.write(prevChr +"\t" +prevStart +"\t" +prevEnd + "\t")
|
|
|
96 for gsym in genes:
|
|
|
97 outTable.write(gsym + ", " )
|
|
|
98 outTable.write("\t" + prevStatus + "\n")
|
|
|
99 prevChr = ''
|
|
|
100 prevStatus = ''
|
|
|
101 genes = []
|
|
|
102
|
|
|
103 if len(genes) > 0:
|
|
|
104 outTable.write(prevChr +"\t" +prevStart +"\t" +prevEnd + "\t")
|
|
|
105 for gsym in genes:
|
|
|
106 outTable.write(gsym + ", ")
|
|
|
107 outTable.write("\t" + prevStatus + "\n")
|
|
|
108
|
|
|
109 filteredTable.close()
|
|
|
110 bedOut.close()
|
|
|
111 outTable.close()
|
