Mercurial > repos > fcaramia > contra

comparison Contra/scripts/cn_apply_threshold.py @ 0:7564f3b1e675

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author fcaramia
date Thu, 13 Sep 2012 02:31:43 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:7564f3b1e675
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()