Mercurial > repos > mvdbeek > dapars

comparison dapars.py @ 2:47b0044bc7f8 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/mvdbeek/dapars commit a02ced0b40d8946aa78a58321bd8f9771148391e-dirty
author mvdbeek
date Wed, 28 Oct 2015 05:56:52 -0400
parents 1b20ba32b4c5
children aca85eb3eb5b
comparison
equal deleted inserted replaced
1:1b20ba32b4c5 2:47b0044bc7f8
31 help="Start search for breakpoint n nucleotides downstream of UTR start") 31 help="Start search for breakpoint n nucleotides downstream of UTR start")
32 parser.add_argument("-ct", "--coverage_threshold", required=False, type=float, default=20, 32 parser.add_argument("-ct", "--coverage_threshold", required=False, type=float, default=20,
33 help="minimum coverage in each aligment to be considered for determining breakpoints") 33 help="minimum coverage in each aligment to be considered for determining breakpoints")
34 parser.add_argument("-b", "--breakpoint_bed", required=False, type=argparse.FileType('w'), 34 parser.add_argument("-b", "--breakpoint_bed", required=False, type=argparse.FileType('w'),
35 help="Write bedfile with coordinates of breakpoint positions to supplied path.") 35 help="Write bedfile with coordinates of breakpoint positions to supplied path.")
36 parser.add_argument("-v", "--version", action='version', version='%(prog)s 0.1.4') 36 parser.add_argument("-v", "--version", action='version', version='%(prog)s 0.1.5')
37 return parser.parse_args() 37 return parser.parse_args()
38 38
39 39
40 class UtrFinder(): 40 class UtrFinder():
41 """ 41 """
125 utr_d["new_start"] = utr_d["start"] + end_shift 125 utr_d["new_start"] = utr_d["start"] + end_shift
126 if utr_d["new_start"] + 50 < utr_d["new_end"]: 126 if utr_d["new_start"] + 50 < utr_d["new_end"]:
127 utr_dict[gene] = utr_d 127 utr_dict[gene] = utr_d
128 return utr_dict 128 return utr_dict
129 129
130 def get_utr_coverage(self):
131 """
132 Returns a dict:
133 { UTR : [coverage_aligment1, ...]}
134 """
135 utr_coverages = {}
136 for utr, utr_d in self.utr_dict.iteritems():
137 if utr_d["chr"] in self.available_chromosomes:
138 if utr_d["strand"] == "+":
139 is_reverse = False
140 else:
141 is_reverse = True
142 utr_coverage = []
143 for bam in self.all_alignments:
144 bp_coverage = get_bp_coverage(bam, utr_d["chr"], utr_d["new_start"], utr_d["new_end"], is_reverse)
145 utr_coverage.append(bp_coverage)
146 utr_coverages[utr] = utr_coverage
147 return utr_coverages
148
149 def get_coverage_weights(self): 130 def get_coverage_weights(self):
150 """ 131 """
151 Return weights for normalizing coverage. 132 Return weights for normalizing coverage.
152 utr_coverage is still confusing. 133 utr_coverage is still confusing.
153 """ 134 """
226 num_control, num_treatment) 207 num_control, num_treatment)
227 if treatment_breakpoint: 208 if treatment_breakpoint:
228 breakpoint_to_result(res_treatment, utr, utr_d, treatment_breakpoint, "treatment_breakpoint", treatment_abundance, is_reverse, 209 breakpoint_to_result(res_treatment, utr, utr_d, treatment_breakpoint, "treatment_breakpoint", treatment_abundance, is_reverse,
229 num_samples, num_control, num_treatment) 210 num_samples, num_control, num_treatment)
230 return res_control, res_treatment 211 return res_control, res_treatment
231
232
233 212
234 213
235 def breakpoint_to_result(res, utr, utr_d, breakpoint, breakpoint_type, 214 def breakpoint_to_result(res, utr, utr_d, breakpoint, breakpoint_type,
236 abundances, is_reverse, num_samples, num_control, num_treatment): 215 abundances, is_reverse, num_samples, num_control, num_treatment):
237 """ 216 """