Initial commit to public repo

.gitignore

@@ -0,0 +1,14 @@
+#compile python modules
+*.pyc
+
+#Setuptools distribution folder
+/dist/
+
+#Python egg metadata, regenerated from source files by setuptools
+/*.egg-info
+/*.egg
+
+#
+.DS_Store
+
+

LICENSE

@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2015 Stanford University
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+

MANIFEST.in

@@ -0,0 +1,5 @@
+include README.md
+include LICENSE
+include nucleoatac/standard_vplot.VMat
+include example/*
+

README.md

@@ -0,0 +1,57 @@
+# NucleoATAC
+Python package for calling nucleosomes using ATAC-Seq data.
+
+Non-Python Dependencies:  
+
+* bedGraphToBigWig (from USCS Genome Browser:  http://hgdownload.cse.ucsc.edu/admin/exe/)
+(must be on PATH)
+
+Python Dependencies:
+
+* python 2.7
+* bx-python
+* scipy
+* numpy
+* pysam
+* cython
+* matplotlib
+
+Tested only on Ubuntu OS.  
+
+Installation:
+
+1. Clone repository.
+2. Run `python setup.py install` within main directory
+(use --user flag for user specific install.  See https://docs.python.org/2/install/ for more details)
+
+Necessary files prior to calling nucleosomes:
+* BAM file with aligned data.  Generally this will be filtered for high mapping quality, etc.
+* Bed file with open chromatin regions.  This file must be sorted, and must not contain any overlapping regions!
+* Tn5 bias track in bigWig format (or alternatively, gDNA insertions).  This must span all regions in bed file (with several hundred additional flanking bases)!
+
+
+How to call nucleosomes with mostly default parameters:
+
+```
+nucleoatac run --bed <bed_file> --bam <bam_file> --bias <Tn5_bias_track> --out <output_basename> --cores <num_cores_to_use>
+```
+
+This command is actually running 3 different commands.  For greater versatility, run each separately.  To see options, use help:  
+
+```
+nucleoatac occ --help
+```
+
+```
+nucleoatac vprocess --help
+```
+```
+nucleoatac nuc --help
+```
+
+
+There is an additional nucleoatac call that can be made (to make a V-plot for specific sites):
+```
+nucleoatac vplot --help
+```
+

bias_model/Human.PWM.txt

@@ -0,0 +1,4 @@
+0.840670406462	0.9666090742495	0.929684509606	1.071886557995	0.7997846167775	1.185153593475	0.822059024821	0.7357916346195	0.806541836888	0.643302166961	1.24334576292	1.28021185724	1.097136548525	1.27563558011	0.6086877368615	1.107767288192	0.8366320562835	0.924266259849	0.783703051031	0.8228817857685	1.076371510055
+1.05083196339	0.953504938524	1.29824892591	0.9388735168795	1.61692228526	0.7012930269375	1.122145787175	1.146965873215	1.2830066904	1.43811004046	0.665358216826	0.672207685113	0.8559204963725	0.8365523557765	1.69890306358	0.876213491886	0.9074774433175	1.066666539765	1.1151256444	1.350103850425	1.068801376385
+1.068801376385	1.350103850425	1.1151256444	1.066666539765	0.9074774433175	0.876213491886	1.69890306358	0.8365523557765	0.8559204963725	0.672207685113	0.66710672673	1.43811004046	1.2830066904	1.146965873215	1.122145787175	0.7012930269375	1.61692228526	0.9388735168795	1.29824892591	0.953504938524	1.05083196339
+1.076371510055	0.8228817857685	0.783703051031	0.924266259849	0.8366320562835	1.107767288192	0.6086877368615	1.27563558011	1.097136548525	1.28021185724	1.21947260924	0.643302166961	0.806541836888	0.7357916346195	0.822059024821	1.185153593475	0.7997846167775	1.071886557995	0.929684509606	0.9666090742495	0.840670406462

bias_model/README.txt

@@ -0,0 +1,6 @@
+Human.PWM.txt and Yeast.PWM.txt were derived from mapped read starts from genomic samples from Human and Yeast nextera sequences, respectively.
+
+ScoreSequence.py was used to generate scores for individual chromosomes using the PWM files. 
+
+
+

bias_model/ScoreSequence.py

@@ -0,0 +1,139 @@
+#!/usr/bin/env python
+
+##### IMPORT MODULES #####
+# import necessary for python
+from __future__ import division
+import os
+import sys
+import string
+import datetime
+from optparse import OptionParser
+import numpy as np
+import gzip
+from Bio import SeqIO
+from Bio.Seq import Seq, reverse_complement
+
+#### OPTIONS ####
+# read options from command line
+opts = OptionParser()
+opts.add_option("-i", "--input", default="", 
+help="<inputFile> Accepts FASTA file with single Sequence")
+opts.add_option("-p", "--pwm", default="",
+help="<pwmFile> File with normalized Tn5 transposon preference. Insertion point should be 11th position")
+opts.add_option("-c", "--chrN", default="",
+help="<chrN> name of chromosome (i.e. chr1)")
+opts.add_option("-o", "--out", default="",
+help="<outFile> Base name for Output file")
+options, arguments = opts.parse_args()
+
+##### INPUTS AND OUTPUTS #####
+in_file = options.input
+pwm_file = options.pwm
+chrN=options.chrN
+a=10#bases around transposition center included in motif
+
+if options.input == "":
+    sys.exit("Requires an input file, use --input to designate. Accepts FASTA files.")
+
+if options.out=="":
+	out = os.path.splitext(os.path.basename(in_file))[0]
+else:
+	out=options.out
+
+
+#open PWM file
+pwm=np.loadtxt(pwm_file)
+w=len(pwm[0])
+pwm=np.log(pwm)
+
+##write function to determine the pwm score...
+def scorer_mono(seq,m):
+	#m should be array of size 4 * len(seq)
+	S=0
+	for i in range(len(seq)):
+		if seq[i]=='A':
+			S=S+m[0][i]
+		elif seq[i]=='C':
+			S=S+m[1][i]
+		elif seq[i]=='G':
+			S=S+m[2][i]
+		elif seq[i]=='T':
+			S=S+m[3][i]
+	#logS=np.log(S)
+	return S
+
+def scorer_di(seq,m):
+	#m should be array of size 16 * len(seq)
+	S=0
+	for i in range((len(seq)-1)):
+		if seq[i:(i+2)]=='AA':
+			S=S+m[0][i]
+		elif seq[i:(i+2)]=='AC':
+			S=S+m[1][i]
+		elif seq[i:(i+2)]=='AG':
+			S=S+m[2][i]
+		elif seq[i:(i+2)]=='AT':
+			S=S+m[3][i]
+		elif seq[i:(i+2)]=='CA':
+			S=S+m[4][i]
+		elif seq[i:(i+2)]=='CC':
+			S=S+m[5][i]
+		elif seq[i:(i+2)]=='CG':
+			S=S+m[6][i]
+		elif seq[i:(i+2)]=='CT':
+			S=S+m[7][i]
+		elif seq[i:(i+2)]=='GA':
+			S=S+m[8][i]
+		elif seq[i:(i+2)]=='GC':
+			S=S+m[9][i]
+		elif seq[i:(i+2)]=='GG':
+			S=S+m[10][i]
+		elif seq[i:(i+2)]=='GT':
+			S=S+m[11][i]
+		elif seq[i:(i+2)]=='TA':
+			S=S+m[12][i]
+		elif seq[i:(i+2)]=='TC':
+			S=S+m[13][i]
+		elif seq[i:(i+2)]=='TG':
+			S=S+m[14][i]
+		elif seq[i:(i+2)]=='TT':
+			S=S+m[15][i]
+	#logS=np.log(S)
+	return S
+
+def Scorer(seq,m):
+	if len(m)==4:
+		return scorer_mono(seq,m)
+	elif len(m)==16:
+		return scorer_di(seq,m)
+
+#Get sequence to score 
+
+#if no chromosome is given, just use first sequence in fasta input
+if chrN=='':
+	rec=SeqIO.parse(in_file,'fasta').next()
+	source=rec.seq
+	chrN=rec.id
+else:
+	for rec in SeqIO.parse(in_file,'fasta'):
+		if rec.id==chrN:
+			source=rec.seq
+
+Seq_length=len(source)
+
+#open output file for writing scores
+if options.chrN!='':
+	output=open(out+'.'+chrN+'.Scores.bed','w')
+else:
+	output=open(out+'.Scores.bed','w')
+
+
+for i in range(a,Seq_length-w+a):
+	forwardseq=str(source[i-a:i+w-a]).upper()
+	f=Scorer(forwardseq,pwm)
+	output.write('\t'.join([chrN,str(i),str(i+1),str(f)])+'\n')
+
+
+output.close()
+
+

bias_model/Yeast.PWM.txt

@@ -0,0 +1,4 @@
+0.8249902267275	0.963678911394	0.977394731096	1.094331796415	0.838984693388	1.140060661135	0.8634411405355	0.7351558672495	0.8570707149595	0.657202162157	1.24925773037	1.325312968185	1.11984130425	1.322947131275	0.617562709625	1.052434771076	0.8703150237785	0.9013347604265	0.784602391176	0.79790349428	1.073410201865
+1.080695856195	0.957660337654	1.292056428965	0.8915020935565	1.6430442671	0.718367258387	1.14008567413	1.09189268792	1.22147931095	1.34084133077	0.614913565933	0.687329339637	0.815718787164	0.814494923033	1.696092676885	0.971494577519	0.8253162076685	1.11547971422	1.09140085467	1.426465114795	1.082995843845
+1.082995843845	1.426465114795	1.09140085467	1.11547971422	0.8253162076685	0.971494577519	1.696092676885	0.814494923033	0.815718787164	0.687329339637	0.614126251469	1.34084133077	1.22147931095	1.09189268792	1.14008567413	0.718367258387	1.6430442671	0.8915020935565	1.292056428965	0.957660337654	1.080695856195
+1.073410201865	0.79790349428	0.784602391176	0.9013347604265	0.8703150237785	1.052434771076	0.617562709625	1.322947131275	1.11984130425	1.325312968185	1.22925904295	0.657202162157	0.8570707149595	0.7351558672495	0.8634411405355	1.140060661135	0.838984693388	1.094331796415	0.977394731096	0.963678911394	0.8249902267275