Merge branch 'gamma' to change NFR fit to be based on gamma distribution rather than exponential distribution. This is more robust and corresponds better with observed distributions.

Author: Alicia Schep

example/README.txt

@@ -1,8 +1,14 @@
 
-Example files to test out functions
+#This folder contains Example files to test out functions
+
+#example.bam created by intersecting a bam file with example.slopped.bed (which has 1000bp flank on each side of example.bed)
+
+
+#To test out nucleoatac installation, try:
+
+nucleoatac run --bed example.bed --bam example.bam --fasta sacCer3.fa --out test_example 
+
 
-example.bam created by intersecting a bam file with example.slopped.bed (which has 1000bp flank on each side of example.bed)
 
-example.Scores.bw has scores for ranges in example.slopped.bed
 
 

example/sacCer3.fa

@@ -0,0 +1,17 @@
+chrI	230218	6	50	51
+chrII	813184	234836	50	51
+chrIII	316620	1064292	50	51
+chrIV	1531933	1387252	50	51
+chrV	576874	2949830	50	51
+chrVI	270161	3538249	50	51
+chrVII	1090940	3813822	50	51
+chrVIII	562643	4926590	50	51
+chrIX	439888	5500493	50	51
+chrX	745751	5949185	50	51
+chrXI	666816	6709859	50	51
+chrXII	1078177	7390020	50	51
+chrXIII	924431	8489770	50	51
+chrXIV	784333	9432698	50	51
+chrXV	1091291	10232725	50	51
+chrXVI	948066	11345850	50	51
+chrM	85779	12312884	50	51

example/sacCer3.fa.fai

@@ -15,6 +15,7 @@
 from pyatac.utils import smooth, call_peaks, read_chrom_sizes_from_fasta
 from pyatac.chunkmat2d import FragmentMat2D, BiasMat2D
 from pyatac.bias import InsertionBiasTrack, PWM
+from scipy.special import gamma
 
 
 class FragmentMixDistribution:
@@ -25,24 +26,24 @@ def __init__(self,  lower = 0, upper =2000):
     def getFragmentSizes(self, bamfile, chunklist = None):
         self.fragmentsizes = FragmentSizes(self.lower, self.upper)
         self.fragmentsizes.calculateSizes(bamfile, chunks = chunklist)
-    def modelNFR(self, boundary = 115):
-        """Model NFR distribution with exponential distribution"""
-        b = np.where(self.fragmentsizes.get(self.lower,boundary) == max(self.fragmentsizes.get(self.lower,boundary)))[0][0]+10 + self.lower
-        def exp_pdf(x,*p): #defines the PDF
-            k=p[0]
-            a=p[1]
-            x=x-b
-            return a*k*np.exp(-k*x)
-        x = np.array(range(b,boundary))
-        p0 = (.1,1)
-        coeff, var_matrix = optimize.curve_fit(exp_pdf,x, self.fragmentsizes.get(b,boundary),
+    def modelNFR(self, boundaries = (45,115)):
+        """Model NFR distribution with gamma distribution"""
+        smoothed = smooth(self.fragmentsizes.get(),19, window = "gaussian",sd = 3, mode='same')
+        def gamma_pdf(x,*p): #defines the PDF
+            k = p[0]
+            theta =p[1]
+            a = p[2]
+            return a * x**(k-1) * np.exp(-x/theta) / (theta**k * gamma(k))
+        x = np.arange(boundaries[0],boundaries[1])        
+        p0 = (3, 25, 1)                        
+        coeff, var_matrix = optimize.curve_fit(gamma_pdf,x, self.fragmentsizes.get(boundaries[0],boundaries[1]),
                                                p0=p0)
-        nfr = np.concatenate((self.fragmentsizes.get(self.lower,boundary), exp_pdf(np.array(range(boundary,self.upper)),*coeff)))
+        nfr = np.concatenate((self.fragmentsizes.get(self.lower,boundaries[1]), gamma_pdf(np.array(range(boundaries[1],self.upper)),*coeff)))
         nfr[nfr==0] = min(nfr[nfr!=0])*0.01
         self.nfr_fit = FragmentSizes(self.lower,self.upper, vals = nfr)
-        nuc = np.concatenate((np.zeros(boundary-self.lower),
-                            self.fragmentsizes.get(boundary,self.upper) -
-                            self.nfr_fit.get(boundary,self.upper)))
+        nuc = np.concatenate((np.zeros(boundaries[1]-self.lower),
+                            self.fragmentsizes.get(boundaries[1],self.upper) -
+                            self.nfr_fit.get(boundaries[1],self.upper)))
         nuc[nuc<=0]=min(min(nfr)*0.1,min(nuc[nuc>0])*0.001)
         self.nuc_fit = FragmentSizes(self.lower, self.upper, vals = nuc)
     def plotFits(self,filename=None):

nucleoatac/Occupancy.py

@@ -8,7 +8,7 @@
 
 
 setup(name='NucleoATAC',
-    version='0.2.3',
+    version='0.3.0',
     description='python package for calling nucleosomes with ATAC-Seq',
     classifiers=[
         'Development Status :: 3 - Alpha',

setup.py

@@ -1,7 +1,7 @@
 from unittest import TestCase
 import numpy as np
 import nucleoatac.NucleosomeCalling as Nuc
-import nucleoatac.VMat as V
+import pyatac.VMat as V
 from pyatac.chunkmat2d import BiasMat2D
 from pyatac.chunk import ChunkList
 from pyatac.bias import InsertionBiasTrack

tests/test_var.py

@@ -2,7 +2,7 @@
 import numpy as np
 
 import nucleoatac.NucleosomeCalling as Nuc
-import nucleoatac.VMat as V
+import pyatac.VMat as V
 from pyatac.chunkmat2d import FragmentMat2D
 from pyatac.chunk import ChunkList