neat_genreads: utilities/genMutModel.py annotate

annotate utilities/genMutModel.py @ 10:7d10b55965c9 draft default tip

planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty

author	thondeboer
date	Wed, 16 May 2018 17:02:51 -0400
parents	6e75a84e9338
children

rev	line source
0 6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	1 #!/usr/bin/env python
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	2
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	3 import sys
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	4 import os
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	5 import re
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	6 import bisect
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	7 import pickle
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	8 import argparse
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	9 import numpy as np
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	10 #matplotlib is not used as far as i can see
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	11 #import matplotlib.pyplot as mpl
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	12
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	13 # absolute path to the directory above this script
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	14 SIM_PATH = '/'.join(os.path.realpath(__file__).split('/')[:-2])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	15 sys.path.append(SIM_PATH+'/py/')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	16
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	17 from refFunc import indexRef
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	18
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	19 REF_WHITELIST = [str(n) for n in xrange(1,30)] + ['x','y','X','Y','mt','Mt','MT']
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	20 REF_WHITELIST += ['chr'+n for n in REF_WHITELIST]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	21 VALID_NUCL = ['A','C','G','T']
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	22 VALID_TRINUC = [VALID_NUCL[i]+VALID_NUCL[j]+VALID_NUCL[k] for i in xrange(len(VALID_NUCL)) for j in xrange(len(VALID_NUCL)) for k in xrange(len(VALID_NUCL))]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	23 # if parsing a dbsnp vcf, and no CAF= is found in info tag, use this as default val for population freq
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	24 VCF_DEFAULT_POP_FREQ = 0.00001
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	25
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	26
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	27 #########################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	28 # VARIOUS HELPER FUNCTIONS #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	29 #########################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	30
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	31
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	32 # given a reference index, grab the sequence string of a specified reference
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	33 def getChrFromFasta(refPath,ref_inds,chrName):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	34 for i in xrange(len(ref_inds)):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	35 if ref_inds[i][0] == chrName:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	36 ref_inds_i = ref_inds[i]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	37 break
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	38 refFile = open(refPath,'r')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	39 refFile.seek(ref_inds_i[1])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	40 myDat = ''.join(refFile.read(ref_inds_i[2]-ref_inds_i[1]).split('\n'))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	41 return myDat
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	42
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	43 # cluster a sorted list
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	44 def clusterList(l,delta):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	45 outList = [[l[0]]]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	46 prevVal = l[0]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	47 currentInd = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	48 for n in l[1:]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	49 if n-prevVal <= delta:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	50 outList[currentInd].append(n)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	51 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	52 currentInd += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	53 outList.append([])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	54 outList[currentInd].append(n)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	55 prevVal = n
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	56 return outList
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	57
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	58 def list_2_countDict(l):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	59 cDict = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	60 for n in l:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	61 if n not in cDict:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	62 cDict[n] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	63 cDict[n] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	64 return cDict
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	65
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	66 def getBedTracks(fn):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	67 f = open(fn,'r')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	68 trackDict = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	69 for line in f:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	70 splt = line.strip().split('\t')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	71 if splt[0] not in trackDict:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	72 trackDict[splt[0]] = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	73 trackDict[splt[0]].extend([int(splt[1]),int(splt[2])])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	74 f.close()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	75 return trackDict
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	76
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	77 def getTrackLen(trackDict):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	78 totSum = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	79 for k in trackDict.keys():
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	80 for i in xrange(0,len(trackDict[k]),2):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	81 totSum += trackDict[k][i+1] - trackDict[k][i] + 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	82 return totSum
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	83
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	84 def isInBed(track,ind):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	85 myInd = bisect.bisect(track,ind)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	86 if myInd&1:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	87 return True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	88 if myInd < len(track):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	89 if track[myInd-1] == ind:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	90 return True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	91 return False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	92
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	93 ## return the mean distance to the median of a cluster
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	94 #def mean_dist_from_median(c):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	95 # centroid = np.median([n for n in c])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	96 # dists = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	97 # for n in c:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	98 # dists.append(abs(n-centroid))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	99 # return np.mean(dists)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	100 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	101 ## get median value from counting dictionary
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	102 #def quick_median(countDict):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	103 # midPoint = sum(countDict.values())/2
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	104 # mySum = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	105 # myInd = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	106 # sk = sorted(countDict.keys())
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	107 # while mySum < midPoint:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	108 # mySum += countDict[sk[myInd]]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	109 # if mySum >= midPoint:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	110 # break
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	111 # myInd += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	112 # return myInd
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	113 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	114 ## get median deviation from median of counting dictionary
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	115 #def median_deviation_from_median(countDict):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	116 # myMedian = quick_median(countDict)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	117 # deviations = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	118 # for k in sorted(countDict.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	119 # d = abs(k-myMedian)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	120 # deviations[d] = countDict[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	121 # return quick_median(deviations)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	122
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	123
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	124 #################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	125 # PARSE INPUT OPTIONS #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	126 #################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	127
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	128
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	129 parser = argparse.ArgumentParser(description='genMutModel.py')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	130 parser.add_argument('-r', type=str, required=True, metavar='<str>', help="* ref.fa")
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	131 parser.add_argument('-m', type=str, required=True, metavar='<str>', help="* mutations.tsv [.vcf]")
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	132 parser.add_argument('-o', type=str, required=True, metavar='<str>', help="* output.p")
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	133 parser.add_argument('-bi', type=str, required=False, metavar='<str>', default=None, help="only_use_these_regions.bed")
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	134 parser.add_argument('-be', type=str, required=False, metavar='<str>', default=None, help="exclude_these_regions.bed")
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	135 parser.add_argument('--save-trinuc', required=False,action='store_true', default=False, help='save trinuc counts for ref')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	136 parser.add_argument('--no-whitelist',required=False,action='store_true', default=False, help='allow any non-standard ref')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	137 parser.add_argument('--skip-common', required=False,action='store_true', default=False, help='do not save common snps + high mut regions')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	138 args = parser.parse_args()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	139 (REF, TSV, OUT_PICKLE, SAVE_TRINUC, NO_WHITELIST, SKIP_COMMON) = (args.r, args.m, args.o, args.save_trinuc, args.no_whitelist, args.skip_common)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	140
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	141 MYBED = None
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	142 if args.bi != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	143 print 'only considering variants in specified bed regions...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	144 MYBED = (getBedTracks(args.bi),True)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	145 elif args.be != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	146 print 'only considering variants outside of specified bed regions...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	147 MYBED = (getBedTracks(args.be),False)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	148
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	149 if TSV[-4:] == '.vcf':
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	150 IS_VCF = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	151 elif TSV[-4:] == '.tsv':
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	152 IS_VCF = False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	153 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	154 print '\nError: Unknown format for mutation input.\n'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	155 exit(1)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	156
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	157
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	158 #####################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	159 # main() #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	160 #####################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	161
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	162
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	163 def main():
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	164
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	165 ref_inds = indexRef(REF)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	166 refList = [n[0] for n in ref_inds]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	167
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	168 # how many times do we observe each trinucleotide in the reference (and input bed region, if present)?
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	169 TRINUC_REF_COUNT = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	170 TRINUC_BED_COUNT = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	171 printBedWarning = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	172 # [(trinuc_a, trinuc_b)] = # of times we observed a mutation from trinuc_a into trinuc_b
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	173 TRINUC_TRANSITION_COUNT = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	174 # total count of SNPs
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	175 SNP_COUNT = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	176 # overall SNP transition probabilities
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	177 SNP_TRANSITION_COUNT = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	178 # total count of indels, indexed by length
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	179 INDEL_COUNT = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	180 # tabulate how much non-N reference sequence we've eaten through
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	181 TOTAL_REFLEN = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	182 # detect variants that occur in a significant percentage of the input samples (pos,ref,alt,pop_fraction)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	183 COMMON_VARIANTS = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	184 # tabulate how many unique donors we've encountered (this is useful for identifying common variants)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	185 TOTAL_DONORS = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	186 # identify regions that have significantly higher local mutation rates than the average
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	187 HIGH_MUT_REGIONS = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	188
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	189 # load and process variants in each reference sequence individually, for memory reasons...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	190 for refName in refList:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	191
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	192 if (refName not in REF_WHITELIST) and (not NO_WHITELIST):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	193 print refName,'is not in our whitelist, skipping...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	194 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	195
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	196 print 'reading reference "'+refName+'"...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	197 refSequence = getChrFromFasta(REF,ref_inds,refName).upper()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	198 TOTAL_REFLEN += len(refSequence) - refSequence.count('N')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	199
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	200 # list to be used for counting variants that occur multiple times in file (i.e. in multiple samples)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	201 VDAT_COMMON = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	202
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	203
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	204 """ ##########################################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	205 ### COUNT TRINUCLEOTIDES IN REF ###
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	206 ########################################################################## """
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	207
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	208
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	209 if MYBED != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	210 if printBedWarning:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	211 print "since you're using a bed input, we have to count trinucs in bed region even if you specified a trinuc count file for the reference..."
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	212 printBedWarning = False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	213 if refName in MYBED[0]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	214 refKey = refName
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	215 elif ('chr' in refName) and (refName not in MYBED[0]) and (refName[3:] in MYBED[0]):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	216 refKey = refName[3:]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	217 elif ('chr' not in refName) and (refName not in MYBED[0]) and ('chr'+refName in MYBED[0]):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	218 refKey = 'chr'+refName
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	219 if refKey in MYBED[0]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	220 subRegions = [(MYBED[0][refKey][n],MYBED[0][refKey][n+1]) for n in xrange(0,len(MYBED[0][refKey]),2)]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	221 for sr in subRegions:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	222 for i in xrange(sr[0],sr[1]+1-2):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	223 trinuc = refSequence[i:i+3]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	224 if not trinuc in VALID_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	225 continue # skip if trinuc contains invalid characters, or not in specified bed region
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	226 if trinuc not in TRINUC_BED_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	227 TRINUC_BED_COUNT[trinuc] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	228 TRINUC_BED_COUNT[trinuc] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	229
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	230 if not os.path.isfile(REF+'.trinucCounts'):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	231 print 'counting trinucleotides in reference...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	232 for i in xrange(len(refSequence)-2):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	233 if i%1000000 == 0 and i > 0:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	234 print i,'/',len(refSequence)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	235 #break
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	236 trinuc = refSequence[i:i+3]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	237 if not trinuc in VALID_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	238 continue # skip if trinuc contains invalid characters
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	239 if trinuc not in TRINUC_REF_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	240 TRINUC_REF_COUNT[trinuc] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	241 TRINUC_REF_COUNT[trinuc] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	242 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	243 print 'skipping trinuc counts (for whole reference) because we found a file...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	244
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	245
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	246 """ ##########################################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	247 ### READ INPUT VARIANTS ###
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	248 ########################################################################## """
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	249
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	250
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	251 print 'reading input variants...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	252 f = open(TSV,'r')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	253 isFirst = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	254 for line in f:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	255
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	256 if IS_VCF and line[0] == '#':
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	257 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	258 if isFirst:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	259 if IS_VCF:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	260 # hard-code index values based on expected columns in vcf
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	261 (c1,c2,c3,m1,m2,m3) = (0,1,1,3,3,4)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	262 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	263 # determine columns of fields we're interested in
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	264 splt = line.strip().split('\t')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	265 (c1,c2,c3) = (splt.index('chromosome'),splt.index('chromosome_start'),splt.index('chromosome_end'))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	266 (m1,m2,m3) = (splt.index('reference_genome_allele'),splt.index('mutated_from_allele'),splt.index('mutated_to_allele'))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	267 (d_id) = (splt.index('icgc_donor_id'))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	268 isFirst = False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	269 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	270
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	271 splt = line.strip().split('\t')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	272 # we have -1 because tsv/vcf coords are 1-based, and our reference string index is 0-based
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	273 [chrName,chrStart,chrEnd] = [splt[c1],int(splt[c2])-1,int(splt[c3])-1]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	274 [allele_ref,allele_normal,allele_tumor] = [splt[m1].upper(),splt[m2].upper(),splt[m3].upper()]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	275 if IS_VCF:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	276 if len(allele_ref) != len(allele_tumor):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	277 # indels in tsv don't include the preserved first nucleotide, so lets trim the vcf alleles
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	278 [allele_ref,allele_normal,allele_tumor] = [allele_ref[1:],allele_normal[1:],allele_tumor[1:]]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	279 if not allele_ref: allele_ref = '-'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	280 if not allele_normal: allele_normal = '-'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	281 if not allele_tumor: allele_tumor = '-'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	282 # if alternate alleles are present, lets just ignore this variant. I may come back and improve this later
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	283 if ',' in allele_tumor:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	284 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	285 vcf_info = ';'+splt[7]+';'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	286 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	287 [donor_id] = [splt[d_id]]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	288
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	289 # if we encounter a multi-np (i.e. 3 nucl --> 3 different nucl), let's skip it for now...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	290 if ('-' not in allele_normal and '-' not in allele_tumor) and (len(allele_normal) > 1 or len(allele_tumor) > 1):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	291 print 'skipping a complex variant...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	292 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	293
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	294 # to deal with '1' vs 'chr1' references, manually change names. this is hacky and bad.
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	295 if 'chr' not in chrName:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	296 chrName = 'chr'+chrName
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	297 if 'chr' not in refName:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	298 refName = 'chr'+refName
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	299 # skip irrelevant variants
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	300 if chrName != refName:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	301 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	302
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	303 # if variant is outside the regions we're interested in (if specified), skip it...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	304 if MYBED != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	305 refKey = refName
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	306 if not refKey in MYBED[0] and refKey[3:] in MYBED[0]: # account for 1 vs chr1, again...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	307 refKey = refKey[3:]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	308 if refKey not in MYBED[0]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	309 inBed = False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	310 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	311 inBed = isInBed(MYBED[0][refKey],chrStart)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	312 if inBed != MYBED[1]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	313 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	314
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	315 # we want only snps
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	316 # so, no '-' characters allowed, and chrStart must be same as chrEnd
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	317 if '-' not in allele_normal and '-' not in allele_tumor and chrStart == chrEnd:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	318 trinuc_ref = refSequence[chrStart-1:chrStart+2]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	319 if not trinuc_ref in VALID_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	320 continue # skip ref trinuc with invalid characters
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	321 # only consider positions where ref allele in tsv matches the nucleotide in our reference
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	322 if allele_ref == trinuc_ref[1]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	323 trinuc_normal = refSequence[chrStart-1] + allele_normal + refSequence[chrStart+1]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	324 trinuc_tumor = refSequence[chrStart-1] + allele_tumor + refSequence[chrStart+1]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	325 if not trinuc_normal in VALID_TRINUC or not trinuc_tumor in VALID_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	326 continue # skip if mutation contains invalid char
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	327 key = (trinuc_normal,trinuc_tumor)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	328 if key not in TRINUC_TRANSITION_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	329 TRINUC_TRANSITION_COUNT[key] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	330 TRINUC_TRANSITION_COUNT[key] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	331 SNP_COUNT += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	332 key2 = (allele_normal,allele_tumor)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	333 if key2 not in SNP_TRANSITION_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	334 SNP_TRANSITION_COUNT[key2] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	335 SNP_TRANSITION_COUNT[key2] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	336
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	337 if IS_VCF:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	338 myPopFreq = VCF_DEFAULT_POP_FREQ
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	339 if ';CAF=' in vcf_info:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	340 cafStr = re.findall(r";CAF=.*?(?=;)",vcf_info)[0]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	341 if ',' in cafStr:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	342 myPopFreq = float(cafStr[5:].split(',')[1])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	343 VDAT_COMMON.append((chrStart,allele_ref,allele_normal,allele_tumor,myPopFreq))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	344 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	345 VDAT_COMMON.append((chrStart,allele_ref,allele_normal,allele_tumor))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	346 TOTAL_DONORS[donor_id] = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	347 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	348 print '\nError: ref allele in variant call does not match reference.\n'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	349 exit(1)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	350
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	351 # now let's look for indels...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	352 if '-' in allele_normal: len_normal = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	353 else: len_normal = len(allele_normal)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	354 if '-' in allele_tumor: len_tumor = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	355 else: len_tumor = len(allele_tumor)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	356 if len_normal != len_tumor:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	357 indel_len = len_tumor - len_normal
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	358 if indel_len not in INDEL_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	359 INDEL_COUNT[indel_len] = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	360 INDEL_COUNT[indel_len] += 1
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	361
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	362 if IS_VCF:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	363 myPopFreq = VCF_DEFAULT_POP_FREQ
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	364 if ';CAF=' in vcf_info:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	365 cafStr = re.findall(r";CAF=.*?(?=;)",vcf_info)[0]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	366 if ',' in cafStr:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	367 myPopFreq = float(cafStr[5:].split(',')[1])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	368 VDAT_COMMON.append((chrStart,allele_ref,allele_normal,allele_tumor,myPopFreq))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	369 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	370 VDAT_COMMON.append((chrStart,allele_ref,allele_normal,allele_tumor))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	371 TOTAL_DONORS[donor_id] = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	372 f.close()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	373
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	374 # if we didn't find anything, skip ahead along to the next reference sequence
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	375 if not len(VDAT_COMMON):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	376 print 'Found no variants for this reference, moving along...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	377 continue
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	378
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	379 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	380 # identify common mutations
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	381 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	382 percentile_var = 95
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	383 if IS_VCF:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	384 minVal = np.percentile([n[4] for n in VDAT_COMMON],percentile_var)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	385 for k in sorted(VDAT_COMMON):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	386 if k[4] >= minVal:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	387 COMMON_VARIANTS.append((refName,k[0],k[1],k[3],k[4]))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	388 VDAT_COMMON = {(n[0],n[1],n[2],n[3]):n[4] for n in VDAT_COMMON}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	389 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	390 N_DONORS = len(TOTAL_DONORS)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	391 VDAT_COMMON = list_2_countDict(VDAT_COMMON)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	392 minVal = int(np.percentile(VDAT_COMMON.values(),percentile_var))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	393 for k in sorted(VDAT_COMMON.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	394 if VDAT_COMMON[k] >= minVal:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	395 COMMON_VARIANTS.append((refName,k[0],k[1],k[3],VDAT_COMMON[k]/float(N_DONORS)))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	396
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	397 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	398 # identify areas that have contained significantly higher random mutation rates
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	399 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	400 dist_thresh = 2000
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	401 percentile_clust = 97
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	402 qptn = 1000
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	403 # identify regions with disproportionately more variants in them
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	404 VARIANT_POS = sorted([n[0] for n in VDAT_COMMON.keys()])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	405 clustered_pos = clusterList(VARIANT_POS,dist_thresh)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	406 byLen = [(len(clustered_pos[i]),min(clustered_pos[i]),max(clustered_pos[i]),i) for i in xrange(len(clustered_pos))]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	407 #byLen = sorted(byLen,reverse=True)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	408 #minLen = int(np.percentile([n[0] for n in byLen],percentile_clust))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	409 #byLen = [n for n in byLen if n[0] >= minLen]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	410 candidate_regions = []
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	411 for n in byLen:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	412 bi = int((n[1]-dist_thresh)/float(qptn))*qptn
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	413 bf = int((n[2]+dist_thresh)/float(qptn))*qptn
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	414 candidate_regions.append((n[0]/float(bf-bi),max([0,bi]),min([len(refSequence),bf])))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	415 minVal = np.percentile([n[0] for n in candidate_regions],percentile_clust)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	416 for n in candidate_regions:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	417 if n[0] >= minVal:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	418 HIGH_MUT_REGIONS.append((refName,n[1],n[2],n[0]))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	419 # collapse overlapping regions
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	420 for i in xrange(len(HIGH_MUT_REGIONS)-1,0,-1):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	421 if HIGH_MUT_REGIONS[i-1][2] >= HIGH_MUT_REGIONS[i][1] and HIGH_MUT_REGIONS[i-1][0] == HIGH_MUT_REGIONS[i][0]:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	422 avgMutRate = 0.5HIGH_MUT_REGIONS[i-1][3]+0.5HIGH_MUT_REGIONS[i][3] # not accurate, but I'm lazy
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	423 HIGH_MUT_REGIONS[i-1] = (HIGH_MUT_REGIONS[i-1][0], HIGH_MUT_REGIONS[i-1][1], HIGH_MUT_REGIONS[i][2], avgMutRate)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	424 del HIGH_MUT_REGIONS[i]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	425
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	426 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	427 # if we didn't count ref trinucs because we found file, read in ref counts from file now
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	428 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	429 if os.path.isfile(REF+'.trinucCounts'):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	430 print 'reading pre-computed trinuc counts...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	431 f = open(REF+'.trinucCounts','r')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	432 for line in f:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	433 splt = line.strip().split('\t')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	434 TRINUC_REF_COUNT[splt[0]] = int(splt[1])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	435 f.close()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	436 # otherwise, save trinuc counts to file, if desired
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	437 elif SAVE_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	438 if MYBED != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	439 print 'unable to save trinuc counts to file because using input bed region...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	440 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	441 print 'saving trinuc counts to file...'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	442 f = open(REF+'.trinucCounts','w')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	443 for trinuc in sorted(TRINUC_REF_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	444 f.write(trinuc+'\t'+str(TRINUC_REF_COUNT[trinuc])+'\n')
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	445 f.close()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	446
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	447 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	448 # if using an input bed region, make necessary adjustments to trinuc ref counts based on the bed region trinuc counts
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	449 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	450 if MYBED != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	451 if MYBED[1] == True: # we are restricting our attention to bed regions, so ONLY use bed region trinuc counts
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	452 TRINUC_REF_COUNT = TRINUC_BED_COUNT
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	453 else: # we are only looking outside bed regions, so subtract bed region trinucs from entire reference trinucs
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	454 for k in TRINUC_REF_COUNT.keys():
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	455 if k in TRINUC_BED_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	456 TRINUC_REF_COUNT[k] -= TRINUC_BED_COUNT[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	457
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	458 # if for some reason we didn't find any valid input variants, exit gracefully...
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	459 totalVar = SNP_COUNT + sum(INDEL_COUNT.values())
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	460 if totalVar == 0:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	461 print '\nError: No valid variants were found, model could not be created. (Are you using the correct reference?)\n'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	462 exit(1)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	463
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	464 """ ##########################################################################
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	465 ### COMPUTE PROBABILITIES ###
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	466 ########################################################################## """
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	467
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	468
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	469 #for k in sorted(TRINUC_REF_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	470 # print k, TRINUC_REF_COUNT[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	471 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	472 #for k in sorted(TRINUC_TRANSITION_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	473 # print k, TRINUC_TRANSITION_COUNT[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	474
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	475 # frequency that each trinuc mutated into anything else
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	476 TRINUC_MUT_PROB = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	477 # frequency that a trinuc mutates into another trinuc, given that it mutated
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	478 TRINUC_TRANS_PROBS = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	479 # frequency of snp transitions, given a snp occurs.
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	480 SNP_TRANS_FREQ = {}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	481
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	482 for trinuc in sorted(TRINUC_REF_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	483 myCount = 0
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	484 for k in sorted(TRINUC_TRANSITION_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	485 if k[0] == trinuc:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	486 myCount += TRINUC_TRANSITION_COUNT[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	487 TRINUC_MUT_PROB[trinuc] = myCount / float(TRINUC_REF_COUNT[trinuc])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	488 for k in sorted(TRINUC_TRANSITION_COUNT.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	489 if k[0] == trinuc:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	490 TRINUC_TRANS_PROBS[k] = TRINUC_TRANSITION_COUNT[k] / float(myCount)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	491
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	492 for n1 in VALID_NUCL:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	493 rollingTot = sum([SNP_TRANSITION_COUNT[(n1,n2)] for n2 in VALID_NUCL if (n1,n2) in SNP_TRANSITION_COUNT])
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	494 for n2 in VALID_NUCL:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	495 key2 = (n1,n2)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	496 if key2 in SNP_TRANSITION_COUNT:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	497 SNP_TRANS_FREQ[key2] = SNP_TRANSITION_COUNT[key2] / float(rollingTot)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	498
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	499 # compute average snp and indel frequencies
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	500 SNP_FREQ = SNP_COUNT/float(totalVar)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	501 AVG_INDEL_FREQ = 1.-SNP_FREQ
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	502 INDEL_FREQ = {k:(INDEL_COUNT[k]/float(totalVar))/AVG_INDEL_FREQ for k in INDEL_COUNT.keys()}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	503 if MYBED != None:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	504 if MYBED[1] == True:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	505 AVG_MUT_RATE = totalVar/float(getTrackLen(MYBED[0]))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	506 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	507 AVG_MUT_RATE = totalVar/float(TOTAL_REFLEN - getTrackLen(MYBED[0]))
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	508 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	509 AVG_MUT_RATE = totalVar/float(TOTAL_REFLEN)
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	510
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	511 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	512 # if values weren't found in data, appropriately append null entries
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	513 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	514 printTrinucWarning = False
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	515 for trinuc in VALID_TRINUC:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	516 trinuc_mut = [trinuc[0]+n+trinuc[2] for n in VALID_NUCL if n != trinuc[1]]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	517 if trinuc not in TRINUC_MUT_PROB:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	518 TRINUC_MUT_PROB[trinuc] = 0.
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	519 printTrinucWarning = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	520 for trinuc2 in trinuc_mut:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	521 if (trinuc,trinuc2) not in TRINUC_TRANS_PROBS:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	522 TRINUC_TRANS_PROBS[(trinuc,trinuc2)] = 0.
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	523 printTrinucWarning = True
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	524 if printTrinucWarning:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	525 print 'Warning: Some trinucleotides transitions were not encountered in the input dataset, probabilities of 0.0 have been assigned to these events.'
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	526
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	527 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	528 # print some stuff
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	529 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	530 for k in sorted(TRINUC_MUT_PROB.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	531 print 'p('+k+' mutates) =',TRINUC_MUT_PROB[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	532
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	533 for k in sorted(TRINUC_TRANS_PROBS.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	534 print 'p('+k[0]+' --> '+k[1]+' \| '+k[0]+' mutates) =',TRINUC_TRANS_PROBS[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	535
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	536 for k in sorted(INDEL_FREQ.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	537 if k > 0:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	538 print 'p(ins length = '+str(abs(k))+' \| indel occurs) =',INDEL_FREQ[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	539 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	540 print 'p(del length = '+str(abs(k))+' \| indel occurs) =',INDEL_FREQ[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	541
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	542 for k in sorted(SNP_TRANS_FREQ.keys()):
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	543 print 'p('+k[0]+' --> '+k[1]+' \| SNP occurs) =',SNP_TRANS_FREQ[k]
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	544
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	545 #for n in COMMON_VARIANTS:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	546 # print n
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	547
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	548 #for n in HIGH_MUT_REGIONS:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	549 # print n
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	550
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	551 print 'p(snp) =',SNP_FREQ
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	552 print 'p(indel) =',AVG_INDEL_FREQ
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	553 print 'overall average mut rate:',AVG_MUT_RATE
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	554 print 'total variants processed:',totalVar
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	555
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	556 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	557 # save variables to file
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	558 #
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	559 if SKIP_COMMON:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	560 OUT_DICT = {'AVG_MUT_RATE':AVG_MUT_RATE,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	561 'SNP_FREQ':SNP_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	562 'SNP_TRANS_FREQ':SNP_TRANS_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	563 'INDEL_FREQ':INDEL_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	564 'TRINUC_MUT_PROB':TRINUC_MUT_PROB,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	565 'TRINUC_TRANS_PROBS':TRINUC_TRANS_PROBS}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	566 else:
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	567 OUT_DICT = {'AVG_MUT_RATE':AVG_MUT_RATE,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	568 'SNP_FREQ':SNP_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	569 'SNP_TRANS_FREQ':SNP_TRANS_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	570 'INDEL_FREQ':INDEL_FREQ,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	571 'TRINUC_MUT_PROB':TRINUC_MUT_PROB,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	572 'TRINUC_TRANS_PROBS':TRINUC_TRANS_PROBS,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	573 'COMMON_VARIANTS':COMMON_VARIANTS,
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	574 'HIGH_MUT_REGIONS':HIGH_MUT_REGIONS}
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	575 pickle.dump( OUT_DICT, open( OUT_PICKLE, "wb" ) )
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	576
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	577
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	578 if __name__ == "__main__":
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	579 main()
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	580
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	581
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	582
6e75a84e9338 planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty thondeboer parents: diff changeset	583

Mercurial > repos > thondeboer > neat_genreads

annotate utilities/genMutModel.py @ 10:7d10b55965c9 draft default tip