Mercurial > repos > drosofff > msp_sr_signature

comparison signature.py @ 2:2b30861d95f4 draft

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