annotate piRNAsignature.py @ 7:77dd5960dc3d draft default tip

Uploaded
author drosofff
date Mon, 23 Jun 2014 05:27:03 -0400
parents 8d4ca527888b
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
1 #!/usr/bin/python
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
2 # script for computing overlap signatures from a bowtie output
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
3 # Christophe Antoniewski <drosofff@gmail.com>
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
4 # Usage piRNAsignature.py <1:input> <2:format of input> <3:minsize query> <4:maxsize query> <5:minsize target> <6:maxsize target>
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
5 # <7:minscope> <8:maxscope> <9:output> <10:bowtie index> <11:procedure option> <12: graph (global or lattice)>
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
6 # <13: R code>
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
7
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
8 import sys, subprocess
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
9 from smRtools import *
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
10 from collections import defaultdict # test whether it is required
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
11
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
12 if sys.argv[11] == "--extract_index":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
13 if sys.argv[2] == "tabular":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
14 Genome = HandleSmRNAwindows (sys.argv[1],"tabular",sys.argv[10],"bowtieIndex")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
15 elif sys.argv[2] == "sam":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
16 Genome = HandleSmRNAwindows (sys.argv[1],"sam",sys.argv[10],"bowtieIndex")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
17 else:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
18 Genome = HandleSmRNAwindows (sys.argv[1],"bam",sys.argv[10],"bowtieIndex")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
19 else:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
20 if sys.argv[2] == "tabular":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
21 Genome = HandleSmRNAwindows (sys.argv[1],"tabular",sys.argv[10],"fastaSource")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
22 elif sys.argv[2] == "sam":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
23 Genome = HandleSmRNAwindows (sys.argv[1],"sam",sys.argv[10],"fastaSource")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
24 else:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
25 Genome = HandleSmRNAwindows (sys.argv[1],"bam",sys.argv[10],"fastaSource")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
26 # this decisional tree may be simplified if sam and bam inputs are treated the same way by pysam
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
27
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
28 # replace objDic by Genome.instanceDict or... objDic = Genome.instanceDict
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
29 objDic = Genome.instanceDict
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
30
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
31 minquery = int(sys.argv[3])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
32 maxquery = int(sys.argv[4])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
33 mintarget = int(sys.argv[5])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
34 maxtarget = int(sys.argv[6])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
35 minscope = int(sys.argv[7])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
36 maxscope = int(sys.argv[8]) + 1
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
37 general_frequency_table = dict ([(i,0) for i in range(minscope,maxscope)])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
38 general_percent_table = dict ([(i,0) for i in range(minscope,maxscope)])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
39 OUT = open (sys.argv[9], "w")
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
40
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
41 if sys.argv[12] == "global":
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
42 ###### for normalized summing of local_percent_table(s)
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
43 readcount_dic = {}
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
44 Total_read_in_objDic = 0
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
45 for item in objDic:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
46 readcount_dic[item] = objDic[item].readcount(minquery, maxquery)
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
47 Total_read_in_objDic += readcount_dic[item]
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
48 ######
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
49 for x in (objDic):
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
50 local_frequency_table = objDic[x].signature( minquery, maxquery, mintarget, maxtarget, range(minscope,maxscope) )
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
51 local_percent_table = objDic[x].hannon_signature( minquery, maxquery, mintarget, maxtarget, range(minscope,maxscope) )
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
52 try:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
53 for overlap in local_frequency_table.keys():
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
54 general_frequency_table[overlap] = general_frequency_table.get(overlap, 0) + local_frequency_table[overlap]
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
55 except:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
56 pass
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
57 try:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
58 for overlap in local_percent_table.keys():
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
59 general_percent_table[overlap] = general_percent_table.get(overlap, 0) + (1./Total_read_in_objDic*readcount_dic[x]*local_percent_table[overlap])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
60 except:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
61 pass
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
62 print >> OUT, "overlap\tnum of pairs\tprobability"
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
63 for classe in sorted(general_frequency_table):
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
64 print >> OUT, "%i\t%i\t%f" % (classe, general_frequency_table[classe], general_percent_table[classe])
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
65
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
66 else:
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
67 print >> OUT, "overlap\tnum of pairs\tprobability\titem"
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
68 for x in (objDic):
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
69 local_frequency_table = objDic[x].signature( minquery, maxquery, mintarget, maxtarget, range(minscope,maxscope) )
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
70 local_percent_table = objDic[x].hannon_signature( minquery, maxquery, mintarget, maxtarget, range(minscope,maxscope) )
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
71 for classe in range(minscope,maxscope):
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
72 print >> OUT, "%i\t%i\t%f\t%s" % (classe, local_frequency_table[classe], local_percent_table[classe], x)
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
73
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
74 OUT.close()
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
75
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
76 ## Run the R script that is defined in the xml using the Rscript binary provided with R.
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
77 R_command="Rscript "+ sys.argv[13]
8d4ca527888b Uploaded
drosofff
parents:
diff changeset
78 process = subprocess.Popen(R_command.split())