Mercurial > repos > drosofff > msp_sr_signature
comparison signature.py @ 0:a2f293717ce3 draft
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/msp_sr_signature commit fe40dec87779c1fcfbd03330e653aa886f4a2cda
author | drosofff |
---|---|
date | Wed, 21 Oct 2015 11:35:25 -0400 |
parents | |
children | 6218b518cd16 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:a2f293717ce3 |
---|---|
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 | |
9 import sys | |
10 import subprocess | |
11 import argparse | |
12 from smRtools import * | |
13 from collections import defaultdict # test whether it is required | |
14 | |
15 | |
16 def Parser(): | |
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 | |
46 | |
47 args = Parser() | |
48 | |
49 if args.extract_index: | |
50 GenomeFormat = "bowtieIndex" | |
51 else: | |
52 GenomeFormat = "fastaSource" | |
53 | |
54 if args.inputFormat == "tabular": | |
55 Genome = HandleSmRNAwindows( | |
56 args.input, args.inputFormat, args.referenceGenome, GenomeFormat) | |
57 elif args.inputFormat == "sam": | |
58 Genome = HandleSmRNAwindows( | |
59 args.input, args.inputFormat, args.referenceGenome, GenomeFormat) | |
60 else: | |
61 Genome = HandleSmRNAwindows( | |
62 args.input, args.inputFormat, args.referenceGenome, GenomeFormat) | |
63 | |
64 # replace objDic by Genome.instanceDict or... objDic = Genome.instanceDict | |
65 objDic = Genome.instanceDict | |
66 | |
67 args.maxscope += 1 | |
68 | |
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") | |
74 | |
75 if args.graph == "global": | |
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]) | |
105 | |
106 else: | |
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) | |
116 | |
117 OUT.close() | |
118 | |
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 | |
122 process = subprocess.Popen(R_command.split()) |