0
+ − 1 #!/usr/bin/env python
+ − 2 """
+ − 3 Snoop thru a fasta file looking for microsatellite repeats of given periods
+ − 4 Output format: length_of_repeat left_flank_length right_flank_length repeat_motif hamming_distance read_name read_sequence read_quality (additional columns)
+ − 5
+ − 6 If --r option turned on, output format will have additional columns behind:
+ − 7 read_name read_chr pre_s pre_e tr_s tr_e suf_s suf_e tr_len tr_ref_seq
+ − 8
+ − 9 pre_s where the read start
+ − 10 pre_e the last position before microsatellite
+ − 11 tr_s where microsatellite start
+ − 12 tr_e where microsatellite end
+ − 13 suf_s first base after microsatellite
+ − 14 tr_ref_seq reference sequence corresponding to microsatellite
+ − 15
+ − 16 * output positions are 0 based
+ − 17
+ − 18 :Author: Chen Sun (cxs1031@cse.psu.edu); Bob Harris (rsharris@bx.psu.edu)
+ − 19
+ − 20 modifing log:
+ − 21
+ − 22 09/27/2013
+ − 23 replace function dense_intervals with function non_negative_intervals, which do not need to import such file.
+ − 24
+ − 25 10/18/2013
+ − 26 modify function find_repeat_element to get a quick speed, under the condition that hamming_distance = 0, which means do not allowed any mutation/indel
+ − 27
+ − 28 02/25/2014
+ − 29 add function that can deal with mapped reads
+ − 30 with additional output
+ − 31
+ − 32 02/28/2014
+ − 33 modify the 0-based end point, as in 0-base area, it is half-open [ )
+ − 34 so the 0-based site, should always be added by 1
+ − 35
+ − 36 03/05/2014
+ − 37 deal with multi-fasta
+ − 38 """
+ − 39 from sys import argv,stdin,stderr,exit
+ − 40 from string import maketrans
+ − 41 from md5 import new as md5_new
+ − 42 import re
+ − 43 #from pyfracluster import dense_intervals
+ − 44
+ − 45 def usage(s=None):
+ − 46 message = """
+ − 47 usage: microsat_snoop [fasta_file] [options]
+ − 48 <fasta_file> Name of file to read sequences from; if absent,
+ − 49 sequences are read from stdin
+ − 50 --fasta Input file is in fasta format
+ − 51 (this is the default)
+ − 52 --fastq Input file is in fastq format
+ − 53 (default is fasta unless filename is .fastq)
+ − 54 --fastq:noquals Input file is in fastq format, but discard quals
+ − 55 --sam Input file is SAM file
+ − 56 --r Indicate additional output information, if indicated,
+ − 57 --ref option is mendatory
+ − 58 --ref=<filepath> Reference file (absolute) path
+ − 59 --period=<length> (mandatory,cumulative) repeat length(s) to be
+ − 60 searched for
+ − 61 <length> is expected to be small, less than 10
+ − 62 <length> can also be a comma-separated list, or
+ − 63 a range <low>..<high>
+ − 64 --rate=<fraction> control the candidate repeat interval detector;
+ − 65 it will consider intervals with at least
+ − 66 <fraction> of matches when shifted by the period;
+ − 67 <fraction> is between 0 and 1 and can be either a
+ − 68 real number or <n>/<d>
+ − 69 (default is 6/7)
+ − 70 --minlength=<length> minimum length of intervals reported, in bp
+ − 71 (default is 20)
+ − 72 --progress=<count> how often to report the sequence we're searching
+ − 73 (default is no progress report)
+ − 74 --allowduplicates process all input sequences
+ − 75 (this is the default)
+ − 76 --noduplicates ignore any input sequence that's the same as an
+ − 77 earlier sequence
+ − 78 --nonearduplicates ignore any input sequence that has the same first
+ − 79 100 bp as an earlier sequence
+ − 80 --nonearduplicate=<length> ignore any input sequence that has the same first
+ − 81 <length> bp as an earlier sequence
+ − 82 --hamming=<count> Don't report candidate repeat intervals that have
+ − 83 more than <count> mismatches
+ − 84 (default is to do no such filtering)
+ − 85 --prefix=<length> Don't report candidate repeat intervals that
+ − 86 start within <length> of the sequence start
+ − 87 (default is to do no such filtering)
+ − 88 --suffix=<length> Don't report candidate repeat intervals that
+ − 89 end within <length> of the sequence end
+ − 90 (default is to do no such filtering)
+ − 91 --subsample=<k>/<n> Process only the <k>th sequence of every group of
+ − 92 <n> sequences; <k> ranges from 1 to <n>
+ − 93 --multipleruns Consider all candidate intervals in a sequence
+ − 94 (default is to consider only the longest)
+ − 95 --partialmotifs Consider microatelites with a partial motif
+ − 96 (default is to consider only whole motifs)
+ − 97 --splitbyvalidity Preprocess sequences, splitting at Ns; this
+ − 98 prevents candidates from including Ns
+ − 99 (default is not to split)
+ − 100 --noflankdisplay Show entire sequence as flanking regions
+ − 101 (this is the default)
+ − 102 --flankdisplay=<length> Limit length of flanking regions shown
+ − 103 --readnamesuffix=<string> Root of suffix to append to read names; e.g. 1
+ − 104 for forward, 2 for reverse; this triggers other
+ − 105 info to be included in the suffix
+ − 106 (default is "1" for fastq; no suffix for fasta)
+ − 107 --head=<number> limit the number of sequences processed
+ − 108 --markend Write a marker line upon completion
+ − 109 (default is not to write a marker)
+ − 110 --help=details Describe the process, and quit"""
+ − 111
+ − 112 if (s == None): exit (message)
+ − 113 else: exit ("%s\n%s" % (s,message))
+ − 114
+ − 115
+ − 116 detailedDescription = """In broad terms, the process works as follows:
+ − 117
+ − 118 (1) Identify intervals that are highly correlated with the interval shifted by
+ − 119 P (the repeat period). These intervals are called "runs" or "candidates".
+ − 120 The level of correlation required is controlled by rateThreshold.
+ − 121 Depending on whether we want to look for more than one microsat, we either
+ − 122 find the longest such run (simple algorithm) or many runs (more complicated
+ − 123 algorithm). The following steps are then performed on each run.
+ − 124
+ − 125 (2) Find the most likely repeat motif in the run. This is done by counting
+ − 126 all kmers (of length P) and choosing the most frequent. If that kmer is
+ − 127 itself covered by a sub-repeat we discard this run. The idea is that we
+ − 128 can ignore a 6-mer like ACGACG because we will find it when we are looking
+ − 129 for 3-mers.
+ − 130
+ − 131 (3) Once we identify the most likely repeat motif, we then modify the
+ − 132 interval, adjusting start and end to find the interval that has the fewest
+ − 133 mismatches vs. a sequence of the motif repeated (hamming distance). Only
+ − 134 whole copies of the motif are considered.
+ − 135
+ − 136 (4) At this point we have a valid microsat interval (in the eyes of the
+ − 137 program). It is subjected to some filtering stages (hamming distance or too
+ − 138 close to an end), and if it satisfies those conditions, it's reported to
+ − 139 the user."""
+ − 140
+ − 141 def main():
+ − 142 global debug
+ − 143
+ − 144 #=== parse the command line ===
+ − 145
+ − 146 inputFilename = None
+ − 147 referenceFileName = None #add by Chen Sun on 02/25
+ − 148 inputFormat = None
+ − 149 repeatPeriods = []
+ − 150 rateThreshold = 6 / 7.0
+ − 151 lengthThreshold = 20
+ − 152 reportProgress = None
+ − 153 discardDuplicates = False
+ − 154 discardNearDuplicates = False
+ − 155 nearDuplicatePrefix = 100
+ − 156 hammingThreshold = 0
+ − 157 prefixThreshold = None
+ − 158 suffixThreshold = None
+ − 159 subsampleK = None
+ − 160 subsampleN = None
+ − 161 reportMultipleRuns = False
+ − 162 allowPartialMotifs = False
+ − 163 splitByValidity = False
+ − 164 flankDisplayLimit = None
+ − 165 readNameSuffix = None
+ − 166 headLimit = None
+ − 167 markEndOfFile = False
+ − 168 additionalInfo = False
+ − 169 debug = []
+ − 170
+ − 171 for arg in argv[1:]:
+ − 172 if (arg == "--fasta"):
+ − 173 inputFormat = "fasta"
+ − 174 elif (arg == "--fastq"):
+ − 175 inputFormat = "fastq"
+ − 176 elif (arg == "--fastq:noquals"):
+ − 177 inputFormat = "fastq:noquals"
+ − 178 elif (arg == "--sam"):
+ − 179 inputFormat = "sam"
+ − 180 elif (arg == "--r"):
+ − 181 additionalInfo = True
+ − 182 elif (arg.startswith("--ref=")):
+ − 183 referenceFileName = arg.split("=",1)[1]
+ − 184 elif (arg.startswith("--period=")):
+ − 185 val = arg.split("=",1)[1]
+ − 186 for period in val.split(","):
+ − 187 if (".." in period):
+ − 188 (lowPeriod,highPeriod) = period.split("..",1)
+ − 189 lowPeriod = int(lowPeriod)
+ − 190 highPeriod = int(highPeriod)
+ − 191 for period in xrange(lowPeriod,highPeriod+1):
+ − 192 repeatPeriods += [period]
+ − 193 else:
+ − 194 repeatPeriods += [int(period)]
+ − 195 elif (arg.startswith("--rate=")):
+ − 196 val = arg.split("=",1)[1]
+ − 197 rateThreshold = float_or_fraction(val)
+ − 198 assert (0.0 < rateThreshold <= 1.0), "%s not a valid rate" % val
+ − 199 elif (arg.startswith("--minlength=")):
+ − 200 val = arg.split("=",1)[1]
+ − 201 lengthThreshold = int(val)
+ − 202 assert (lengthThreshold >= 0)
+ − 203 elif (arg.startswith("--progress=")):
+ − 204 val = arg.split("=",1)[1]
+ − 205 reportProgress = int(val)
+ − 206 elif (arg == "--allowduplicates"):
+ − 207 discardDuplicates = False
+ − 208 discardNearDuplicates = False
+ − 209 elif (arg == "--noduplicates"):
+ − 210 discardDuplicates = True
+ − 211 discardNearDuplicates = False
+ − 212 elif (arg == "--nonearduplicates"):
+ − 213 discardDuplicates = False
+ − 214 discardNearDuplicates = True
+ − 215 elif (arg.startswith("--nonearduplicate=")):
+ − 216 val = arg.split("=",1)[1]
+ − 217 discardDuplicates = False
+ − 218 discardNearDuplicates = True
+ − 219 nearDuplicatePrefix = int(val)
+ − 220 assert (nearDuplicatePrefix > 0)
+ − 221 elif (arg.startswith("--hamming=")):
+ − 222 val = arg.split("=",1)[1]
+ − 223 hammingThreshold = int(val)
+ − 224 assert (hammingThreshold >= 0)
+ − 225 elif (arg.startswith("--prefix=")):
+ − 226 val = arg.split("=",1)[1]
+ − 227 prefixThreshold = int(val)
+ − 228 assert (prefixThreshold >= 0)
+ − 229 elif (arg.startswith("--suffix=")):
+ − 230 val = arg.split("=",1)[1]
+ − 231 suffixThreshold = int(val)
+ − 232 assert (suffixThreshold >= 0)
+ − 233 elif (arg.startswith("--subsample=")):
+ − 234 val = arg.split("=",1)[1]
+ − 235 (k,n) = val.split("/",2)
+ − 236 subsampleK = int(k)
+ − 237 subsampleN = int(n)
+ − 238 assert (0 < subsampleK <= subsampleN)
+ − 239 elif (arg == "--multipleruns"):
+ − 240 reportMultipleRuns = True
+ − 241 elif (arg == "--partialmotifs"):
+ − 242 allowPartialMotifs = True
+ − 243 elif (arg == "--splitbyvalidity"):
+ − 244 splitByValidity = True
+ − 245 elif (arg == "--noflankdisplay"):
+ − 246 flankDisplayLimit = None
+ − 247 elif (arg.startswith("--flankdisplay=")):
+ − 248 val = arg.split("=",1)[1]
+ − 249 flankDisplayLimit = int(val)
+ − 250 assert (flankDisplayLimit >= 0)
+ − 251 elif (arg.startswith("--readnamesuffix")):
+ − 252 readNameSuffix = arg.split("=",1)[1]
+ − 253 elif (arg.startswith("--head=")):
+ − 254 headLimit = int_with_unit(arg.split("=",1)[1])
+ − 255 elif (arg == "--markend"):
+ − 256 markEndOfFile = True
+ − 257 elif (arg == "--help=details"):
+ − 258 exit (detailedDescription)
+ − 259 elif (arg.startswith("--debug=")):
+ − 260 debug += (arg.split("=",1)[1]).split(",")
+ − 261 elif (arg.startswith("--")):
+ − 262 usage("unrecognized option: %s" % arg)
+ − 263 elif (inputFilename == None):
+ − 264 inputFilename = arg
+ − 265 else:
+ − 266 usage("unrecognized option: %s" % arg)
+ − 267
+ − 268 #=== determine periods of interest ===
+ − 269
+ − 270 if (repeatPeriods == []):
+ − 271 usage("you gotta give me a repeat period")
+ − 272
+ − 273 if (additionalInfo == True):
+ − 274 if (referenceFileName == None):
+ − 275 usage("reference file path needed. use --ref=<reference> to indicate")
+ − 276
+ − 277 periodSeed = {}
+ − 278 for period in repeatPeriods:
+ − 279 if (period < 1): usage("period %d is not valid" % period)
+ − 280 periodSeed[period] = True
+ − 281
+ − 282 repeatPeriods = [period for period in periodSeed]
+ − 283 repeatPeriods.sort()
+ − 284
+ − 285 #=== determine input format ===
+ − 286
+ − 287 if (inputFormat == "fasta"): sequence_reader = fasta_sequences
+ − 288 elif (inputFormat == "fastq"): sequence_reader = fastq_sequences
+ − 289 elif (inputFormat == "fastq:noquals"): sequence_reader = fastq_sequences
+ − 290 elif (inputFormat == "sam"): sequence_reader = sam_sequences
+ − 291 elif (inputFilename == None): sequence_reader = fasta_sequences
+ − 292 elif (inputFilename.endswith(".fastq")): sequence_reader = fastq_sequences
+ − 293 elif (inputFilename.endswith(".fq")): sequence_reader = fastq_sequences
+ − 294 elif (inputFilename.endswith(".sam")): sequence_reader = sam_sequences
+ − 295 else: sequence_reader = fasta_sequences
+ − 296
+ − 297 if (inputFilename != None): inputF = file(inputFilename,"rt")
+ − 298 else: inputF = stdin
+ − 299
+ − 300 if (readNameSuffix == None) \
+ − 301 and (sequence_reader == fastq_sequences) \
+ − 302 and (inputFormat != "fastq:noquals"):
+ − 303 readNameSuffix = "1"
+ − 304
+ − 305 #=== process the sequences ===
+ − 306
+ − 307 refSequence = {}
+ − 308 rightName = ""
+ − 309 sequence = ""
+ − 310 if additionalInfo:
+ − 311 firstFasta = True
+ − 312 originalRefF = open(referenceFileName)
+ − 313 for line in originalRefF.readlines():
+ − 314 line = line.replace('\r','')
+ − 315 line = line.replace('\n','')
+ − 316 if line.startswith(">"):
+ − 317 if firstFasta:
+ − 318 firstFasta = False
+ − 319 else:
+ − 320 refSequence[rightName] = sequence
+ − 321 rightName = line[1:]
+ − 322 sequence = ""
+ − 323 continue
+ − 324 sequence += line
+ − 325 originalRefF.close()
+ − 326 refSequence[rightName] = sequence
+ − 327
+ − 328 sequenceSeen = {}
+ − 329
+ − 330 numSequences = 0
+ − 331 for seqInfo in sequence_reader(inputF):
+ − 332 numSequences += 1
+ − 333 if (headLimit != None) and (numSequences > headLimit):
+ − 334 print >>stderr, "limit of %d sequences reached" % headLimit
+ − 335 break
+ − 336
+ − 337 if (sequence_reader == sam_sequences):
+ − 338 #seqName,"".join(seqNucs).upper().translate(nonDnaMap), refName, pre_s, cigar
+ − 339 (name, sequence, refName, pre_s, cigar) = seqInfo
+ − 340 quals = None
+ − 341 elif (sequence_reader == fastq_sequences):
+ − 342 (name,sequence,quals) = seqInfo
+ − 343 if (inputFormat == "fastq:noquals"): quals = None
+ − 344 else:
+ − 345 (name,sequence) = seqInfo
+ − 346 quals = None
+ − 347
+ − 348 if (reportProgress != None) and (numSequences % reportProgress == 0):
+ − 349 print >>stderr, "%s %d" % (name,numSequences)
+ − 350
+ − 351 # if we're subsampling and not interested in this sequence, skip it
+ − 352
+ − 353 if (subsampleN != None):
+ − 354 if ((numSequences-1) % subsampleN != (subsampleK-1)):
+ − 355 continue
+ − 356
+ − 357 # if this sequence is shorter than the length of interest, skip it
+ − 358
+ − 359 seqLen = len(sequence)
+ − 360 if (seqLen < period) or (seqLen < lengthThreshold): continue
+ − 361
+ − 362 # if we're not interested in duplicates and this is one, skip it;
+ − 363 # note that we assume no hash collisions occur, i.e. that all hash
+ − 364 # matches are truly sequence matches
+ − 365
+ − 366 if (discardDuplicates):
+ − 367 h = hash108(sequence)
+ − 368 if (h in sequenceSeen): continue
+ − 369 sequenceSeen[h] = True
+ − 370 elif (discardNearDuplicates):
+ − 371 h = hash108(sequence[:nearDuplicatePrefix])
+ − 372 if (h in sequenceSeen): continue
+ − 373 sequenceSeen[h] = True
+ − 374
+ − 375 # split the sequence into chunks of valid nucleotides
+ − 376
+ − 377 if (splitByValidity):
+ − 378 chunks = [(start,end) for (start,end) in nucleotide_runs(sequence)]
+ − 379 else:
+ − 380 chunks = [(0,len(sequence))]
+ − 381
+ − 382 # evaluate for each period of interest
+ − 383
+ − 384 for period in repeatPeriods:
+ − 385
+ − 386 # operate on each chunk
+ − 387
+ − 388 for (chunkStart,chunkEnd) in chunks:
+ − 389 chunkLen = chunkEnd - chunkStart
+ − 390 if (chunkLen < period) or (chunkLen < lengthThreshold): continue
+ − 391
+ − 392 if ("validity" in debug) or ("correlation" in debug) or ("runs" in debug):
+ − 393 print >>stderr, ">%s_%d_%d" % (name,chunkStart,chunkEnd)
+ − 394
+ − 395 # compute correlation sequence
+ − 396
+ − 397 corr = correlation_sequence(sequence,period,chunkStart,chunkEnd)
+ − 398
+ − 399 if ("correlation" in debug) or ("runs" in debug):
+ − 400 print >>stderr, sequence[chunkStart:chunkEnd]
+ − 401 print >>stderr, corr
+ − 402
+ − 403 # find runs (candidates for being a microsat)
+ − 404
+ − 405 if (reportMultipleRuns):
+ − 406 runs = all_suitable_runs(corr,lengthThreshold-period,rateThreshold, hammingThreshold)
+ − 407 else:
+ − 408 runs = longest_suitable_run(corr,lengthThreshold,rateThreshold)
+ − 409 if (runs == []): continue
+ − 410
+ − 411
+ − 412 if ("runs" in debug):
+ − 413 for (start,end) in runs:
+ − 414 run = [" "] * seqLen
+ − 415 for ix in xrange(start-period,end):
+ − 416 run[ix] = "*"
+ − 417 print >>stderr, "".join(run)
+ − 418
+ − 419 if ("candidates" in debug):
+ − 420 for (start,end) in runs:
+ − 421 print >>stderr, "%s %d %d" % (name,start,end)
+ − 422
+ − 423 # process runs and report those that pass muster
+ − 424
+ − 425 runCount = 0
+ − 426 for (start,end) in runs:
+ − 427 runCount += 1
+ − 428
+ − 429 start = chunkStart + start - period
+ − 430 end = chunkStart + end
+ − 431
+ − 432 (kmer,d,start,end) = find_repeat_element(hammingThreshold, period,sequence,start,end,allowPartials=allowPartialMotifs)
+ − 433 if (kmer == None): continue # (no useful repeat kmer was found)
+ − 434
+ − 435 rptExtent = end - start
+ − 436 prefixLen = start
+ − 437 suffixLen = seqLen - end
+ − 438 if (rptExtent <= period): continue
+ − 439 if (hammingThreshold != None) and (d > hammingThreshold): continue
+ − 440 if (prefixThreshold != None) and (prefixLen < prefixThreshold): continue
+ − 441 if (suffixThreshold != None) and (suffixLen < suffixThreshold): continue
+ − 442
+ − 443 if (flankDisplayLimit == None):
+ − 444 seq = sequence[:start] \
+ − 445 + sequence[start:end].lower() \
+ − 446 + sequence[end:]
+ − 447 else:
+ − 448 seq = sequence[max(chunkStart,start-flankDisplayLimit):start] \
+ − 449 + sequence[start:end].lower() \
+ − 450 + sequence[end:min(chunkEnd,end+flankDisplayLimit)]
+ − 451 reportName = name
+ − 452 if (readNameSuffix != None):
+ − 453 reportName += "_"+readNameSuffix+"_per"+str(period)+"_"+str(runCount)
+ − 454 if (quals == None or quals == "." or quals == "\t."): quals = "\t."
+ − 455 else: quals = "\t" + quals
+ − 456 if not additionalInfo:
+ − 457 print "%d\t%d\t%d\t%s\t%d\t%s\t%s%s" \
+ − 458 % (rptExtent,prefixLen,suffixLen,kmer,d,reportName,seq,quals)
+ − 459 else:
+ − 460 #pre_e = pre_s + prefixLen - 1
+ − 461 refPoint = pre_s
+ − 462 donorPoint = 0
+ − 463
+ − 464 donorBeforeStart = prefixLen - 1 #pre_e
+ − 465 donorMicroStart = prefixLen #tr_s
+ − 466 donorMicroEnd = donorMicroStart + rptExtent - 1 #tr_e
+ − 467 donorAfterMicro = donorMicroEnd + 1 #suf_s
+ − 468 donorEnd = len(seq) - 1 #suf_e
+ − 469
+ − 470 set_pre_e = False
+ − 471 set_tr_s = False
+ − 472 set_tr_e = False
+ − 473 set_suf_s = False
+ − 474 set_suf_e = False
+ − 475
+ − 476 pre_e = 0
+ − 477 tr_s = 0
+ − 478 tr_e = 0
+ − 479 suf_s = 0
+ − 480 suf_e = 0
+ − 481
+ − 482 matchList = re.findall('(\d+)([IDM])', cigar)
+ − 483 unCognitiveCigar = False
+ − 484 for matchN, matchType in matchList:
+ − 485 matchNum = int(matchN)
+ − 486 if matchType == "M":
+ − 487 donorPoint = donorPoint + matchNum
+ − 488 refPoint = refPoint + matchNum
+ − 489 elif matchType == "D":
+ − 490 refPoint = refPoint + matchNum
+ − 491 continue
+ − 492 elif matchType == "I":
+ − 493 donorPoint = donorPoint + matchNum
+ − 494 else:
+ − 495 unCognitiveCigar = True
+ − 496 break
+ − 497
+ − 498 if not set_pre_e:
+ − 499 if donorPoint >= donorBeforeStart:
+ − 500 pre_e = refPoint - (donorPoint - donorBeforeStart)
+ − 501 set_pre_e = True
+ − 502 else:
+ − 503 continue
+ − 504
+ − 505 if not set_tr_s:
+ − 506 if donorPoint >= donorMicroStart:
+ − 507 tr_s = refPoint - (donorPoint - donorMicroStart)
+ − 508 set_tr_s = True
+ − 509 else:
+ − 510 continue
+ − 511
+ − 512 if not set_tr_e:
+ − 513 if donorPoint >= donorMicroEnd:
+ − 514 tr_e = refPoint - (donorPoint - donorMicroEnd)
+ − 515 set_tr_e = True
+ − 516 else:
+ − 517 continue
+ − 518
+ − 519 if not set_suf_s:
+ − 520 if donorPoint >= donorAfterMicro:
+ − 521 suf_s = refPoint - (donorPoint - donorAfterMicro)
+ − 522 set_suf_s = True
+ − 523 else:
+ − 524 continue
+ − 525
+ − 526 if not set_suf_e:
+ − 527 if donorPoint >= donorEnd:
+ − 528 suf_e = refPoint - (donorPoint - donorEnd)
+ − 529 set_suf_e = True
+ − 530 else:
+ − 531 continue
+ − 532
+ − 533 if unCognitiveCigar:
+ − 534 break
+ − 535 tr_len = tr_e - tr_s + 1
+ − 536
+ − 537 if refName not in refSequence:
+ − 538 tr_ref_seq = "."
+ − 539 else:
+ − 540 if refSequence[refName] == "":
+ − 541 tr_ref_seq = "."
+ − 542 elif len(refSequence[refName]) <= tr_e:
+ − 543 tr_ref_seq = "."
+ − 544 else:
+ − 545 tr_ref_seq = refSequence[refName][tr_s:tr_e+1]
+ − 546
+ − 547 pre_e += 1
+ − 548 tr_e += 1
+ − 549 suf_e += 1
+ − 550 print "%d\t%d\t%d\t%s\t%d\t%s\t%s%s\t%s\t%s\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s" \
+ − 551 % (rptExtent,prefixLen,suffixLen,kmer,d,reportName,seq,quals,reportName,refName,pre_s,pre_e,tr_s,tr_e,suf_s,suf_e,tr_len,tr_ref_seq)
+ − 552
+ − 553 if (markEndOfFile):
+ − 554 print "# microsat_snoop end-of-file"
+ − 555
+ − 556 if (inputF != stdin):
+ − 557 inputF.close()
+ − 558
+ − 559 # non_negative_intervals
+ − 560 # find intervals with exactly + and no -
+ − 561 # from string like this : +++++++++---+++++++++
+ − 562 def non_negative_intervals(seq, minLength=None):
+ − 563
+ − 564 start = -1
+ − 565 end = -1
+ − 566 firstPlus = 1
+ − 567 #print seq
+ − 568 for ix in range(len(seq)): # for every char in seq
+ − 569 ch = seq[ix]
+ − 570 if(ch == "+"):
+ − 571 if(firstPlus):
+ − 572 firstPlus = 0
+ − 573 start = ix
+ − 574 else:
+ − 575 continue
+ − 576 elif(ch == "-"):
+ − 577 if(start >= 0):
+ − 578 end = ix-1
+ − 579 if((end - start + 1) >= minLength):
+ − 580 yield (start,end+1)
+ − 581 start = -1
+ − 582 firstPlus = 1
+ − 583 if(start > 0):
+ − 584 if((ix - start + 1) >= minLength):
+ − 585 yield (start, ix+1)
+ − 586
+ − 587
+ − 588 ###################################################################
+ − 589 # modified by Chen Sun on 7/11/2014
+ − 590 # We do not want other modules, so parse these functions inside
+ − 591 #
+ − 592 ###################################################################
+ − 593
+ − 594 # parse a string of the form {positives}/{positives_and_neutrals}
+ − 595
+ − 596 def parse_spec(s):
+ − 597 if ("/" not in s): raise ValueError
+ − 598 (n,d) = s.split("/",1)
+ − 599 if (not n.startswith("{")) or (not n.endswith("}")): raise ValueError
+ − 600 if (not d.startswith("{")) or (not d.endswith("}")): raise ValueError
+ − 601
+ − 602 positives = n[1:-1]
+ − 603 d = d[1:-1]
+ − 604
+ − 605 for ch in positives:
+ − 606 if (ch not in d): raise ValueError
+ − 607
+ − 608 neutrals = [ch for ch in d if (ch not in positives)]
+ − 609 return (positives,neutrals)
+ − 610
+ − 611
+ − 612 # convert a string to a number, allowing fractions
+ − 613
+ − 614 def float_or_fraction(s):
+ − 615 if ("/" in s):
+ − 616 (numer,denom) = s.split("/",1)
+ − 617 return float(numer)/float(denom)
+ − 618 else:
+ − 619 return float(s)
+ − 620
+ − 621
+ − 622 # dense_intervals--
+ − 623 # Find all non-overlapping runs with a good enough rate (of positives), and
+ − 624 # which meet our length threshold.
+ − 625 #
+ − 626 # The algorithm used is adapted from Zhang, Berman, Miller, "Post-processing
+ − 627 # long pairwise alignments", Bioinformatics Vol. 15 no. 12 1999.
+ − 628 #
+ − 629 # $$$ we use the denominator as the threshold, but we really should use the
+ − 630 # $$$ .. numerator, comparing it to minLength*rate
+ − 631
+ − 632 def dense_intervals(seq,rate,positives,neutrals,blockers="",minLength=None):
+ − 633
+ − 634 if (blockers == None):
+ − 635 blockers = "".join([chr(n) for n in range(1,256)
+ − 636 if (chr(n) not in positives)
+ − 637 and (chr(n) not in neutrals)])
+ − 638
+ − 639 stackLeft = [None] # stack with each entry containing five
+ − 640 stackRight = [None] # .. elements; note that entry zero is not
+ − 641 stackLeftScore = [None] # .. used
+ − 642 stackRightScore = [None]
+ − 643 stackLower = [None]
+ − 644 top = 0
+ − 645 score = 0
+ − 646
+ − 647 for ix in range(len(seq)):
+ − 648 ch = seq[ix]
+ − 649 if (ch in blockers):
+ − 650 # emit intervals
+ − 651
+ − 652 for sp in range(1,top+1):
+ − 653 left = stackLeft [sp] + 1
+ − 654 right = stackRight[sp]
+ − 655
+ − 656 while (left < right) and (seq[left] not in positives): left += 1
+ − 657 while (right > left) and (seq[right] not in positives): right -= 1
+ − 658
+ − 659 right += 1
+ − 660 if (minLength == None) or (right - left >= minLength):
+ − 661 yield (left,right)
+ − 662
+ − 663 #empty stack
+ − 664
+ − 665 stackLeft = [None]
+ − 666 stackRight = [None]
+ − 667 stackLeftScore = [None]
+ − 668 stackRightScore = [None]
+ − 669 stackLower = [None]
+ − 670 top = 0
+ − 671 score = 0
+ − 672 continue
+ − 673
+ − 674 if (ch in positives): weight = 1-rate
+ − 675 elif (ch in neutrals): weight = -rate
+ − 676 else: raise ValueError
+ − 677
+ − 678 score += weight
+ − 679 #if ("algorithm" in debug):
+ − 680 # print >>sys.stderr, "%3d: %c %5.2f" % (ix, ch, score),
+ − 681
+ − 682 if (weight < 0):
+ − 683 #if ("algorithm" in debug):
+ − 684 # print >>sys.stderr
+ − 685 continue
+ − 686
+ − 687 if (top > 0) and (stackRight[top] == ix-1):
+ − 688 # add this site to the interval on top of the stack
+ − 689
+ − 690 stackRight [top] = ix
+ − 691 stackRightScore[top] = score
+ − 692
+ − 693 #if ("algorithm" in debug):
+ − 694 # print >>sys.stderr, \
+ − 695 # " extending [%d] %d-%d %4.1f %4.1f" \
+ − 696 # % (top,
+ − 697 # stackLeft [top], stackRight [top],
+ − 698 # stackLeftScore[top], stackRightScore[top]),
+ − 699
+ − 700 else:
+ − 701 # create a one site interval
+ − 702
+ − 703 top += 1
+ − 704 if (top >= len(stackLeft)):
+ − 705 stackLeft += [None]
+ − 706 stackRight += [None]
+ − 707 stackLeftScore += [None]
+ − 708 stackRightScore += [None]
+ − 709 stackLower += [None]
+ − 710
+ − 711 stackLeft [top] = ix - 1
+ − 712 stackLeftScore [top] = score - weight
+ − 713 stackRight [top] = ix
+ − 714 stackRightScore[top] = score
+ − 715 stackLower [top] = top - 1
+ − 716
+ − 717 while (stackLower[top] > 0) \
+ − 718 and (stackLeftScore[stackLower[top]] > stackLeftScore[top]):
+ − 719 stackLower[top] = stackLower[stackLower[top]]
+ − 720
+ − 721 #if ("algorithm" in debug):
+ − 722 # print >>sys.stderr, \
+ − 723 # " creating [%d] %d-%d %4.1f %4.1f -> %d" \
+ − 724 # % (top,
+ − 725 # stackLeft [top], stackRight [top],
+ − 726 # stackLeftScore[top], stackRightScore[top],
+ − 727 # stackLower [top]),
+ − 728
+ − 729 # merge intervals; if there is a previous interval with a no-higher
+ − 730 # left score and no-higher right score, merge this interval (and all
+ − 731 # intervening ones) into that one
+ − 732
+ − 733 while (top > 1) \
+ − 734 and (stackLower[top] > 0) \
+ − 735 and (stackRightScore[stackLower[top]] <= stackRightScore[top]):
+ − 736 stackRight [stackLower[top]] = stackRight [top]
+ − 737 stackRightScore[stackLower[top]] = stackRightScore[top]
+ − 738 top = stackLower[top]
+ − 739
+ − 740 #if ("algorithm" in debug):
+ − 741 # print >>sys.stderr, \
+ − 742 # "\n%*s merging [%d] %d-%d %4.1f %4.1f" \
+ − 743 # % (13, "", top,
+ − 744 # stackLeft[top], stackRight [top],
+ − 745 # stackLeftScore[top], stackRightScore[top]),
+ − 746
+ − 747 #if ("algorithm" in debug):
+ − 748 # print >>sys.stderr
+ − 749
+ − 750 # emit intervals
+ − 751
+ − 752 for sp in range(1,top+1):
+ − 753 left = stackLeft [sp] + 1
+ − 754 right = stackRight[sp]
+ − 755
+ − 756 while (left < right) and (seq[left] not in positives): left += 1
+ − 757 while (right > left) and (seq[right] not in positives): right -= 1
+ − 758
+ − 759 right += 1
+ − 760 if (minLength == None) or (right - left >= minLength):
+ − 761 yield (left,right)
+ − 762
+ − 763
+ − 764 ###################################################################
+ − 765 # modified by Chen Sun on 7/11/2014
+ − 766 #
+ − 767 ###################################################################
+ − 768
+ − 769 # correlation_sequence--
+ − 770 # Compute the correlation sequence for a given period. This is a sequence
+ − 771 # of + and - indicating whether the base at a given position matches the one
+ − 772 # P positions earlier (where P is the period). The first P positions are
+ − 773 # blank. Positions with single character runs longer than the period are
+ − 774 # considered as non-matches, unless the period is 1.
+ − 775
+ − 776 def correlation_sequence(sequence,period,start=None,end=None):
+ − 777 if (start == None): start = 0
+ − 778 if (end == None): end = len(sequence)
+ − 779
+ − 780 prevCh = sequence[start]
+ − 781 run = 1
+ − 782 for ix in xrange(start+1,start+period):
+ − 783 ch = sequence[ix]
+ − 784 if (ch != prevCh): run = 1
+ − 785 else: run += 1
+ − 786 prevCh = ch
+ − 787
+ − 788 corr = [" "] * period
+ − 789 for ix in xrange(start+period,end):
+ − 790 rptCh = sequence[ix-period]
+ − 791 ch = sequence[ix]
+ − 792 if (ch != prevCh): run = 1
+ − 793 else: run += 1
+ − 794 if (ch in "ACGT") \
+ − 795 and (ch == rptCh) \
+ − 796 and ((period == 1) or (run < period)):
+ − 797 corr += ["+"]
+ − 798 else:
+ − 799 corr += ["-"]
+ − 800 prevCh = ch
+ − 801
+ − 802 return "".join(corr)
+ − 803
+ − 804
+ − 805 # longest_suitable_run--
+ − 806 # Find longest run with a good enough rate (of positives).
+ − 807 #
+ − 808 # We score a "+" as 1-r and anything else as -r. This is based on the fol-
+ − 809 # lowing derivation (p is the number of "+"s, n is the number of non-"+"s):
+ − 810 # p/(p+n) >= r
+ − 811 # ==> p >= rp + rn
+ − 812 # ==> (1-r)p - rn >= 0
+ − 813 #
+ − 814 # We adapt an algorithm from "Programming Pearls", pg. 81 (2000 printing).
+ − 815 #
+ − 816 # $$$ we use the denominator as the threshold, but we really should use the
+ − 817 # $$$ .. numerator, comparing it to minLength*rate
+ − 818 #
+ − 819 # $$$ this needs to account for $$$ this situation:
+ − 820 # $$$ sequence: ACGACGACGACGTTATTATTATTA
+ − 821 # $$$ matches: +++++++++---+++++++++
+ − 822 # $$$ this is currently considered to be one interval (if rate <= 6/7), but it
+ − 823 # $$$ ought to be two; we can't just post-process, though, because some other
+ − 824 # $$$ interval might be longer than the longest half of this; maybe what we
+ − 825 # $$$ need to do is consider matches at distances -P and -2P, or if we match
+ − 826 # $$$ -P but that itself was a mismatch, we should carry the mismatch forward
+ − 827
+ − 828 def longest_suitable_run(seq,minLength,rate):
+ − 829 maxEndingHere = 0
+ − 830 maxSoFar = 0
+ − 831 start = None
+ − 832
+ − 833 for ix in xrange(len(seq)):
+ − 834 if (seq[ix] == "+"): s = 1-rate
+ − 835 else: s = -rate
+ − 836
+ − 837 if (maxEndingHere+s < 0):
+ − 838 maxEndingHere = 0
+ − 839 block = ix
+ − 840 else:
+ − 841 maxEndingHere += s
+ − 842 if (maxEndingHere >= maxSoFar):
+ − 843 maxSoFar = maxEndingHere
+ − 844 start = block + 1
+ − 845 end = ix + 1
+ − 846
+ − 847 if (start == None) or (end - start < minLength):
+ − 848 return []
+ − 849 else:
+ − 850 return [(start,end)]
+ − 851
+ − 852
+ − 853 # all_suitable_runs--
+ − 854 # Find all non-overlapping runs with a good enough rate (of positives), and
+ − 855 # which meet our length threshold.
+ − 856 # $$$ this needs to post-process the intervals, splitting them to account for
+ − 857 # $$$ this situation:
+ − 858 # $$$ sequence: ACGACGACGACGTTATTATTATTA
+ − 859 # $$$ matches: +++++++++---+++++++++
+ − 860 # $$$ this is currently reported as one interval (if rate <= 6/7), but it
+ − 861 # $$$ ought to be two
+ − 862
+ − 863 def all_suitable_runs(seq,minCorrLength,rate, hammingThreshold):
+ − 864
+ − 865 ################################################################
+ − 866 # modified by Chen Sun on 07/11/2014
+ − 867 #
+ − 868 ################################################################
+ − 869
+ − 870 if hammingThreshold > 0:
+ − 871 return [(start,end) for (start,end) in dense_intervals(seq,rate,"+","-",blockers=None,minLength=minCorrLength)]
+ − 872 elif hammingThreshold == 0:
+ − 873 return [(start,end) for (start,end) in non_negative_intervals(seq, minLength=minCorrLength)]
+ − 874
+ − 875
+ − 876 # find_repeat_element--
+ − 877 # Find the most plausible repeat element for a run, and nudge the ends of
+ − 878 # the run if needed. Note that we will not consider kmers that represent
+ − 879 # shorter repeats. For example, we won't report ACTACT as a 6-mer since we
+ − 880 # consider this to have a shorter period than 6.
+ − 881
+ − 882 def find_repeat_element(hammingThreshold, period,seq,start,end,allowPartials=False):
+ − 883
+ − 884 if hammingThreshold > 0:
+ − 885 (kmer,bestD,bestStart,bestEnd) = find_hamming_repeat_element(period,seq,start,end,allowPartials)
+ − 886 return (kmer,bestD,bestStart,bestEnd)
+ − 887 # count the number of occurences of each k-mer; note that we can't
+ − 888 # reject kmers containing smaller repeats yet, since for a sequence like
+ − 889 # ACACACACACAAACACACACACACACACAC we must first discover ACACAC as the best
+ − 890 # 6-mer, and THEN reject it; if we reject ACACAC while counting, we'd end
+ − 891 # up reporting something like ACACAA as the best motif
+ − 892
+ − 893 if ("element" in debug):
+ − 894 print >>stderr, "find_repeat_element(%d,%d,%d)" % (period,start,end)
+ − 895
+ − 896 if ("partial" in debug):
+ − 897 print period, seq, start, end, allowPartials;
+ − 898 print seq[start:end]
+ − 899
+ − 900 kmerToCount = {}
+ − 901 kmerToFirst = {}
+ − 902 for ix in xrange(start,end-(period-1)):
+ − 903 kmer = seq[ix:ix+period]
+ − 904 if ("N" in kmer): continue
+ − 905 if (kmer not in kmerToCount):
+ − 906 kmerToCount[kmer] = 1
+ − 907 kmerToFirst[kmer] = ix
+ − 908 else:
+ − 909 kmerToCount[kmer] += 1
+ − 910 #if ("element" in debug):
+ − 911 # print >>stderr, " %d: %s" % (ix,kmer)
+ − 912
+ − 913 # choose the best k-mer; this is simply the most frequently occurring one,
+ − 914 # with ties broken by whichever one came first
+ − 915
+ − 916 kmers = [(-kmerToCount[kmer],kmerToFirst[kmer],kmer) for kmer in kmerToCount]
+ − 917 if (kmers == []): return (None,None,start,end)
+ − 918 kmers.sort()
+ − 919
+ − 920 if ("element" in debug):
+ − 921 for (count,first,kmer) in kmers:
+ − 922 print >>stderr, " %s: %d" % (kmer,-count)
+ − 923
+ − 924 (count,first,kmer) = kmers[0]
+ − 925 if (contains_repeat(kmer)): return (None,None,start,end)
+ − 926
+ − 927 # determine the hamming distance between the run and a simple repeat, for
+ − 928 # each "plausible" start and end; we compute the distance for each such
+ − 929 # interval, and choose the one with the lowest hamming distance; ties are
+ − 930 # broken in a deterministic-but-unspecified manner
+ − 931
+ − 932 bestD = bestStart = bestEnd = None
+ − 933 ###################################################################################
+ − 934 # modified by Chen Sun(cxs1031@cse.psu.edu) on 10/18/2013
+ − 935 # since we do not allow hamming_distance > 0, which means we do not allow mutation,
+ − 936 # we do not need this section to produce bestStart and End
+ − 937 ###################################################################################
+ − 938
+ − 939 #for (s,e) in plausible_intervals(start,end,period,len(seq),allowPartials=allowPartials):
+ − 940 # d = hamming_distance(seq,s,e,kmer)
+ − 941 # if (d == None): continue
+ − 942 # if (bestD == None) or (d <= bestD):
+ − 943 # (bestD,bestStart,bestEnd) = (d,s,e)
+ − 944
+ − 945
+ − 946
+ − 947 bestStart = start
+ − 948
+ − 949 if(allowPartials):
+ − 950 bestEnd = end
+ − 951 elif(not allowPartials):
+ − 952 bestEnd = start
+ − 953 pattern = seq[start:start+period]
+ − 954 if ("partial" in debug):
+ − 955 print "kmer:", kmer
+ − 956 if(pattern != kmer):
+ − 957 print "pattern:", pattern
+ − 958
+ − 959 while(bestEnd <= end-period):
+ − 960 bestEnd += period
+ − 961
+ − 962 # bestD will always be 0, as we do not allow mutation
+ − 963 bestD = 0
+ − 964
+ − 965 if ("partial" in debug):
+ − 966 print bestD, bestStart, bestEnd
+ − 967
+ − 968 ###################################################################################
+ − 969 # modified by Chen Sun(cxs1031@cse.psu.edu) on 10/10
+ − 970 #
+ − 971 ###################################################################################
+ − 972 return (kmer,bestD,bestStart,bestEnd)
+ − 973
+ − 974
+ − 975 def find_hamming_repeat_element(period,seq,start,end,allowPartials=False):
+ − 976
+ − 977 # count the number of occurences of each k-mer; note that we can't
+ − 978 # reject kmers containing smaller repeats yet, since for a sequence like
+ − 979 # ACACACACACAAACACACACACACACACAC we must first discover ACACAC as the best
+ − 980 # 6-mer, and THEN reject it; if we reject ACACAC while counting, we'd end
+ − 981 # up reporting something like ACACAA as the best motif
+ − 982
+ − 983 if ("element" in debug):
+ − 984 print >>stderr, "find_repeat_element(%d,%d,%d)" % (period,start,end)
+ − 985
+ − 986 kmerToCount = {}
+ − 987 kmerToFirst = {}
+ − 988 for ix in xrange(start,end-(period-1)):
+ − 989 kmer = seq[ix:ix+period]
+ − 990 if ("N" in kmer): continue
+ − 991 if (kmer not in kmerToCount):
+ − 992 kmerToCount[kmer] = 1
+ − 993 kmerToFirst[kmer] = ix
+ − 994 else:
+ − 995 kmerToCount[kmer] += 1
+ − 996 #if ("element" in debug):
+ − 997 # print >>stderr, " %d: %s" % (ix,kmer)
+ − 998
+ − 999 # choose the best k-mer; this is simply the most frequently occurring one,
+ − 1000 # with ties broken by whichever one came first
+ − 1001
+ − 1002 kmers = [(-kmerToCount[kmer],kmerToFirst[kmer],kmer) for kmer in kmerToCount]
+ − 1003 if (kmers == []): return (None,None,start,end)
+ − 1004 kmers.sort()
+ − 1005
+ − 1006 if ("element" in debug):
+ − 1007 for (count,first,kmer) in kmers:
+ − 1008 print >>stderr, " %s: %d" % (kmer,-count)
+ − 1009
+ − 1010 (count,first,kmer) = kmers[0]
+ − 1011 if (contains_repeat(kmer)): return (None,None,start,end)
+ − 1012
+ − 1013 # determine the hamming distance between the run and a simple repeat, for
+ − 1014 # each "plausible" start and end; we compute the distance for each such
+ − 1015 # interval, and choose the one with the lowest hamming distance; ties are
+ − 1016 # broken in a deterministic-but-unspecified manner
+ − 1017
+ − 1018 bestD = bestStart = bestEnd = None
+ − 1019
+ − 1020 for (s,e) in plausible_intervals(start,end,period,len(seq),allowPartials=allowPartials):
+ − 1021 d = hamming_distance(seq,s,e,kmer)
+ − 1022 if (d == None): continue
+ − 1023 if (bestD == None) or (d <= bestD):
+ − 1024 (bestD,bestStart,bestEnd) = (d,s,e)
+ − 1025
+ − 1026 return (kmer,bestD,bestStart,bestEnd)
+ − 1027
+ − 1028 # plausible_intervals--
+ − 1029 # Yield all plausible intervals intersecting with a run. We generate all
+ − 1030 # starts within P bp of the run's start. For each of these, we either (a) try
+ − 1031 # all ends within P bp of run's end, or (b) trim the new interval to a whole
+ − 1032 # multiple of the period, and report this short interval and the longer
+ − 1033 # interval with one more period appended. Case (a) allows partial motifs,
+ − 1034 # while case (b) only allows whole motifs.
+ − 1035
+ − 1036 def plausible_intervals(start,end,period,seqLen,allowPartials=False):
+ − 1037
+ − 1038 # generate intervals that allow a partial copy of the motif
+ − 1039
+ − 1040 if (allowPartials):
+ − 1041 for candStart in xrange(start-(period-1),start+period):
+ − 1042 if (candStart < 0): continue
+ − 1043 for candEnd in xrange(end-(period-1),end+period):
+ − 1044 if (candEnd > seqLen): continue
+ − 1045 if (candEnd <= candStart+period): continue
+ − 1046 yield (candStart,candEnd)
+ − 1047
+ − 1048 # -OR- generate intervals that allow only whole copies of the motif
+ − 1049
+ − 1050 else:
+ − 1051 for candStart in xrange(start-(period-1),start+period):
+ − 1052 if (candStart < 0): continue
+ − 1053 candEnd = candStart + ((end-candStart)/period)*period
+ − 1054 yield (candStart,candEnd)
+ − 1055 candEnd += period
+ − 1056 if (candEnd <= seqLen): yield (candStart,candEnd)
+ − 1057
+ − 1058
+ − 1059 # hamming_distance--
+ − 1060 # Determine the hamming distance between the run and a simple repeat.
+ − 1061 # $$$ improve this by allowing gaps, and stopping when we reach a threshold
+ − 1062
+ − 1063 kmerToDiffs = {} # (this is used for memo-ization)
+ − 1064
+ − 1065 def hamming_distance(seq,start,end,kmer):
+ − 1066 period = len(kmer)
+ − 1067 if (end < start + period): return None
+ − 1068
+ − 1069 wholeEnd = start + ((end-start)/period)*period
+ − 1070
+ − 1071 if (kmer not in kmerToDiffs):
+ − 1072 kmerToDiffs[kmer] = { kmer:0 }
+ − 1073
+ − 1074 d = 0
+ − 1075 for ix in xrange(start,wholeEnd,period):
+ − 1076 qmer = seq[ix:ix+period] # same size as the kmer motif
+ − 1077 if (qmer in kmerToDiffs[kmer]):
+ − 1078 d += kmerToDiffs[kmer][qmer]
+ − 1079 continue
+ − 1080 diffs = 0
+ − 1081 for iy in xrange(0,period):
+ − 1082 if (qmer[iy] != kmer[iy]): diffs += 1
+ − 1083 kmerToDiffs[kmer][qmer] = diffs
+ − 1084 d += diffs
+ − 1085
+ − 1086 if (end > wholeEnd):
+ − 1087 qmer = seq[wholeEnd:end] # shorter than the kmer motif
+ − 1088 if (qmer in kmerToDiffs[kmer]):
+ − 1089 d += kmerToDiffs[kmer][qmer]
+ − 1090 else:
+ − 1091 diffs = 0
+ − 1092 for iy in xrange(0,len(qmer)):
+ − 1093 if (qmer[iy] != kmer[iy]): diffs += 1
+ − 1094 kmerToDiffs[kmer][qmer] = diffs
+ − 1095 d += diffs
+ − 1096
+ − 1097 return d
+ − 1098
+ − 1099
+ − 1100 # fasta_sequences--
+ − 1101 # Read the fasta sequences from a file. Note that we convert to upper case,
+ − 1102 # and convert any letter other than ACGT to N.
+ − 1103
+ − 1104 nonDnaMap = maketrans("BDEFHIJKLMOPQRSUVWXYZ","NNNNNNNNNNNNNNNNNNNNN")
+ − 1105
+ − 1106 def fasta_sequences(f):
+ − 1107 seqName = None
+ − 1108 seqNucs = None
+ − 1109
+ − 1110 for line in f:
+ − 1111 line = line.strip()
+ − 1112 if (line.startswith(">")):
+ − 1113 if (seqName != None):
+ − 1114 yield (seqName,"".join(seqNucs))
+ − 1115 seqName = sequence_name(line)
+ − 1116 seqNucs = []
+ − 1117 elif (seqName == None):
+ − 1118 assert (False), "first sequence has no header"
+ − 1119 else:
+ − 1120 seqNucs += [line]
+ − 1121
+ − 1122 if (seqName != None):
+ − 1123 yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap))
+ − 1124
+ − 1125
+ − 1126 # fastq_sequences--
+ − 1127 # Read the fastq sequences from a file. Note that we convert to upper case,
+ − 1128 # and convert any letter other than ACGT to N.
+ − 1129
+ − 1130 def fastq_sequences(f):
+ − 1131 lineNum = 0
+ − 1132 for line in f:
+ − 1133 lineNum += 1
+ − 1134 line = line.strip()
+ − 1135
+ − 1136 if (lineNum % 4 == 1):
+ − 1137 assert (line.startswith("@")), \
+ − 1138 "bad read name at line %d" % lineNum
+ − 1139 seqName = line[1:]
+ − 1140 continue
+ − 1141
+ − 1142 if (lineNum % 4 == 2):
+ − 1143 seqNucs = line
+ − 1144 continue
+ − 1145
+ − 1146 if (lineNum % 4 == 3):
+ − 1147 assert (line.startswith("+")), \
+ − 1148 "can't understand line %d:\n%s" % (lineNum,line)
+ − 1149 continue
+ − 1150
+ − 1151 quals = line
+ − 1152 assert (len(quals) == len(seqNucs)), \
+ − 1153 "length mismatch read vs. qualities at line %d" % lineNum
+ − 1154 yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap),quals)
+ − 1155
+ − 1156 assert (lineNum % 4 == 0), \
+ − 1157 "incomplete read at end of file"
+ − 1158
+ − 1159 def sam_sequences(f):
+ − 1160 lineNum = 0
+ − 1161 for line in f:
+ − 1162 lineNum += 1
+ − 1163 line = line.strip()
+ − 1164
+ − 1165 if line.startswith("@"):
+ − 1166 continue
+ − 1167
+ − 1168 columns = line.split("\t")
+ − 1169 seqName = columns[0]
+ − 1170 refName = columns[2]
+ − 1171 pre_s = int(columns[3]) - 1
+ − 1172 cigar = columns[5]
+ − 1173 seqNucs = columns[9]
+ − 1174
+ − 1175 yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap), refName, pre_s, cigar)
+ − 1176
+ − 1177 # sequence_name--
+ − 1178 # Extract the sequence name from a fasta header.
+ − 1179 # $$$ this may need to be improved $$$
+ − 1180
+ − 1181 def sequence_name(s):
+ − 1182 s = s[1:].strip()
+ − 1183 if (s == ""): return ""
+ − 1184 else: return s.split()[0]
+ − 1185
+ − 1186
+ − 1187 # nucleotide_runs--
+ − 1188 # Yield (start,end) for all runs of valid nucleotides in a sequence.
+ − 1189
+ − 1190 def nucleotide_runs(s):
+ − 1191 runs = []
+ − 1192 start = None
+ − 1193 for (ix,nuc) in enumerate(s):
+ − 1194 if (nuc in "ACGT"):
+ − 1195 if (start == None):
+ − 1196 start = ix
+ − 1197 else:
+ − 1198 if (start != None):
+ − 1199 yield (start,ix)
+ − 1200 start = None
+ − 1201
+ − 1202 if (start != None): yield (start,len(s))
+ − 1203
+ − 1204
+ − 1205 # contains_repeat--
+ − 1206 # Determine whether a short sequence contains a repeated element, such as a
+ − 1207 # 6-mer containing a repeated 2-mer (ACACAC) or 3-mer (ACTACT). The repeat
+ − 1208 # must cover the entire sequence, without mismatches.
+ − 1209
+ − 1210 def contains_repeat(kmer):
+ − 1211 kmerLength = len(kmer)
+ − 1212 hasRepeat = False
+ − 1213 rptLen = 1
+ − 1214 while (not hasRepeat) and (2 * rptLen <= kmerLength):
+ − 1215 if (kmerLength % rptLen != 0):
+ − 1216 rptLen += 1
+ − 1217 continue
+ − 1218 isRepeat = True
+ − 1219 for i in xrange(rptLen,kmerLength,rptLen):
+ − 1220 if (kmer[i:i+rptLen] != kmer[:rptLen]):
+ − 1221 isRepeat = False
+ − 1222 break
+ − 1223 if (isRepeat):
+ − 1224 hasRepeat = True
+ − 1225 break
+ − 1226 rptLen += 1
+ − 1227 return hasRepeat
+ − 1228
+ − 1229
+ − 1230 # hash108--
+ − 1231 # Return a 108-bit hash "value" of a string
+ − 1232
+ − 1233 def hash108(s):
+ − 1234 m = md5_new()
+ − 1235 m.update(s)
+ − 1236 return m.hexdigest()[:27]
+ − 1237
+ − 1238
+ − 1239 # float_or_fraction--
+ − 1240 # Convert a string to a number, allowing fractions
+ − 1241
+ − 1242 def float_or_fraction(s):
+ − 1243 if ("/" in s):
+ − 1244 (numer,denom) = s.split("/",1)
+ − 1245 return float(numer)/float(denom)
+ − 1246 else:
+ − 1247 return float(s)
+ − 1248
+ − 1249
+ − 1250 # int_with_unit--
+ − 1251 # Parse a string as an integer, allowing unit suffixes
+ − 1252
+ − 1253 def int_with_unit(s):
+ − 1254 if (s.endswith("K")):
+ − 1255 multiplier = 1000
+ − 1256 s = s[:-1]
+ − 1257 elif (s.endswith("M")):
+ − 1258 multiplier = 1000 * 1000
+ − 1259 s = s[:-1]
+ − 1260 elif (s.endswith("G")):
+ − 1261 multiplier = 1000 * 1000 * 1000
+ − 1262 s = s[:-1]
+ − 1263 else:
+ − 1264 multiplier = 1
+ − 1265
+ − 1266 try: return int(s) * multiplier
+ − 1267 except ValueError: return int(math.ceil(float(s) * multiplier))
+ − 1268
+ − 1269
+ − 1270 if __name__ == "__main__": main()
+ − 1271
