|
0
|
1 #!/usr/bin/python
|
|
|
2 # script for computing overlap signatures from a bowtie output
|
|
|
3 # Christophe Antoniewski <drosofff@gmail.com>
|
|
|
4 # Usage signature.py <1:input> <2:format of input> <3:minsize query> <4:maxsize query> <5:minsize target> <6:maxsize target>
|
|
|
5 # <7:minscope> <8:maxscope> <9:output> <10:bowtie index> <11:procedure option> <12: graph (global or lattice)>
|
|
|
6 # <13: R code>
|
|
|
7 # version 2.0.0
|
|
|
8
|
|
2
|
9 import sys
|
|
|
10 import subprocess
|
|
|
11 import argparse
|
|
0
|
12 from smRtools import *
|
|
2
|
13 from collections import defaultdict # test whether it is required
|
|
|
14
|
|
0
|
15
|
|
|
16 def Parser():
|
|
2
|
17 the_parser = argparse.ArgumentParser()
|
|
|
18 the_parser.add_argument(
|
|
|
19 '--input', action="store", type=str, help="input alignment file")
|
|
|
20 the_parser.add_argument('--inputFormat', action="store", type=str, choices=[
|
|
|
21 "tabular", "bam", "sam"], help="format of alignment file (tabular/bam/sam)")
|
|
|
22 the_parser.add_argument(
|
|
|
23 '--minquery', type=int, help="Minimum readsize of query reads (nt) - must be an integer")
|
|
|
24 the_parser.add_argument(
|
|
|
25 '--maxquery', type=int, help="Maximum readsize of query reads (nt) - must be an integer")
|
|
|
26 the_parser.add_argument(
|
|
|
27 '--mintarget', type=int, help="Minimum readsize of target reads (nt) - must be an integer")
|
|
|
28 the_parser.add_argument(
|
|
|
29 '--maxtarget', type=int, help="Maximum readsize of target reads (nt) - must be an integer")
|
|
|
30 the_parser.add_argument(
|
|
|
31 '--minscope', type=int, help="Minimum overlap analyzed (nt) - must be an integer")
|
|
|
32 the_parser.add_argument(
|
|
|
33 '--maxscope', type=int, help="Maximum overlap analyzed (nt) - must be an integer")
|
|
|
34 the_parser.add_argument(
|
|
|
35 '--outputOverlapDataframe', action="store", type=str, help="Overlap dataframe")
|
|
|
36 the_parser.add_argument('--referenceGenome', action='store',
|
|
|
37 help="path to the bowtie-indexed or fasta reference")
|
|
|
38 the_parser.add_argument('--extract_index', action='store_true',
|
|
|
39 help="specify if the reference is an indexed Bowtie reference")
|
|
|
40 the_parser.add_argument('--graph', action='store', choices=[
|
|
|
41 "global", "lattice"], help="small RNA signature is computed either globally or by item (global-lattice)")
|
|
|
42 the_parser.add_argument(
|
|
|
43 '--rcode', type=str, help="R code to be passed to the python script")
|
|
|
44 args = the_parser.parse_args()
|
|
|
45 return args
|
|
0
|
46
|
|
|
47 args = Parser()
|
|
|
48
|
|
|
49 if args.extract_index:
|
|
2
|
50 GenomeFormat = "bowtieIndex"
|
|
0
|
51 else:
|
|
2
|
52 GenomeFormat = "fastaSource"
|
|
0
|
53
|
|
|
54 if args.inputFormat == "tabular":
|
|
2
|
55 Genome = HandleSmRNAwindows(
|
|
|
56 args.input, args.inputFormat, args.referenceGenome, GenomeFormat)
|
|
0
|
57 elif args.inputFormat == "sam":
|
|
2
|
58 Genome = HandleSmRNAwindows(
|
|
|
59 args.input, args.inputFormat, args.referenceGenome, GenomeFormat)
|
|
0
|
60 else:
|
|
2
|
61 Genome = HandleSmRNAwindows(
|
|
|
62 args.input, args.inputFormat, args.referenceGenome, GenomeFormat)
|
|
0
|
63
|
|
|
64 # replace objDic by Genome.instanceDict or... objDic = Genome.instanceDict
|
|
|
65 objDic = Genome.instanceDict
|
|
|
66
|
|
|
67 args.maxscope += 1
|
|
|
68
|
|
2
|
69 general_frequency_table = dict(
|
|
|
70 [(i, 0) for i in range(args.minscope, args.maxscope)])
|
|
|
71 general_percent_table = dict(
|
|
|
72 [(i, 0) for i in range(args.minscope, args.maxscope)])
|
|
|
73 OUT = open(args.outputOverlapDataframe, "w")
|
|
0
|
74
|
|
|
75 if args.graph == "global":
|
|
2
|
76 # for normalized summing of local_percent_table(s)
|
|
|
77 readcount_dic = {}
|
|
|
78 Total_read_in_objDic = 0
|
|
|
79 for item in objDic:
|
|
|
80 readcount_dic[item] = objDic[item].readcount(
|
|
|
81 args.minquery, args.maxquery)
|
|
|
82 Total_read_in_objDic += readcount_dic[item]
|
|
|
83 ######
|
|
|
84 for x in (objDic):
|
|
|
85 local_frequency_table = objDic[x].signature(
|
|
|
86 args.minquery, args.maxquery, args.mintarget, args.maxtarget, range(args.minscope, args.maxscope))
|
|
|
87 local_percent_table = objDic[x].hannon_signature(
|
|
|
88 args.minquery, args.maxquery, args.mintarget, args.maxtarget, range(args.minscope, args.maxscope))
|
|
|
89 try:
|
|
|
90 for overlap in local_frequency_table.keys():
|
|
|
91 general_frequency_table[overlap] = general_frequency_table.get(
|
|
|
92 overlap, 0) + local_frequency_table[overlap]
|
|
|
93 except:
|
|
|
94 pass
|
|
|
95 try:
|
|
|
96 for overlap in local_percent_table.keys():
|
|
|
97 general_percent_table[overlap] = general_percent_table.get(
|
|
|
98 overlap, 0) + (1. / Total_read_in_objDic * readcount_dic[x] * local_percent_table[overlap])
|
|
|
99 except:
|
|
|
100 pass
|
|
|
101 print >> OUT, "overlap\tnum of pairs\tprobability"
|
|
|
102 for classe in sorted(general_frequency_table):
|
|
|
103 print >> OUT, "%i\t%i\t%f" % (
|
|
|
104 classe, general_frequency_table[classe], general_percent_table[classe])
|
|
0
|
105
|
|
|
106 else:
|
|
2
|
107 print >> OUT, "overlap\tnum of pairs\tprobability\titem"
|
|
|
108 for x in (objDic):
|
|
|
109 local_frequency_table = objDic[x].signature(
|
|
|
110 args.minquery, args.maxquery, args.mintarget, args.maxtarget, range(args.minscope, args.maxscope))
|
|
|
111 local_percent_table = objDic[x].hannon_signature(
|
|
|
112 args.minquery, args.maxquery, args.mintarget, args.maxtarget, range(args.minscope, args.maxscope))
|
|
|
113 for classe in range(args.minscope, args.maxscope):
|
|
|
114 print >> OUT, "%i\t%i\t%f\t%s" % (
|
|
|
115 classe, local_frequency_table[classe], local_percent_table[classe], x)
|
|
0
|
116
|
|
|
117 OUT.close()
|
|
|
118
|
|
2
|
119 # Run the R script that is defined in the xml using the Rscript binary
|
|
|
120 # provided with R.
|
|
|
121 R_command = "Rscript " + args.rcode
|
|
0
|
122 process = subprocess.Popen(R_command.split())
|
