| Previous changeset 4:d3dbf9e8a0e2 (2018-01-16) Next changeset 6:aec658f828df (2018-03-29) |
|
Commit message:
Uploaded |
|
added:
peak_calling_script.py |
| b |
| diff -r d3dbf9e8a0e2 -r 6242a111983d peak_calling_script.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/peak_calling_script.py Tue Jan 16 15:15:36 2018 -0500 |
| [ |
| @@ -0,0 +1,83 @@ + +from __future__ import print_function +import sys +import numpy +import math +import random +import csv +import matplotlib.pyplot as plt +import pystache +import json +from sklearn import mixture + +x = [] +y = [] + +toolInput = sys.argv[1] +toolOutput = sys.argv[2] +toolWebsite = sys.argv[3] + +with open(sys.argv[1], 'rb') as csvfile: + spamreader = csv.reader(csvfile, delimiter='\t') + for i, row in enumerate(spamreader): + if i != 0: + x.append(int(row[0])) + y.append(int(row[1])) + +# you have to set this manually to weed out all the noise. Every bit of noise should be below it. +threshold = 20 +rightLimit = 200 + +# unravelling histogram into samples. +samples = [] +for no, value in enumerate([int(round(i)) for i in y]): + if value > threshold and no < rightLimit: + for _ in range(value): + samples.append(no) + +# total number of reads +totalAmp = len(samples) + +# reshaping numpy arrays to indicate that we pass a lot of samples, not a lot of features. +xArray = numpy.array(x).reshape(1, -1) +samplesArray = numpy.array(samples).reshape(-1, 1) + +# learning a gaussian mixture model. +gmm2 = mixture.BayesianGaussianMixture(n_components=2).fit(samplesArray) + +# getting the mean of each gaussian +means = [x[int(round(i[0]))] for i in gmm2.means_] + +# rounding errors +roundErr = [i[0] - int(round(i[0])) for i in gmm2.means_] + +# getting the coverage of each gaussian +weights = gmm2.weights_ + +sampleID = toolOutput + ".html" + +with open(toolOutput, "w") as f: + print("sampleID", file=f, end="\t") + print("Al1", file=f, end="\t") + print("Al2", file=f, end="\t") + print("frac1", file=f, end="\t") + print("frac2", file=f, end="\t") + print(file=f) + print(sampleID, file=f, end="\t") + print(means[0], file=f, end="\t") + print(means[1], file=f, end="\t") + print(weights[0], file=f, end="\t") + print(weights[1], file=f, end="\t") + +template_dir = { + "sampleID": sampleID, + "al1": means[0], + "al2": means[1], + "freq1": weights[0], + "freq2": weights[1], + "x": json.dumps(x), + "y": json.dumps(y) + } +with open(toolWebsite) as wt: + with open(sampleID, "w") as wr: + wr.write(pystache.render(wt.read(), template_dir)) |