| Previous changeset 4:4a4d2b78eb55 (2016-08-06) Next changeset 6:afb015481d85 (2016-08-06) |
|
Commit message:
Uploaded cnvsim library package |
|
added:
cnvsim/ART/README.md cnvsim/ART/art_illumina cnvsim/Wessim/Prep_BlatSearch.py cnvsim/Wessim/Prep_Probe2Fa.py cnvsim/Wessim/README.md cnvsim/Wessim/Wessim1.py cnvsim/Wessim/Wessim2.py cnvsim/Wessim/__sub_wessim1.py cnvsim/Wessim/__sub_wessim2.py cnvsim/Wessim/lib/S0293689_Probes.txt cnvsim/Wessim/lib/mvnTable.txt cnvsim/Wessim/models/ill100v4_p.gzip cnvsim/Wessim/models/ill100v4_s.gzip cnvsim/Wessim/models/ill100v5_p.gzip cnvsim/Wessim/models/ill100v5_s.gzip cnvsim/Wessim/models/r454ti_s.gzip cnvsim/__init__.py cnvsim/__init__.pyc cnvsim/exome_simulator.py cnvsim/exome_simulator.pyc cnvsim/fileio.py cnvsim/fileio.pyc cnvsim/genome_simulator.py cnvsim/genome_simulator.pyc |
|
removed:
cnv-sim.py cnv_sim.xml |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnv-sim.py --- a/cnv-sim.py Sat Aug 06 15:25:34 2016 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,79 +0,0 @@ -#!/usr/bin/python - -__author__ = 'Abdelrahman Hosny' - -import os.path -import datetime -import argparse -import shutil - -from cnvsim.fileio import * -from cnvsim.exome_simulator import * -from cnvsim.genome_simulator import * - -def log(message): - print '[CNV SIM {:%Y-%m-%d %H:%M:%S}'.format(datetime.datetime.now()) + "] " + message - -def main(): - parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument("simulation_type", type=str, choices=['genome', 'exome'], \ - help="simulate copy number variations in whole genome or exome regions") - parser.add_argument("genome", type=file, \ - help="path to the referece genome file in FASTA format ") - parser.add_argument("target", type=file, nargs='?', default=None, \ - help="path to the target regions file in BED format (if using exome)") - - parser.add_argument("-o", "--output_dir_name",type=str, default="test", \ - help="a name to be used to create the output directory (overrides existing directory with the same name).") - parser.add_argument("-n", "--n_reads", type=int, default=10000, \ - help="total number of reads without variations") - parser.add_argument("-l", "--read_length", type=int, default=100, \ - help="read length (bp)") - parser.add_argument("--cnv_list", type=file, default=None, \ - help="path to a CNV list file in BED format chr | start | end | variation. If not passed, it is randomly generated using CNV list parameters below") - - cnv_sim_group = parser.add_argument_group('CNV list parameters', "parameters to be used if CNV list is not passed") - cnv_sim_group.add_argument("-g", "--regions_count", type=int, default=30, \ - help="number of CNV regions to be randomly generated") - cnv_sim_group.add_argument("-a", "--amplifications", type=float, default=0.30, \ - help="percentage of amplifications in range [0.0: 1.0].") - cnv_sim_group.add_argument("-d", "--deletions", type=float, default=0.20, \ - help="percentage of deletions in range [0.0: 1.0].") - cnv_sim_group.add_argument("-min", "--minimum", type=float, default=3, \ - help="minimum number of amplifications/deletions introduced") - cnv_sim_group.add_argument("-max", "--maximum", type=float, default=10, \ - help="maximum number of amplifications/deletions introduced") - - args = parser.parse_args() - - simulation_parameters = {} - simulation_parameters['type'] = args.simulation_type - simulation_parameters['genome_file'] = args.genome.name - if args.target is not None: - simulation_parameters['target_file'] = args.target.name - else: - simulation_parameters['target_file'] = None - simulation_parameters['output_dir'] = os.path.join(os.getcwd(), args.output_dir_name) - simulation_parameters['number_of_reads'] = args.n_reads - simulation_parameters['read_length'] = args.read_length - if args.cnv_list is not None: - simulation_parameters['cnv_list_file'] = args.cnv_list.name - else: - simulation_parameters['cnv_list_file'] = None - simulation_parameters['tmp_dir'] = os.path.join(os.getcwd(), args.output_dir_name , "tmp") - - cnv_list_parameters = {} - cnv_list_parameters['regions_count'] = args.regions_count - cnv_list_parameters['amplifications'] = args.amplifications - cnv_list_parameters['deletions'] = args.deletions - cnv_list_parameters['minimum_variations'] = args.minimum - cnv_list_parameters['maximum_variations'] = args.maximum - - if simulation_parameters['type'] == 'genome': - simulate_genome_cnv(simulation_parameters, cnv_list_parameters) - else: - simulate_exome_cnv(simulation_parameters, cnv_list_parameters) - - -if __name__ == '__main__': - main() \ No newline at end of file |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnv_sim.xml --- a/cnv_sim.xml Sat Aug 06 15:25:34 2016 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,56 +0,0 @@ -<tool id="cnvsim" name="Simulate CNV" version="0.9.0"> - <description>in NGS short reads </description> - <command interpreter="python" detect_errors="exit_code"> - #if $type.simulation_type=="whole_genome" - cnv-sim.py -n $number_of_reads -l $read_length -g $regions_count -a $amplifications -d $deletions -min $minimum -max $maximum genome $reference - #else - cnv-sim.py -n $number_of_reads -l $read_length -g $regions_count -a $amplifications -d $deletions -min $minimum -max $maximum exome $reference $target - #end if - </command> - <inputs> - <conditional name="type"> - <param name="simulation_type" type="select" label="Simulation Type"> - <option value="whole_genome">CNV Simulation in Whole Genome</option> - <option value="whole_exome">CNV Simulation in Exome Regions</option> - </param> - <when value="whole_genome"> - <param format="fasta" name="reference" type="data" label="Reference Genome" help="Reference genome to generate reads from"/> - </when> - <when value="whole_exome"> - <param format="fasta" name="reference" type="data" label="Reference Genome" help="Reference genome to generate reads from"/> - <param format="bed" name="target" type="data" label="Target Regions" help="a list of exonic regions"/> - </when> - </conditional> - <param name="number_of_reads" type="integer" value="100000" label="Number of Reads" help="determines the number of reads to be generated for the control and simulated CNV (approximately)"/> - <param name="read_length" type="integer" value="100" label="Read Length (bp)" help="determines the read length fot the generated FASTQ files"/> - - <param name="regions_count" type="integer" value="30" label="Regions Count" help="determines how many randomly-generated regions will show CNVs"/> - <param name="amplifications" type="float" value="0.5" label="Percentage of amplifications" help="determines what fraction of the regions will show amplifications (range: 0.0-1.0)"/> - <param name="deletions" type="float" value="0.2" label="Percentage of deletions" help="determines what fraction of the regions will show deletions (range: 0.0-1.0)"/> - <param name="minimum" type="integer" value="3" label="Variation minimum" help="determines the minumum number of amplifications/deletions introduced in each region"/> - <param name="maximum" type="integer" value="10" label="Variation maximum" help="determines the maximum number of amplifications/deletions introduced in each region"/> - - </inputs> - <outputs> - <data format="bed" name="cnv_list" from_work_dir="test/CNVList.bed" label="CNV List from ${tool.name} on ${on_string}"/> - <data format="fastq" name="control_1" from_work_dir="test/control_1.fastq" label="Control reads 1 from ${tool.name} on ${on_string}"/> - <data format="fastq" name="control_2" from_work_dir="test/control_2.fastq" label="Control reads 2 from ${tool.name} on ${on_string}"/> - <data format="fastq" name="cnv_1" from_work_dir="test/cnv_1.fastq" label="CNV reads 1 from ${tool.name} on ${on_string}"/> - <data format="fastq" name="cnv_2" from_work_dir="test/cnv_2.fastq" label="CNV reads 2 from ${tool.name} on ${on_string}"/> - </outputs> - <help><![CDATA[ - .. class:: infomark - '''TIP''' This tool requires *fasta* format. - ---- - **CNV Simulator** - In genomics, Copy Number Variations (CNVs) is a type of structural variation in a genome where sections of the genome are repeated. The number if repetitions (duplications) varies between individuals in the human population. - - The Copy Number Variation Simulator (CNV Sim) is a tool used to generate a set of artificial DNA fragments for Next Generation Sequencing (NGS) read simulation. When aligned back to the reference genome, the artificial generated reads show variations in the CNV regions. Variations can be either amplifications of deletions. - - CNV-Sim offers two types of simulation: - 1. CNV simulation in whole genome. CNV-Sim wraps the functionality of ART to introduce variations in the genome. - 2. CNV simulation in whole exome. CNV-Sim wraps the functionality of Wessim to introduce variations in the targets. - - Homepage: http://nabavilab.github.io/CNV-Sim/ - ]]></help> -</tool> \ No newline at end of file |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/ART/README.md --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/ART/README.md Sat Aug 06 15:27:30 2016 -0400 |
| b |
| @@ -0,0 +1,4 @@ +# ART: a next-generation sequencing read simulator +The binary file `art_illumina` should be download and put here in this directory. + +Do so through the appropriate setup script in the home directory \ No newline at end of file |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/ART/art_illumina |
| b |
| Binary file cnvsim/ART/art_illumina has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/Prep_BlatSearch.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/Prep_BlatSearch.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,44 @@ +import sys +import os +import argparse + +MINIDENTITY = "90" +MINSCORE = "100" + +def main(argv): + parser = argparse.ArgumentParser(description='Blat Search for probe hybridization in Wessim2', prog='Prep_BlatSearch', formatter_class=argparse.RawTextHelpFormatter) + group1 = parser.add_argument_group('Mandatory input files') + group1.add_argument('-R', metavar = 'FILE', dest = 'reference', required=True, help = '2bit formatted reference file') + group1.add_argument('-P', metavar = 'FILE', dest = 'probe', required=True, help = 'FASTA format probe sequence file generated from Prep_Probe2Fa') + group2 = parser.add_argument_group('Search options') + group2.add_argument('-i', metavar = 'INT', dest = 'minidentity', required=False, help = 'Min-identity for blat match [90]', default="90") + group2.add_argument('-s', metavar = 'INT', dest = 'minscore', required=False, help = 'Min-Score for blat match [100]', default="100") + group3 = parser.add_argument_group('Output options') + group3.add_argument('-o', metavar = 'FILE', dest = 'outfile', required=False, help = 'Output file name [Probe_File_Name.psl]', default='') + + args = parser.parse_args() + ref = args.reference + probefile = args.probe + outfile = args.outfile + if outfile=='': + outfile = probefile + ".psl" + MINIDENTITY = args.minidentity + MINSCORE = args.minscore + + serverStopCommand = "gfServer stop localhost 6666" + serverStartCommand = "gfServer start -canStop localhost 6666 " + ref + blatCommand = "gfClient localhost 6666 / -minIdentity=" + MINIDENTITY + " -minScore=" + MINSCORE + " " + probefile + " " + outfile +# print serverStopCommand +# os.system(serverStopCommand) +# print serverStartCommand +# os.system(serverStartCommand) + print blatCommand + os.system(blatCommand) + + +def usage(): + print ">python Prep_BlatSearch.py ref.2bit probes.txt.fa output.psl" + sys.exit() + +if __name__=="__main__": + main(sys.argv[1:]) |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/Prep_Probe2Fa.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/Prep_Probe2Fa.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,31 @@ +import sys + +def main(argv): + if len(argv) < 1: + usage() + probefile = argv[0] + f = open(probefile) + w = open(probefile + ".fa", 'w') + line = f.readline() + line = f.readline() + while line: + values = line.split("\t") + if len(values) < 3: + line = f.readline() + continue + target = values[0] + seqid = values[1] + seq = values[2] + w.write(">" + seqid + "-" + target + "\n") + w.write(seq + "\n") + line = f.readline() + f.close() + w.close() + +def usage(): + print ">Python Prep_Probe2Fa.py probe.txt" + sys.exit() + + +if __name__=="__main__": + main(sys.argv[1:]) |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/README.md --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/README.md Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,108 @@ +### Introduction +**Wessim** is a simulator for a targeted resequencing as generally known as exome sequencing. Wessim basically generates a set of *artificial* DNA fragments for next generation sequencing (NGS) read simulation. In the targeted resequencing, we constraint the genomic regions that are used to generated DNA fragments to be only a part of the entire genome; they are usually exons and/or a few introns and untranslated regions (UTRs). + +### Install Wessim +Download Wessim using the links in this page, or go to https://github.com/sak042/Wessim +To run Wessim, Python 2.7 or later is required. To install Python, go to http://python.org/ + +### Requirements +The following programs are required to run Wessim or to prepare input files: +* **pysam** library: go to http://code.google.com/p/pysam/ to install pysam +* **numpy** library: go to http://numpy.scipy.org/ to install numpy +* **gfServer** and **gfClient**: In probe hybridization mode, Wessim runs more than 100,000 queries against the reference genome. This essentially requires a local blat server. gfServer and gfClient are pre-compiled programs for establishing private blat server on your computer. go to http://hgdownload.cse.ucsc.edu/admin/exe/ to download gfServer and gfClient (and set your local path to access the two programs anywhere). For more details about the tools, please refer to http://genome.ucsc.edu/FAQ/FAQblat.html#blat5 +* **faToTwoBit**: go to http://hgdownload.cse.ucsc.edu/admin/exe/ and download faToTwoBit. This is required to convert your FASTA file to .2bit +* **samtools**: samtools is needed to index your sample genome FASTA file (samtools faidx). +* **GemSim** error models: Wessim uses GemSim's empirical error models for NGS read generation. Go to GemSim's project page (http://sourceforge.net/projects/gemsim/) to download GemSim. You will find several model files (e.g. ill100v4_p.gzip) under 'models' directory. Save them and remember their location. + +### Preparing Input Files +Wessim requires two major inputs. One is the sample genome sequence, and the other is the target region information. +* **Sample genome sequence**: This is a FASTA file (e.g. ref.fa). You will need to index the file and generate .2bit +<pre><code> +>samtools faidx ref.fa +>faToTwoBit ref.fa ref.2bit +</code></pre> +* **Target region information**: Target regions can be specified by two different ways. + 1. **Ideal targets**: In ideal target mode, you will provide a list of genomic coordinates in a BED file (e.g. chr1 798833 799125). Ideal targets of major exome capture platforms are freely available from vendor's website. For Agilent's SureSelect platforms, go to https://earray.chem.agilent.com/suredesign/ . You must register at their site. After logging in, go to Find Designs and select Agilent Catalog at the menu tab. You will be able to download all information of currently available platforms including ideal target BED files and probe sequence text files. For NimbleGen's SeqCap go to http://www.nimblegen.com/products/seqcap/index.html and find BED files under Design and Annotation Files. + 2. **Probe sequences**: Probe sequences are available for SureSelect platforms in the SureDesign homepage (https://earray.chem.agilent.com/suredesign/) (see above). Usually those files are named "[platform]_probe.txt" + +### Running Wessim +There are two main scripts in the package - Wessim1.py and Wessim2.py. You will use Wessim1 if you are using a BED file for target regions (ideal target approach). However, it is highly recommended to use Wessim2 (probe hybridization approach) when the probe sequence is available; it is much more realistic and recovers the statistics of real data. Two other scripts that start with 'Prep' are used to preparing Wessim2 inputs. You can ignore remaining scripts that start with '__sub'; main Wessim programs will execute these sub scripts automatically. + +The basic synopsis of Wessim1 is like below: +<pre><code> +# Run Wessim1 in ideal target mode +>python Wessim1.py -R ref.fa -B target.bed -n 1000000 -l 100 -M model.gzip -z -o result -t 4 +</code></pre> +This will generate *result.fastq.gz* (single-end mode / gzip compressed) using 4 threads (CPU cores). + +For Wessim2: +<pre><code> +# Generate a FASTA file of probe sequence +>python Prep_Probe2Fa.py probe.txt (this generates probe.txt.fa) +# Establish your local blat server +>gfServer start -canStop localhost 6666 ref.2bit (this will consume one whole thread, you need to use a separated thread to continue the following steps) +# Run blat search to generate the match list +>python Prep_BlatSearch.py -R ref.2bit -P probe.txt.fa -o probe.txt.fa.psl (Note that the path to ref.2bit is not based on your local machine. It should be used without path, because the gfServer has it in its root) +# Run Wessim2 in probe hybridization mode. +>python Wessim2.py -R ref.fa -P probe.txt.fa -B probe_match.txt.fa.psl -n 1000000 -l 76 -M model.gzip -pz -o result +</code></pre> +This will generate *result_1.fastq.gz* and *result_2.fastq.gz* (paired-end mode / gzip compressed). + +### Running in metagenomic mode +You can use more than one genome as your template. To run Wessim in metagenomic mode, you can just write a simple description file that ends with (.meta). Use the meta description file at the place of your reference FASTA file (-R option) +<pre><code> +>python Wessim2.py -R reference.meta -P probe.txt.fa -B probe_match.txt.fa.psl -n 1000000 -l 76 -M model.gzip -pz -o result +</code></pre> + +The format of a tab-delimited .meta file is like below, +<pre><code> +genome1.fasta <tab> abundance of genome1 +genome2.fasta <tab> abundance of genome2 +... +genomeN.fasta <tab> abundance of genomeN +</code></pre> +For example, you can use, +<pre><code> +/data/reference/ref1.fa 0.2 +/data/reference/ref2.fa 0.4 +/data/reference/ref3.fa 0.4 +</code></pre> +Please make sure the overall abundances add up to 1 + +### Wessim Options +You can use '-h' for detailed help in command line. + +``` +Mandatory input files (for Wessim1 and Wessim2 in common): + -R FILE faidx-indexed (R)eference genome FASTA file +For Wessim1 only: + -B FILE Target region .(B)ED file +For Wessim2 only: + -P FILE (P)robe sequence FASTA file + -B FILE (B)lat matched probe regions .PSL file + +Parameters for exome capture: + -f INT mean (f)ragment size. this corresponds to insert size when sequencing in paired-end mode. [200] + -d INT standard (d)eviation of fragment size [50] + -m INT (m)inimum fragment length [read_length + 20] + -x INT slack margin of the given boundaries [0] (only for Wessim1) + -w INT penalty (w)eight for indel in the hybridization [2] (only for Wessim2) + +Parameters for sequencing: + -p generate paired-end reads [single] + -n INT total (n)umber of reads + -l INT read (l)ength (bp) + -M FILE GemSim (M)odel file (.gzip) + -t INT number of (t)hreaded subprocesses [1] + +Output options: + -o FILE (o)utput file header. ".fastq.gz" or ".fastq" will be attached automatically. Output will be splitted into two files in paired-end mode + -z compress output with g(z)ip [false] + -q INT (q)uality score offset [33] + -v (v)erbose; print out intermediate messages. +``` + +### Support or Contact +For GitHub use, check out the documentation at http://help.github.com/pages or contact support@github.com and well help you sort it out. + + |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/Wessim1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/Wessim1.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,197 @@ +import sys +import random +import bisect +import pysam +import gzip +import cPickle +import numpy +from time import time, localtime, strftime +import argparse +from multiprocessing import Process +import os +import math +import pysam + +inds={'A':0,'T':1,'G':2,'C':3,'N':4,'a':0,'t':1,'g':2,'c':3,'n':4} + +def subprogram(command, name): + os.system(command) + print "exiting subprocess " + str(name) + +def main(argv): + t0 = time() + arguline = " ".join(argv) + parser = argparse.ArgumentParser(description='Wessim1: Whole Exome Sequencing SIMulator 1 (Ideal target region-based version)', prog='Wessim1', formatter_class=argparse.RawTextHelpFormatter) + group1 = parser.add_argument_group('Mandatory input files') + group1.add_argument('-R', metavar = 'FILE', dest='reference', required=True, help='faidx-indexed (R)eference genome FASTA file') + group1.add_argument('-B', metavar = 'FILE', dest='region', required=True, help='Target region .(B)ED file') + + group2 = parser.add_argument_group('Parameters for exome capture') + group2.add_argument('-f', metavar = 'INT', type=int, dest='fragsize', required=False, help='mean (f)ragment size. this corresponds to insert size when sequencing in paired-end mode. [200]', default=200) + group2.add_argument('-d', metavar = 'INT', type=int, dest='fragsd', required=False, help='standard (d)eviation of fragment size [50]', default=50) + group2.add_argument('-m', metavar = 'INT', type=int, dest='fragmin', required=False, help='(m)inimum fragment length [read_length + 20]') + group2.add_argument('-x', metavar = 'INT',type=int, dest='slack', required=False, help='slack margin of the given boundaries [0]', default=0) + + group3 = parser.add_argument_group('Parameters for sequencing') + group3.add_argument('-p', action='store_true', help='generate paired-end reads [single]') + group3.add_argument('-n', metavar = 'INT', type=int, dest='readnumber', required=True, help='total (n)umber of reads') + group3.add_argument('-l', metavar = 'INT', type=int, dest='readlength', required=True, help='read (l)ength (bp)') + group3.add_argument('-M', metavar = 'FILE', dest='model', required=True, help='GemSim (M)odel file (.gzip)') + group3.add_argument('-t', metavar = 'INT', type=int, dest='threadnumber', required=False, help='number of (t)hreaded subprocesses [1]', default=1) + + group4 = parser.add_argument_group('Output options') + group4.add_argument('-o', metavar = 'FILE', dest='outfile', help='(o)utput file header. ".fastq.gz" or ".fastq" will be attached automatically. Output will be splitted into two files in paired-end mode', required=True) + group4.add_argument('-z', action='store_true', help='compress output with g(z)ip [false]') + group4.add_argument('-q', metavar = 'INT', type=int, dest='qualbase', required=False, help='(q)uality score offset [33]', default=33) + group4.add_argument('-v', action='store_true', help='(v)erbose; print out intermediate messages.') + + args = parser.parse_args() + reffile = args.reference + regionfile = args.region + + isize = args.fragsize + isd = args.fragsd + imin = args.fragmin + slack = args.slack + + getRegionVector(reffile, regionfile, slack) + paired = args.p + readlength = args.readlength + readnumber = args.readnumber + threadnumber = args.threadnumber + + if imin==None: + if paired: + imin = readlength + 20 + else: + imin = readlength + 20 + if isize < imin: + print "too small mean fragment size (" + str(isize) + ") compared to minimum length (" + str(imin) + "). Increase it and try again." + sys.exit(0) + model = args.model + + outfile = args.outfile + compress = args.z + qualbase = args.qualbase + verbose = args.v + + print "-------------------------------------------" + print "Reference:", reffile + print "Region file:", regionfile + print "Fragment:",isize, "+-", isd, ">", imin + print "Paired-end mode?", paired + print "Sequencing model:", model + print "Read length:", readlength, "Read number:", readnumber + print "Output File:", outfile + print "Gzip compress?", compress + print "Quality base:", qualbase + print "Thread number:", threadnumber + print "Job started at:", strftime("%Y-%m-%d %H:%M:%S", localtime()) + print "-------------------------------------------" + + processes = [] + for t in range(0, threadnumber): + readstart = int(float(readnumber) / float(threadnumber) * t) + 1 + readend = int(float(readnumber) / float(threadnumber) * (t+1)) + command = "python __sub_wessim1.py " + arguline + " -1 " + str(readstart) + " -2 " + str(readend) + " -i " + str(t+1) + p = Process(target=subprogram, args=(command, t+1)) + p.start() + processes.append(p) + for p in processes: + p.join() + t1 = time() + print "Done generating " + str(readnumber) + " reads in %f secs" % (t1 - t0) + print "Merging subresults..." + wread = None + wread2 = None + if paired and compress: + wread = gzip.open(outfile + "_1.fastq.gz", 'wb') + wread2 = gzip.open(outfile + "_2.fastq.gz", 'wb') + elif paired and not compress: + wread = open(outfile + "_1.fastq", 'w') + wread2 = open(outfile + "_2.fastq", 'w') + elif not paired and compress: + wread = gzip.open(outfile + ".fastq.gz", 'wb') + else: + wread = open(outfile + ".fastq", 'w') + if not paired: + for t in range(0, threadnumber): + suboutfile = outfile + "-" + str(t+1) + fread = None + if compress: + suboutfile += ".fastq.gz" + fread = gzip.open(suboutfile, 'rb') + else: + suboutfile += ".fastq" + fread = open(suboutfile, 'r') + line = fread.readline() + while line: + wread.write(line) + line = fread.readline() + fread.close() + os.remove(suboutfile) + wread.close() + else: + for t in range(0, threadnumber): + suboutfile1 = outfile + "-" + str(t+1) + "_1" + suboutfile2 = outfile + "-" + str(t+1) + "_2" + fread1 = None + fread2 = None + if compress: + suboutfile1 += ".fastq.gz" + suboutfile2 += ".fastq.gz" + fread1 = gzip.open(suboutfile1, "rb") + fread2 = gzip.open(suboutfile2, "rb") + else: + suboutfile1 += ".fastq" + suboutfile2 += ".fastq" + fread1 = open(suboutfile1, "r") + fread2 = open(suboutfile2, "r") + line1 = fread1.readline() + line2 = fread2.readline() + while line1 and line2: + wread.write(line1) + wread2.write(line2) + line1 = fread1.readline() + line2 = fread2.readline() + fread1.close() + fread2.close() + os.remove(suboutfile1) + os.remove(suboutfile2) + wread.close() + wread2.close() + sys.exit(0) + +def getRegionVector(fastafile, regionfile, slack): + print "Generating fasta file for given regions..." + faoutfile = regionfile + ".fa" + abdoutfile = regionfile + ".abd" + ref = pysam.Fastafile(fastafile) + f = open(regionfile) + wfa = open(faoutfile, 'w') + wabd = open(abdoutfile, 'w') + i = f.readline() + abd = 0 + while i: + i = f.readline() + values = i.split("\t") + if i.startswith("#") or len(values)<3: + continue + chrom = values[0] + start = max(int(values[1]) - slack, 1) + end = int(values[2]) + slack + header = ">" + chrom + "_" + str(start) + "_" + str(end) + x = ref.fetch(chrom, start, end) + length = len(x) + abd += length + wfa.write(header + "\n") + wfa.write(x + "\n") + wabd.write(str(abd) + "\n") + f.close() + wfa.close() + wabd.close() + +if __name__=="__main__": + main(sys.argv[1:]) |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/Wessim2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/Wessim2.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,170 @@ +import sys +import random +import bisect +import pysam +import gzip +import cPickle +import numpy +from time import time, localtime, strftime +import argparse +from multiprocessing import Process +import os +import math + +inds={'A':0,'T':1,'G':2,'C':3,'N':4,'a':0,'t':1,'g':2,'c':3,'n':4} + +def subprogram(command, name): + os.system(command) + print "exiting subprocess " + str(name) + +def main(argv): + t0 = time() + arguline = " ".join(argv) + parser = argparse.ArgumentParser(description='Wessim2: Whole Exome Sequencing SIMulator 2 (Probe-based version)', prog='Wessim2', formatter_class=argparse.RawTextHelpFormatter) + group1 = parser.add_argument_group('Mandatory input files') + group1.add_argument('-R', metavar = 'FILE', dest='reference', required=True, help='faidx-indexed (R)eference genome FASTA file or meta description file (.meta)') + group1.add_argument('-P', metavar = 'FILE', dest='probe', required=True, help='(P)robe sequence FASTA file') + group1.add_argument('-B', metavar = 'FILE', dest='probeblat', required=True, help='(B)lat matched probe regions .PSL file') + + group2 = parser.add_argument_group('Parameters for exome capture') + group2.add_argument('-f', metavar = 'INT', type=int, dest='fragsize', required=False, help='mean (f)ragment size. this corresponds to insert size when sequencing in paired-end mode. [200]', default=200) + group2.add_argument('-d', metavar = 'INT', type=int, dest='fragsd', required=False, help='standard (d)eviation of fragment size [50]', default=50) + group2.add_argument('-m', metavar = 'INT', type=int, dest='fragmin', required=False, help='(m)inimum fragment length [read_length + 20]') + group2.add_argument('-y', metavar = 'PERCENT',type=int, dest='bind', required=False, help='minimum required fraction of probe match to be h(y)bridized [50]', default=50) + group2.add_argument('-w', metavar = 'INT', type=int, dest='weight', required=False, help='penalty (w)eight for indel in the hybridization [2]', default=2) + + group3 = parser.add_argument_group('Parameters for sequencing') + group3.add_argument('-p', action='store_true', help='generate paired-end reads [single]') + group3.add_argument('-n', metavar = 'INT', type=int, dest='readnumber', required=True, help='total (n)umber of reads') + group3.add_argument('-l', metavar = 'INT', type=int, dest='readlength', required=True, help='read (l)ength (bp)') + group3.add_argument('-M', metavar = 'FILE', dest='model', required=True, help='GemSim (M)odel file (.gzip)') + group3.add_argument('-t', metavar = 'INT', type=int, dest='threadnumber', required=False, help='number of (t)hreaded subprocesses [1]', default=1) + + group4 = parser.add_argument_group('Output options') + group4.add_argument('-o', metavar = 'FILE', dest='outfile', help='(o)utput file header. ".fastq.gz" or ".fastq" will be attached automatically. Output will be splitted into two files in paired-end mode', required=True) + group4.add_argument('-z', action='store_true', help='compress output with g(z)ip [false]') + group4.add_argument('-q', metavar = 'INT', type=int, dest='qualbase', required=False, help='(q)uality score offset [33]', default=33) + group4.add_argument('-v', action='store_true', help='(v)erbose; print out intermediate messages.') + + args = parser.parse_args() + reffile = args.reference + probefile = args.probe + alignfile = args.probeblat + + isize = args.fragsize + isd = args.fragsd + imin = args.fragmin + bind = args.bind + + paired = args.p + readlength = args.readlength + readnumber = args.readnumber + threadnumber = args.threadnumber + if imin==None: + if paired: + imin = readlength + 20 + else: + imin = readlength + 20 + if isize < imin: + print "too small mean fragment size (" + str(isize) + ") compared to minimum length (" + str(imin) + "). Increase it and try again." + sys.exit(0) + model = args.model + + outfile = args.outfile + compress = args.z + qualbase = args.qualbase + verbose = args.v + + print "-------------------------------------------" + print "Reference:", reffile + print "Probeset:", probefile + print "Probematch:", alignfile + print "Fragment:",isize, "+-", isd, ">", imin + print "Paired-end mode?", paired + print "Sequencing model:", model + print "Read length:", readlength, "Read number:", readnumber + print "Output File:", outfile + print "Gzip compress?", compress + print "Quality base:", qualbase + print "Thread number:", threadnumber + print "Job started at:", strftime("%Y-%m-%d %H:%M:%S", localtime()) + print "-------------------------------------------" + + processes = [] + for t in range(0, threadnumber): + readstart = int(float(readnumber) / float(threadnumber) * t) + 1 + readend = int(float(readnumber) / float(threadnumber) * (t+1)) + command = "python __sub_wessim2.py " + arguline + " -1 " + str(readstart) + " -2 " + str(readend) + " -i " + str(t+1) + p = Process(target=subprogram, args=(command, t+1)) + p.start() + processes.append(p) + for p in processes: + p.join() + t1 = time() + print "Done generating " + str(readnumber) + " reads in %f secs" % (t1 - t0) + print "Merging subresults..." + wread = None + wread2 = None + if paired and compress: + wread = gzip.open(outfile + "_1.fastq.gz", 'wb') + wread2 = gzip.open(outfile + "_2.fastq.gz", 'wb') + elif paired and not compress: + wread = open(outfile + "_1.fastq", 'w') + wread2 = open(outfile + "_2.fastq", 'w') + elif not paired and compress: + wread = gzip.open(outfile + ".fastq.gz", 'wb') + else: + wread = open(outfile + ".fastq", 'w') + if not paired: + for t in range(0, threadnumber): + suboutfile = outfile + "-" + str(t+1) + fread = None + if compress: + suboutfile += ".fastq.gz" + fread = gzip.open(suboutfile, 'rb') + else: + suboutfile += ".fastq" + fread = open(suboutfile, 'r') + line = fread.readline() + while line: + wread.write(line) + line = fread.readline() + fread.close() + os.remove(suboutfile) + wread.close() + else: + for t in range(0, threadnumber): + suboutfile1 = outfile + "-" + str(t+1) + "_1" + suboutfile2 = outfile + "-" + str(t+1) + "_2" + fread1 = None + fread2 = None + if compress: + suboutfile1 += ".fastq.gz" + suboutfile2 += ".fastq.gz" + fread1 = gzip.open(suboutfile1, "rb") + fread2 = gzip.open(suboutfile2, "rb") + else: + suboutfile1 += ".fastq" + suboutfile2 += ".fastq" + fread1 = open(suboutfile1, "r") + fread2 = open(suboutfile2, "r") + line1 = fread1.readline() + line2 = fread2.readline() + while line1 and line2: + wread.write(line1) + wread2.write(line2) + line1 = fread1.readline() + line2 = fread2.readline() + fread1.close() + fread2.close() + os.remove(suboutfile1) + os.remove(suboutfile2) + wread.close() + wread2.close() + sys.exit(0) + + +if __name__=="__main__": + main(sys.argv[1:]) |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/__sub_wessim1.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/__sub_wessim1.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| b'@@ -0,0 +1,819 @@\n+import sys\n+import random\n+import bisect\n+import pysam\n+import gzip\n+import cPickle\n+import numpy\n+from time import time\n+import argparse\n+import math\n+\n+inds={\'A\':0,\'T\':1,\'G\':2,\'C\':3,\'N\':4,\'a\':0,\'t\':1,\'g\':2,\'c\':3,\'n\':4}\n+\n+def main(argv):\n+\tt0 = time()\n+\tparser = argparse.ArgumentParser(description=\'sub-wessim: a sub-program for Wessim1. (NOTE!) Do not run this program. Use "Wessim1.py" instead. \', prog=\'wessim1_sub\', formatter_class=argparse.RawTextHelpFormatter)\n+\tgroup1 = parser.add_argument_group(\'Mandatory input files\')\n+\tgroup1.add_argument(\'-R\', metavar = \'FILE\', dest=\'reference\', required=True, help=\'(R)eference genome FASTA file\')\n+\tgroup1.add_argument(\'-B\', metavar = \'FILE\', dest=\'region\', required=True, help=\'Target region .(B)ED file\')\n+\n+\tgroup2 = parser.add_argument_group(\'Parameters for exome capture\')\n+\tgroup2.add_argument(\'-f\', metavar = \'INT\', type=int, dest=\'fragsize\', required=False, help=\'mean (f)ragment size. this corresponds to insert size when sequencing in paired-end mode. [200]\', default=200)\n+\tgroup2.add_argument(\'-d\', metavar = \'INT\', type=int, dest=\'fragsd\', required=False, help=\'standard (d)eviation of fragment size [50]\', default=50)\n+\tgroup2.add_argument(\'-m\', metavar = \'INT\', type=int, dest=\'fragmin\', required=False, help=\'(m)inimum fragment length [read_length + 20 for single-end, 2*read_length + 20 for paired-end]\')\n+\tgroup2.add_argument(\'-y\', metavar = \'PERCENT\',type=int, dest=\'bind\', required=False, help=\'minimum required fraction of probe match to be h(y)bridized [50]\', default=50) \n+\t\n+\tgroup3 = parser.add_argument_group(\'Parameters for sequencing\')\n+\tgroup3.add_argument(\'-p\', action=\'store_true\', help=\'generate paired-end reads [single]\')\n+\tgroup3.add_argument(\'-n\', help=\'do not care\')\n+\tgroup3.add_argument(\'-1\', metavar = \'INT\', type=int, dest=\'readstart\', required=True, help=\'start number of read\')\n+\tgroup3.add_argument(\'-2\', metavar = \'INT\', type=int, dest=\'readend\', required=True, help=\'end number of read\')\n+\tgroup3.add_argument(\'-l\', metavar = \'INT\', type=int, dest=\'readlength\', required=True, help=\'read (l)ength (bp)\')\n+\tgroup3.add_argument(\'-i\', metavar = \'INT\', type=int, dest=\'processid\', required=True, help=\'subprocess (i)d\')\n+\tgroup3.add_argument(\'-M\', metavar = \'FILE\', dest=\'model\', required=True, help=\'GemSim (M)odel file (.gzip)\')\n+\tgroup3.add_argument(\'-t\', help=\'do not care\')\n+\n+\tgroup4 = parser.add_argument_group(\'Output options\')\n+\tgroup4.add_argument(\'-o\', metavar = \'FILE\', dest=\'outfile\', help=\'(o)utput file header. ".fastq.gz" or ".fastq" will be attached automatically. Output will be splitted into two files in paired-end mode\', required=True)\n+\tgroup4.add_argument(\'-z\', action=\'store_true\', help=\'compress output with g(z)ip [false]\')\n+\tgroup4.add_argument(\'-q\', metavar = \'INT\', type=int, dest=\'qualbase\', required=False, help=\'(q)uality score offset [33]\', default=33)\n+\tgroup4.add_argument(\'-v\', action=\'store_true\', help=\'(v)erbose; print out intermediate messages.\')\n+\n+\targs = parser.parse_args()\n+\treffile = args.reference\n+\tregionfile = args.region\n+\tfaoutfile = regionfile + ".fa"\n+\tabdoutfile = regionfile + ".abd"\n+\t\n+\tisize = args.fragsize\n+\tisd = args.fragsd\n+\timin = args.fragmin\n+\tbind = args.bind\n+\tsubid = args.processid\n+\n+\tpaired = args.p\n+\treadlength = args.readlength\n+\treadstart = args.readstart\n+\treadend = args.readend\t\t\n+\tif imin==None:\n+\t\tif paired:\n+\t\t\timin = readlength + 20\n+\t\telse:\n+\t\t\timin = readlength + 20\n+\tif isize < imin:\n+\t\tprint "too small mean fragment size (" + str(isize) + ") compared to minimum length (" + str(imin) + "). Increase it and try again."\n+\t\tsys.exit(0) \n+\tmodel = args.model\n+\n+\tf = open(faoutfile)\n+\ti = f.readline()\n+\tseqlist = []\n+\tabdlist = []\n+\twhile i:\n+\t\theader = i.strip()[1:]\n+\t\tseq = f.readline().strip()\n+\t\tseqlist.append((header, seq))\n+\t\ti = f.readline()\n+\tf.close()\n+\tf = open(abdoutfile)\n+\ti = f.readline()\n+\twhile i:\n+\t\tabd = int(i.strip())\n+\t\tabdlist.append(abd)\n+\t\ti = f.readline()\n+\tf.close()\n+\tl'..b' bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>g:\n+\t\t\tread=read[:pos+3]+\'C\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>t:\n+\t\t\tread=read[:pos+3]+\'G\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>a:\n+\t\t\tread=read[:pos+3]+\'T\'+read[pos+4:] \n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telse:\n+\t\t\tread=read[:pos+3]+\'A\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\tif index in delD:\n+\t\t\tdelete=delD[index]()\n+\t\t\tread=read[:pos+4]+read[pos+delete+4:]\n+\t\tif index in insD:\n+\t\t\tinsert=insD[index]()\n+\t\t\tread=read[:pos+4]+insert+read[pos+4:]\n+\t\t\tfor i in insert:\n+\t\t\t\tiPos=pos-1\n+\t\t\t\twhile iPos>=0:\n+\t\t\t\t\ttry:\n+\t\t\t\t\t\tqualslist.append(iQ[iPos]())\n+\t\t\t\t\t\tsuccess=True\n+\t\t\t\t\t\tbreak\n+\t\t\t\t\texcept:\n+\t\t\t\t\t\tiPos-=1\n+\t\t\t\t\tif success==False:\n+\t\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\t\tpos+=len(insert)\n+\t\tpos+=1\n+\tqualslist.append(qualslist[-1])\n+\treadback = read\n+\tread=read[4:readLen+4]\n+\tquals=\'\'.join(qualslist)[:readLen]\n+\tif len(quals)!=len(read):\n+\t\tprint "unexpected stop"\n+\t\treturn None, None\n+\treturn read,quals\t \n+\n+def generateM(sd, newSD, x,t, gcVector):\n+\tgcSD = numpy.std(gcVector)*(newSD/sd)\n+\ts00 = gcSD*gcSD + newSD*newSD*t*t\n+\ts11 = newSD*newSD\n+\trho = newSD*t/math.sqrt(s00)\n+\tm = numpy.matrix([[s00, rho*math.sqrt(s00*s11)], [rho*math.sqrt(s00*s11), s11]])\n+\tw, v = numpy.linalg.eig(m)\n+\td = numpy.matrix([[math.sqrt(w[0]),0],[0,math.sqrt(w[1])]])\n+\tM = v*d\n+\treturn M, m\n+\n+def generateMatrices(sd,x, gcVector):\n+\tM1, m1 = generateM(sd, sd, x,1/0.9, gcVector)\n+\te1 = numpy.matrix([[1],[0]])\n+\te2 = numpy.matrix([[0],[1]])\n+\tlongAxis1 = M1*e1\n+\tlongAxis2 = M1*e2\n+\tlongAxis = longAxis1\n+\tif norm(longAxis1) < norm(longAxis2):\n+\t\tlongAxis = longAxis2\t\t\n+\tM2 = []\n+\tm2 = []\n+\tnewSD = sd;\n+\tfor i in range(100, 1000):\n+\t\tnewSD = sd*i/100.0\n+\t\tM2, m2= generateM(sd, newSD,x,0.5, gcVector)\n+\t\tif norm(numpy.linalg.inv(M2)*longAxis)<1.0:\n+\t\t\tbreak\n+\tu1 = numpy.linalg.inv(M1)\n+\tu2 = numpy.linalg.inv(M2)\n+\treturn u1, u2, newSD, m1, m2\n+\n+def getProb(l,n,x,sd,gcSD,alpha, mvnpdf):\n+\tp1 = mvnpdf[0][int(cut((l-x)/sd)*100)]\n+\tp2 = mvnpdf[0][int(cut((n-(x/2+(l-x)*alpha))/(l*gcSD/x))*100)]\n+\treturn float(p1)*float(p2)\n+\t\n+\t\n+def H2(l, n, x, sd1, sd2, gcSD, mvnpdf):\n+\tbp = getProb(l,n,x,sd1,gcSD,.5,mvnpdf)\n+\tap = getProb(l,n,x,sd2,gcSD,9/7,mvnpdf)\n+\tv = ap/bp\n+\n+\tr = random.random()\n+\ttoKeep = v > r\n+\treturn toKeep\t\n+\t\n+\t\n+def norm(x):\n+\ty=x[0]*x[0]+x[1]*x[1]\n+\treturn math.sqrt(y)\n+\t\n+def cut(x):\n+\ty = abs(x)\n+\tif y >5.00:\n+\t\ty = 5.00\n+\treturn y\n+\n+def H(l, n, x, u1, u2, mvnpdf):\n+\tu = numpy.matrix([[x/2], [x]])\n+\tnl1 = numpy.matrix([[n],[l]])\n+\tv1 = u1*(nl1-u)\n+\tv2 = u2*(nl1-u)\n+\n+\tp1 = mvnpdf[int(cut(v1[0])*100)][int(cut(v1[1])*100)]\n+\tp2 = mvnpdf[int(cut(v2[0])*100)][int(cut(v2[1])*100)]\n+\tv = float(p1)/float(p2)\n+\t\n+\tr = random.random()\n+\ttoKeep = v > r\n+\treturn toKeep\n+\t\n+def readmvnTable():\n+\tf = open("lib/mvnTable.txt")\n+\tcontext = f.read()\n+\tlines = context.split("\\n")\n+\tmvnTable = []\n+\tfor line in lines:\n+\t\tvalues = line.split("\\t")\n+\t\tif len(values)<500:\n+\t\t\tcontinue\n+\t\tmvnTable.append(values)\n+\tf.close()\n+\treturn mvnTable\n+\n+def getIndex(abdlist, pos):\n+\ti = bisect.bisect_right(abdlist, pos)\n+\treturn i \n+\n+if __name__=="__main__":\n+\tmain(sys.argv[1:])\n+\tsys.exit(0)\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/__sub_wessim2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/__sub_wessim2.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| b'@@ -0,0 +1,900 @@\n+import sys\n+import random\n+import bisect\n+import pysam\n+import gzip\n+import cPickle\n+import numpy\n+from time import time\n+import argparse\n+import math\n+\n+inds={\'A\':0,\'T\':1,\'G\':2,\'C\':3,\'N\':4,\'a\':0,\'t\':1,\'g\':2,\'c\':3,\'n\':4}\n+\n+def main(argv):\n+\tt0 = time()\n+\tparser = argparse.ArgumentParser(description=\'sub-wessim: a sub-program for Wessim2. (NOTE!) Do not run this program. Use "Wessim2.py" instead. \', prog=\'wessim2-sub\', formatter_class=argparse.RawTextHelpFormatter)\n+\tgroup1 = parser.add_argument_group(\'Mandatory input files\')\n+\tgroup1.add_argument(\'-R\', metavar = \'FILE\', dest=\'reference\', required=True, help=\'(R)eference genome FASTA file\')\n+\tgroup1.add_argument(\'-P\', metavar = \'FILE\', dest=\'probe\', required=True, help=\'(P)robe sequence FASTA file\')\n+\tgroup1.add_argument(\'-B\', metavar = \'FILE\', dest=\'probeblat\', required=True, help=\'(B)lat matched probe regions .PSL file\')\n+\n+\tgroup2 = parser.add_argument_group(\'Parameters for exome capture\')\n+\tgroup2.add_argument(\'-f\', metavar = \'INT\', type=int, dest=\'fragsize\', required=False, help=\'mean (f)ragment size. this corresponds to insert size when sequencing in paired-end mode. [200]\', default=200)\n+\tgroup2.add_argument(\'-d\', metavar = \'INT\', type=int, dest=\'fragsd\', required=False, help=\'standard (d)eviation of fragment size [50]\', default=50)\n+\tgroup2.add_argument(\'-m\', metavar = \'INT\', type=int, dest=\'fragmin\', required=False, help=\'(m)inimum fragment length [read_length + 20 for single-end, 2*read_length + 20 for paired-end]\')\n+\tgroup2.add_argument(\'-y\', metavar = \'PERCENT\',type=int, dest=\'bind\', required=False, help=\'minimum required fraction of probe match to be h(y)bridized [50]\', default=50) \n+\tgroup2.add_argument(\'-w\', metavar = \'INT\', type=int, dest=\'weight\', required=False, help=\'penalty (w)eight for indel in the hybridization [2]\', default=2)\n+\t\n+\tgroup3 = parser.add_argument_group(\'Parameters for sequencing\')\n+\tgroup3.add_argument(\'-p\', action=\'store_true\', help=\'generate paired-end reads [single]\')\n+\tgroup3.add_argument(\'-n\', help=\'do not care\')\n+\tgroup3.add_argument(\'-1\', metavar = \'INT\', type=int, dest=\'readstart\', required=True, help=\'start number of read (given by main process)\')\n+\tgroup3.add_argument(\'-2\', metavar = \'INT\', type=int, dest=\'readend\', required=True, help=\'end number of read (given by main process)\')\n+\tgroup3.add_argument(\'-l\', metavar = \'INT\', type=int, dest=\'readlength\', required=True, help=\'read (l)ength (bp)\')\n+\tgroup3.add_argument(\'-i\', metavar = \'INT\', type=int, dest=\'processid\', required=True, help=\'subprocess (i)d\')\n+\tgroup3.add_argument(\'-M\', metavar = \'FILE\', dest=\'model\', required=True, help=\'GemSim (M)odel file (.gzip)\')\n+\tgroup3.add_argument(\'-t\', help=\'do not care\')\n+\n+\tgroup4 = parser.add_argument_group(\'Output options\')\n+\tgroup4.add_argument(\'-o\', metavar = \'FILE\', dest=\'outfile\', help=\'(o)utput file header. ".fastq.gz" or ".fastq" will be attached automatically. Output will be splitted into two files in paired-end mode\', required=True)\n+\tgroup4.add_argument(\'-z\', action=\'store_true\', help=\'compress output with g(z)ip [false]\')\n+\tgroup4.add_argument(\'-q\', metavar = \'INT\', type=int, dest=\'qualbase\', required=False, help=\'(q)uality score offset [33]\', default=33)\n+\tgroup4.add_argument(\'-v\', action=\'store_true\', help=\'(v)erbose; print out intermediate messages.\')\n+\n+\targs = parser.parse_args()\n+\treffile = args.reference\n+\tfref = None\n+\tfrefs = []\n+\tmetap = []\n+\tmetamode = False\n+\tprobefile = args.probe\n+\talignfile = args.probeblat\n+\t\n+\tisize = args.fragsize\n+\tisd = args.fragsd\n+\timin = args.fragmin\n+\tbind = args.bind\n+\tweight = args.weight\n+\tsubid = args.processid\n+\n+\tpaired = args.p\n+\treadlength = args.readlength\n+\treadstart = args.readstart\n+\treadend = args.readend\t\t\n+\tif imin==None:\n+\t\tif paired:\n+\t\t\timin = readlength + 20\n+\t\telse:\n+\t\t\timin = readlength + 20\n+\tif isize < imin:\n+\t\tprint "too small mean fragment size (" + str(isize) + ") compared to minimum length (" + str(imin) + "). Increase it and try again."\n+\t\ts'..b'))\n+\t\telif val>c:\n+\t\t\tread=read[:pos+3]+\'N\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>g:\n+\t\t\tread=read[:pos+3]+\'C\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>t:\n+\t\t\tread=read[:pos+3]+\'G\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telif val>a:\n+\t\t\tread=read[:pos+3]+\'T\'+read[pos+4:] \n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\telse:\n+\t\t\tread=read[:pos+3]+\'A\'+read[pos+4:]\n+\t\t\tbPos=pos-1\n+\t\t\twhile bPos>=0:\n+\t\t\t\ttry:\n+\t\t\t\t\tqualslist.append(bQ[bPos]())\n+\t\t\t\t\tsuccess=True\n+\t\t\t\t\tbreak\n+\t\t\t\texcept:\n+\t\t\t\t\tbPos-1\n+\t\t\t\tif success==False:\n+\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\tif index in delD:\n+\t\t\tdelete=delD[index]()\n+\t\t\tread=read[:pos+4]+read[pos+delete+4:]\n+\t\tif index in insD:\n+\t\t\tinsert=insD[index]()\n+\t\t\tread=read[:pos+4]+insert+read[pos+4:]\n+\t\t\tfor i in insert:\n+\t\t\t\tiPos=pos-1\n+\t\t\t\twhile iPos>=0:\n+\t\t\t\t\ttry:\n+\t\t\t\t\t\tqualslist.append(iQ[iPos]())\n+\t\t\t\t\t\tsuccess=True\n+\t\t\t\t\t\tbreak\n+\t\t\t\t\texcept:\n+\t\t\t\t\t\tiPos-=1\n+\t\t\t\t\tif success==False:\n+\t\t\t\t\t\tqualslist.append(chr(2+qual))\n+\t\t\tpos+=len(insert)\n+\t\tpos+=1\n+\tqualslist.append(qualslist[-1])\n+\treadback = read\n+\tread=read[4:readLen+4]\n+\tquals=\'\'.join(qualslist)[:readLen]\n+\tif len(quals)!=len(read):\n+\t\tprint "unexpected stop"\n+\t\treturn None, None\n+\treturn read,quals\t \n+\n+def generateM(sd, newSD, x,t, gcVector):\n+\tgcSD = numpy.std(gcVector)*(newSD/sd)\n+\ts00 = gcSD*gcSD + newSD*newSD*t*t\n+\ts11 = newSD*newSD\n+\trho = newSD*t/math.sqrt(s00)\n+\tm = numpy.matrix([[s00, rho*math.sqrt(s00*s11)], [rho*math.sqrt(s00*s11), s11]])\n+\tw, v = numpy.linalg.eig(m)\n+\td = numpy.matrix([[math.sqrt(w[0]),0],[0,math.sqrt(w[1])]])\n+\tM = v*d\n+\treturn M, m\n+\n+def generateMatrices(sd,x, gcVector):\n+\tM1, m1 = generateM(sd, sd, x,1/0.9, gcVector)\n+\te1 = numpy.matrix([[1],[0]])\n+\te2 = numpy.matrix([[0],[1]])\n+\tlongAxis1 = M1*e1\n+\tlongAxis2 = M1*e2\n+\tlongAxis = longAxis1\n+\tif norm(longAxis1) < norm(longAxis2):\n+\t\tlongAxis = longAxis2\t\t\n+\tM2 = []\n+\tm2 = []\n+\tnewSD = sd;\n+\tfor i in range(100, 1000):\n+\t\tnewSD = sd*i/100.0\n+\t\tM2, m2= generateM(sd, newSD,x,0.5, gcVector)\n+\t\tif norm(numpy.linalg.inv(M2)*longAxis)<1.0:\n+\t\t\tbreak\n+\tu1 = numpy.linalg.inv(M1)\n+\tu2 = numpy.linalg.inv(M2)\n+\treturn u1, u2, newSD, m1, m2\n+\n+def getProb(l,n,x,sd,gcSD,alpha, mvnpdf):\n+\tp1 = mvnpdf[0][int(cut((l-x)/sd)*100)]\n+\tp2 = mvnpdf[0][int(cut((n-(x/2+(l-x)*alpha))/(l*gcSD/x))*100)]\n+\treturn float(p1)*float(p2)\n+\t\n+\t\n+def H2(l, n, x, sd1, sd2, gcSD, mvnpdf):\n+\tbp = getProb(l,n,x,sd1,gcSD,.5,mvnpdf)\n+\tap = getProb(l,n,x,sd2,gcSD,9/7,mvnpdf)\n+\tv = ap/bp\n+\n+\tr = random.random()\n+\ttoKeep = v > r\n+\treturn toKeep\t\n+\t\n+\t\n+def norm(x):\n+\ty=x[0]*x[0]+x[1]*x[1]\n+\treturn math.sqrt(y)\n+\t\n+def cut(x):\n+\ty = abs(x)\n+\tif y >5.00:\n+\t\ty = 5.00\n+\treturn y\n+\n+def H(l, n, x, u1, u2, mvnpdf):\n+\tu = numpy.matrix([[x/2], [x]])\n+\tnl1 = numpy.matrix([[n],[l]])\n+\tv1 = u1*(nl1-u)\n+\tv2 = u2*(nl1-u)\n+\n+\tp1 = mvnpdf[int(cut(v1[0])*100)][int(cut(v1[1])*100)]\n+\tp2 = mvnpdf[int(cut(v2[0])*100)][int(cut(v2[1])*100)]\n+\tv = float(p1)/float(p2)\n+\t\n+\tr = random.random()\n+\ttoKeep = v > r\n+\treturn toKeep \n+\t\n+def readmvnTable():\n+\tf = open("lib/mvnTable.txt")\n+\tcontext = f.read()\n+\tlines = context.split("\\n")\n+\tmvnTable = []\n+\tfor line in lines:\n+\t\tvalues = line.split("\\t")\n+\t\tif len(values)<500:\n+\t\t\tcontinue\n+\t\tmvnTable.append(values)\n+\tf.close()\n+\treturn mvnTable\n+\t\n+if __name__=="__main__":\n+\tmain(sys.argv[1:])\n+\tsys.exit(0)\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/lib/S0293689_Probes.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/lib/S0293689_Probes.txt Sat Aug 06 15:27:30 2016 -0400 |
| b |
| b'@@ -0,0 +1,508726 @@\n+TargetID\tProbeID\tSequence\tReplication\tStrand\tCoordinates\n+mRNA|AL390972\tA_36_B233385\tGGACTTGCTTCCCCGTCTGACTGTCCTTGACTTCTAGAATGGAAGAAGCTGAGCTGGTGAAGGGAAGACTCCAGGCCATCACAGTAAGTCTGCATACAGTTATAACTCATTAACGTTGTT\tNA\t+\tchr1:100111836-100111955\n+mRNA|AL390972\tA_36_B244963\tCTGGAACTTTTTTATTTCAGGATAAAAGAAAAATACAGGAAGAAATCTCACAGAAGCGTCTGAAAATAGAGGAAGACAAACTAAAGCACCAGCATTTGAAGGTAATTTATGGCTTTAATT\tNA\t+\tchr1:100127855-100127974\n+mRNA|AL390972\tA_36_B857996\tATTGAAATTCATAAGATTTGAAAGCAGTAAGAGATTTTGGAAATTATTGACCAATATAAGAACATATCTTTATGTTTCCTGACATCAGAAAAAGGCCTTGAGGGAGAAATGGCTTCTAGA\tNA\t+\tchr1:100133110-100133229\n+mRNA|AL390972\tA_36_B260148\tGAGATTTTGGAAATTATTGACCAATATAAGAACATATCTTTATGTTTCCTGACATCAGAAAAAGGCCTTGAGGGAGAAATGGCTTCTAGATGGAATCAGCAGCGGAAAAGAACAGGAAGA\tNA\t+\tchr1:100133140-100133259\n+mRNA|AL390972\tA_36_B260149\tGATGAAGAAGCAAAATCAACAAGACCAGCACCAGATCCAGGTTCTAGAACAAAGTATCCTCAGGTATGGCCCTCACTGAGATACTCCACAACTACCTTTATTTTGACTTTATCCCATGTT\tNA\t+\tchr1:100133260-100133379\n+mRNA|AL390972\tA_36_B232404\tCAGGCTTGAGAAAGAGATCCAAGATCTTGAAAAAGCTGAACTGCAAATCTCAACGAAGGAAGAGGCCATTTTAAAGAAACTAAAGTCAATTGAGCGGACAACAGAAGACATTATAAGAGT\tNA\t+\tchr1:100152229-100152348\n+mRNA|AL390972\tA_36_B242879\tTTTTTAAACTGTATTAATCATGGAGAATTCTTTTTTCTTACAGTCTGTGAAAGTGGAAAGAGAAGAAAGAGCAGAAGGTATGTTTTTGAATAACTTATTTTATGTTAAACTAAGTTTTTT\tNA\t+\tchr1:100152443-100152562\n+mRNA|AL390972\tA_36_B232199\tCAGAGTCAATTGAGGACATCTATGCTAATATCCCTGACCTTCCAAAGTCCTACATACCTTCTAGGTTAAGGAAGGAGATAAATGAAGAAAAAGAAGATGATGAACAAAATAGGAAAGGTA\tNA\t+\tchr1:100152629-100152748\n+mRNA|AL390972\tA_36_B240075\tTTTACAAAAGAAAAGTTAACTAAATCAAAGATCAAATTTGTTTCTTAACAGCTTTATATGCCATGGAAATTAAAGTTGAAAAAGACTTGAAGACTGGAGAAAGTACAGTTCTGTCTTCAA\tNA\t+\tchr1:100154280-100154399\n+mRNA|AL390972\tA_36_B856165\tATCAAATTTGTTTCTTAACAGCTTTATATGCCATGGAAATTAAAGTTGAAAAAGACTTGAAGACTGGAGAAAGTACAGTTCTGTCTTCAATACCTCTGCCATCAGATGACTTTAAAGGTA\tNA\t+\tchr1:100154310-100154429\n+mRNA|AL390972\tA_36_B240076\tTACCTCTGCCATCAGATGACTTTAAAGGTACAGGAATAAAAGTTTATGATGATGGGCAAAAGTCAGTGTATGCAGTAAGTTCTAATCACAGTGCAGCATACAATGGCACCGATGGCCTGG\tNA\t+\tchr1:100154400-100154519\n+mRNA|AL390972\tA_36_B856501\tGTGCAGCATACAATGGCACCGATGGCCTGGCACCAGTTGAAGTAGAGGAACTTCTAAGACAAGCCTCAGAGAGAAACTCTAAATCCCCAACAGAGTATCATGAGCCTGTATATGCCAATC\tNA\t+\tchr1:100154490-100154609\n+mRNA|AL390972\tA_36_B240077\tCACCAGTTGAAGTAGAGGAACTTCTAAGACAAGCCTCAGAGAGAAACTCTAAATCCCCAACAGAGTATCATGAGCCTGTATATGCCAATCCCTTTTACAGGCCTACAACCCCACAGAGAG\tNA\t+\tchr1:100154520-100154639\n+mRNA|AL390972\tA_36_B240078\tAAACGGTGACCCCTGGACCAAACTTTCAAGAAAGGATAAAGATTAAAACTAATGGACTGGGTATTGGTGTAAATGAATCCATACACAATATGGGCAATGGTCTTTCAGAGGAAAGGGGAA\tNA\t+\tchr1:100154640-100154759\n+mRNA|AL390972\tA_36_B240079\tACAACTTCAATCACATCAGTCCCATTCCGCCAGTGCCTCATCCCCGATCAGTGATTCAACAAGCAGAAGAGAAGCTTCACACCCCGCAAAAAAGGCTAATGACTCCTTGGGAAGAATCGA\tNA\t+\tchr1:100154760-100154879\n+mRNA|AL390972\tA_36_B240080\tATGTCATGCAGGACAAAGATGCACCCTCTCCAAAGCCAAGGCTGAGCCCCAGAGAGACAATATTTGGGAAATCTGAACACCAGAATTCTTCACCCACTTGTCAGGAGGACGAGGAAGATG\tNA\t+\tchr1:100154880-100154999\n+mRNA|AL390972\tA_36_B240081\tTCAGATATAATATCGTTCATTCCCTGCCTCCAGACATAAATGATACAGAACCGGTGACAATGATTTTCATGGGGTATCAGCAGGCAGAAGACAGTGAAGAAGATAAGAAGTTTCTGACAG\tNA\t+\tchr1:100155000-100155119\n+mRNA|AL390972\tA_36_B240082\tGATATGATGGGATCATCCATGCTGAGCTGGTTGTGATTGATGATGAGGAGGAGGAGGATGAAGGAGAAGCAGAGAAACCGTCCTACCACCCCATAGCTCCCCATAGTCAGGTGTACCAGC\tNA\t+\tchr1:100155120-100155239\n+mRNA|AL390972\tA_36_B240083\tCAGCCAAACCAACACCACTTCCTAGAAAAAGATCAGAAGCTAGTCCTCATGAAAACACAAATCATAAATCCCCCCACAAAAATTCCATATCTCTGAAAGAGCAAGAAGAAAGCTTAGGCA\tNA\t+\tchr1:100155240-100155359\n+mRNA|AL390972\tA_36_B240084\tGCCCTGTCCACCATTCCCCATTTGATGCTCAGACAACTGGAGATGGGACTGAGGATCCATCCTTAACAGGTAATAAAAATGACAAGGCATGCCACTGCTGTTCAGTCATGTGACCAGCTT\tNA\t+\tchr1:100155360-100155479\n+mRNA|AL390972\tA_36_B238064\tTTTTTTTATAACTCTCTTTTTTTTTCTTTAATTTGCAGCTTTAAGGATGAGAATGGCAAAGCTGGGAAAAAAGGTGATCTAAGAGTTGTACCACCTATATAAACATCCTTTGAAGAAGAA\tNA\t+\tchr1:100159537-100159656\n+mRNA|BC156175\tA_36_B257433\tAGAATTCCTCTACAGACATGACCCCAGTCTTCCAAGTGAGAATTCACAAATATGCTCCAGGCAGTCAGGACAGGCTTCAAGTTTCTTTGGTTTGATGATAATTATCACTTGGCCTGCAAA\tNA\t+\tchr1:10'..b'GCATCAGATTTTACCTAATACAGCAGAACTCCTAAAAA\tNA\t+\tchrY:9368026-9368145\n+ens|ENST00000429039\tA_36_B883966\tAGTATGGCTCAGAGAAAATGCGTTTTAACATGAGTTTGTGTTTCTCTAGGGGACTCCCAGTTGTTGAGTTGAATATGATGGAGCATCAGATTTTACCTAATACAGCAGAACTCCTAAAAA\tNA\t+\tchrY:9368026-9368145\n+mRNA|U94385\tA_36_B104530\tCCCATGACCCCAGAGTCTGCACTGGAGGAGCTGCTGGCCGTTCAGGTGGAGCTGGAGCCGGTTAATGCCCAAGCCAGGAAGGCCTTTTCTCGGCAGCGGGAAAAGATGGAGCGGAGGCGC\tNA\t+\tchrY:9386007-9386126\n+mRNA|U94385\tA_36_B104620\tCCCATGACCCCAGAGTCTGCACTGGAGGAGCTGCTGGCCGTTCAGGTGGAGCTGGAGCCGGTTAATGCCCAAGCCAGGAAGGCCTTTTCTCGGCAGCGGGAAAAGATGGAGCGGAGGCGC\tNA\t+\tchrY:9386007-9386126\n+mRNA|U94385\tA_36_B883983\tCCCATGACCCCAGAGTCTGCACCGGAGGAGCTGCTGGCCGTTCAGGTGGAGCTGGAGCCGGTTAATGCCCAAGCCAGGAAGGCCTTTTCTCGGCAGCGGGAAAAGATGGAGCGGAGGCGC\tNA\t+\tchrY:9386007-9386126\n+mRNA|U94385\tA_36_B104531\tAAGCCCCACCTAGACCGCAGAGGCGCCGTCATCCAGAGCGTCCCTGGCTTCTGGGCCAATGTTGTATCCTTCTCAGTGTTTCTTCGGCCTTTCTAGTGGAGAGGTGCTCTCGGGGAAGTG\tNA\t+\tchrY:9386127-9386246\n+mRNA|U94385\tA_36_B104621\tAAGCCCCACCTAGACCGCAGAGGCGCCGTCATCCAGAGCGTCCCTGGCTTCTGGGCCAATGTTGTATCCTTCTCAGTGTTTCTTCGGCCTTTCTAGTGGAGAGGTGCTCTCGGGGAAGTG\tNA\t+\tchrY:9386127-9386246\n+mRNA|U94385\tA_36_B883984\tAAGCCCCACCTAGACCGCAGAGGCGCCGTCATCCAGAGCGTCCCTGGCTTCTGGGCCAATGTTGTATCCTTCTCAGTGTTTCTTCGGCCTTTCTAGTGGAGAGGTGCTCTCGGGGAAGTG\tNA\t+\tchrY:9386127-9386246\n+ens|ENST00000338964\tA_36_B104527\tTGCGGCCCCCCGGCCCTTCGCGCGCAGTCCCTTAGGGGGCGCCTGGAAGCCCGGCGCATGCGCCCTGAGGGCTCGCTGACCTACCGGGTGCCAGAGAGGCTGCGGCAGGGTTTCTGTGGC\tNA\t-\tchrY:9745686-9745805\n+ens|ENST00000338964\tA_36_B104617\tTGCGGCCCCCCGGCTCTTCGCGCGCAGTCCCTTAGGGGGCGCCTGGAAGCCCGGCGCATGCGCCCTGAGGGCTCGCTGACCTACCGGGTGCCAGAGAGGCTGCGGCAGGGTTCCTGTGGC\tNA\t-\tchrY:9745686-9745805\n+ens|ENST00000338964\tA_36_B883980\tTGCGGCCCCCCGGCCCTTCGCGCGCAGTCCCTTAGGGGGCGCCTGGAAGCCCGGCGCATGCGCCCTGAGGGCTCGCTGACCTACCGGGTGCCAGAGAGGCTGCGGCAGGGTTTCTGTGGC\tNA\t-\tchrY:9745686-9745805\n+ens|ENST00000429799\tA_36_B104615\tATGGTTCCATCTTTATATTTTAAAACATTACCCAGGTGCCCTTCTTGTGAGGGAAGCCACCCTCTTGTTCCTCCACGGCTTCCTCTTGCAGATCTCAGACTTCCTGAAGGGCTTCTGTTT\tNA\t-\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104730\tAAACAGAAGCCCTTCAGGAAGTCTGAGATCTGCAAGAGGAAGCCGTGGAGGAACAAGAGGGTGGCTTCCCTCACAAGAAGGGCACCTGGGTAATGTTTTAAAATATAAAGATGGAACCAT\tNA\t+\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104744\tATGGTTCCATCTTTATATTTTAAAACATTACCCAGGTGCCCTTCTTGTGAGGGAAGCCACCCTCTTGTTCCTCCACGGCTTCCTCTTGCAGATCTCAGACTTCCTGAAGGGCTTCTGTTT\tNA\t-\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104847\tAAACAGAAGCCCTTCAGGAAGTCTGAGATCTGCAAGAGGAAGCAGTGGAGGAACAAGAGGGTGGCTTCCCTCACATGAAGGGCACCTGGGTAATGTTTTAAAATATAAAGATGGAACCAT\tNA\t+\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104907\tAAACAGAAGCCCTTCAGGAAGTCTGAGATCTGCAAGAGGAAGCCGTGGAGGAACAAGAGGGTGGCTTCCCTCACATGAAGGGCACCTGGGTAATGTTTTAAAATATAAAGATGGAACCAT\tNA\t+\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104969\tATGGTTCCATCTTTATATTTTAAAACATTACCCAGGTGCCCTTCATGTGAGGGAAGCCACCCTCTTGTTCCTCCACTGCTTCCTCTTGCAGATCTCAGACTTCCTGAAGGGCTTCTGTTT\tNA\t-\tchrY:9867967-9868086\n+ens|ENST00000429799\tA_36_B104616\tCTCGAAGAAGCTGGTGGTCTCCGCCTACCACCACTTTGAAAAGATGGTTTCTTGGCTTGTTCTACTTTTATTGCTTTTCCATGCAAAGACTAGAAGTATTAAGGGTACTATCAATAACAC\tNA\t-\tchrY:9868087-9868206\n+ens|ENST00000429799\tA_36_B104729\tGTGTTATTGATAGTACCCTTAATACTTCTAGTCTTTGCATGGAAAAGCAATAAAAGTAGAACAAGCCAAGAAACCATCTTTTCAAAGTGGTGGTAGGCGGAGACCACCAGCTTCTTCGAG\tNA\t+\tchrY:9868087-9868206\n+ens|ENST00000429799\tA_36_B104745\tCTCGAAGAAGCTGGTGGTCTCCGCCTACCACCACTTTGAAAAGATGGTTTCTTGGCTTGTTCTACTTTTATTGCTTTTCCATGCAAAGACTAGAAGTATTAAGGGTACTATCAATAACAC\tNA\t-\tchrY:9868087-9868206\n+ens|ENST00000429799\tA_36_B104846\tGTGTTATTGATAGTACCCTTAATACTTCTAGTCTTTGCATGGAAAAGCAATAAAAGTAGAACAAGCCAAGAAACCATCTTTTCAAAGTGGTGGTAGGCGGAGACCACCAGCTTCTTCGAG\tNA\t+\tchrY:9868087-9868206\n+ens|ENST00000429799\tA_36_B104906\tGTGTTATTGATAGTACCCTTAATACTTCTAGTCTTTGCATGGAAAAGCAATAAAAGTAGAACAAGCCAAGAAACCATCTTTTCAAAGTGGTGGTAGGCGGAGACCACCAGCTTCTTCGAG\tNA\t+\tchrY:9868087-9868206\n+ens|ENST00000429799\tA_36_B104970\tCTCGAAGAAGCTGGTGGTCTCCGCCTACCACCACTTTGAAAAGATGGTTTCTTGGCTTGTTCTACTTTTATTGCTTTTCCATGCAAAGACTAGAAGTATTAAGGGTACTATCAATAACAC\tNA\t-\tchrY:9868087-9868206\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/lib/mvnTable.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/Wessim/lib/mvnTable.txt Sat Aug 06 15:27:30 2016 -0400 |
| b |
| b'@@ -0,0 +1,1 @@\n+0.3989422804014327\t0.3989223337860822\t0.3988624999236661\t0.3987627967620997\t0.398623254204605\t0.398443914094764\t0.3982248301956069\t0.397966068162751\t0.3976677055116088\t0.3973298315786883\t0.3969525474770118\t0.3965359660456858\t0.3960802117936561\t0.3955854208376874\t0.3950517408346113\t0.3944793309078889\t0.3938683615685408\t0.3932190146304972\t0.3925314831204289\t0.3918059711821211\t0.3910426939754559\t0.3902418775700743\t0.3894037588337904\t0.3885285853158359\t0.3876166151250142\t0.3866681168028492\t0.3856833691918161\t0.3846626612987428\t0.3836062921534786\t0.382514570662924\t0.3813878154605241\t0.3802263547513249\t0.3790305261527017\t0.3778006765308646\t0.3765371618332539\t0.3752403469169379\t0.3739106053731284\t0.3725483193479334\t0.371153879359466\t0.3697276841114324\t0.3682701403033233\t0.3667816624373361\t0.3652626726221539\t0.3637136003737134\t0.3621348824130922\t0.3605269624616479\t0.3588902910335446\t0.3572253252258008\t0.3555325285059971\t0.3538123704977797\t0.3520653267642995\t0.3502918785897258\t0.3484925127589745\t0.3466677213357917\t0.3448180014393333\t0.3429438550193839\t0.3410457886303526\t0.3391243132041922\t0.3371799438223806\t0.3352131994871061\t0.3332246028917997\t0.331214680191153\t0.3291839607707648\t0.3271329770165545\t0.3250622640840822\t0.3229723596679143\t0.3208638037711725\t0.3187371384754015\t0.3165929077108928\t0.3144316570275973\t0.3122539333667613\t0.3100602848334161\t0.3078512604698529\t0.3056274100302099\t0.3033892837563001\t0.3011374321548044\t0.2988724057759528\t0.2965947549938157\t0.2943050297883251\t0.2920037795291414\t0.2896915527614827\t0.2873688969940283\t0.2850363584890072\t0.2826944820545803\t0.2803438108396206\t0.2779848861309965\t0.2756182471534567\t0.2732444308722162\t0.270863971798338\t0.2684774017970024\t0.2660852498987548\t0.2636880421138181\t0.2612863012495532\t0.2588805467311488\t0.2564712944256203\t0.254059056469189\t0.2516443410981171\t0.2492276524830659\t0.2468094905670427\t0.2443903509069996\t0.2419707245191434\t0.2395510977280134\t0.2371319520193796\t0.2347137638970118\t0.2322970047433662\t0.229882140684233\t0.2274696324573859\t0.2250599352852697\t0.2226534987517611\t0.2202507666830333\t0.2178521770325505\t0.2154581617702197\t0.2130691467757178\t0.2106855517360152\t0.2083077900471083\t0.2059362687199747\t0.2035713882907594\t0.2012135427351974\t0.1988631193872759\t0.1965204988621365\t0.194186054983213\t0.1918601547135994\t0.1895431580916402\t0.1872354181707296\t0.1849372809633053\t0.1826490853890219\t0.1803711632270803\t0.1781038390726936\t0.1758474302976624\t0.173602247015033\t0.1713685920478074\t0.1691467609016724\t0.1669370417417138\t0.1647397153730768\t0.1625550552255341\t0.1603833273419196\t0.158224790370383\t0.1560796955604208\t0.1539482867626337\t0.1518308004321616\t0.1497274656357448\t0.1476385040623557\t0.1455641300373476\t0.1435045505400624\t0.1414599652248388\t0.1394305664453603\t0.1374165392822818\t0.1354180615740713\t0.1334353039510023\t0.131468429872231\t0.1295175956658917\t0.1275829505721419\t0.1256646367890882\t0.1237627895215231\t0.1218775370324018\t0.1200090006969856\t0.1181572950595823\t0.1163225278928071\t0.1145048002592924\t0.1127042065757706\t0.1109208346794555\t0.1091547658966474\t0.1074060751134838\t0.1056748308487636\t0.1039610953287642\t0.102264924563978\t0.1005863684276905\t0.09892547073632371\t0.09728226933146751\t0.09565679616352402\t0.09404907737688695\t0.09245913339658068\t0.09088697901628287\t0.089332623487655\t0.08779607061090562\t0.08627731882651153\t0.08477636130802223\t0.08329318605587446\t0.0818277759921428\t0.08038010905615417\t0.07895015830089415\t0.07753789199013399\t0.07614327369620731\t0.0747662623983676\t0.07340681258165689\t0.07206487433621798\t0.07074039345698338\t0.06943331154367419\t0.06814356610104455\t0.06687109063930714\t0.06561581477467658\t0.06437766432996934\t0.06315656143519865\t0.06195242462810516\t0.06076516895456478\t0.05959470606881608\t0.05844094433345147\t0.05730378891911713\t0.05618314190386805\t0.05507890237212577\t0.05399096651318806\t0.05291922771924027\t0.05186357668282057\t0.05082390149369116\t0.04980008773507077\t0.04879201857918276\t0.04779957488207703\t0.04682263527768316\t0.04586'..b'7\t0.002623126024781024\t0.002541150028726521\t0.002461489724640701\t0.00238408820146484\t0.002308889668006496\t0.002235839439688538\t0.002164883925171062\t0.002095970612857942\t0.002029048057299768\t0.001964065865504374\t0.00190097468316608\t0.001839726180824278\t0.001780273039961879\t0.001722568939053678\t0.00166656853957458\t0.001612227471977123\t0.001559502321647691\t0.001508350614850307\t0.001458730804666746\t0.001410602256941385\t0.001363925236238904\t0.001318660891822742\t0.001274771243661833\t0.00123221916847302\t0.001190968385806117\t0.001150983444178485\t0.001112229707265565\t0.001074673340153736\t0.00103828129566141\t0.001003021300734238\t0.0009688618429198459\t0.0009357721569274798\t0.0009037222112775245\t0.0008726826950457602\t0.0008426250047069012\t0.0008135212310818084\t0.0007853441463924686\t0.0007580671914287103\t0.0007316644628303096\t0.0007061107004880362\t0.0006813812750668915\t0.0006574521756546765\t0.0006342999975387576\t0.0006119019301137719\t0.0005902357449227856\t0.0005692797838342526\t0.0005490129473569587\t0.0005294146830949348\t0.0005104649743441856\t0.0004921443288328931\t0.0004744337676066206\t0.0004573148140598567\t0.0004407694831151325\t0.0004247802705507514\t0.0004093301424780788\t0.0003944025249691562\t0.0003799812938353214\t0.000366050764557335\t0.0003525956823674454\t0.0003396012124836542\t0.000327052930496375\t0.0003149368129075216\t0.0003032392278220042\t0.00029194692579146\t0.0002810470308099863\t0.0002705270314615207\t0.0002603747722184425\t0.0002505784448908607\t0.0002411265802259932\t0.0002320080396569424\t0.0002232120072001021\t0.000214727981500367\t0.0002065457680232255\t0.0001986554713927727\t0.0001910474878745976\t0.0001837124980024571\t0.0001766414593475709\t0.0001698255994293436\t0.000163256408766242\t0.0001569256340655323\t0.0001508252715505178\t0.0001449475604238911\t0.0001392849764657599\t0.0001338302257648854\t0.0001285762385816211\t0.0001235161633410232\t0.0001186433607545658\t0.0001139513980688646\t0.0001094340434398006\t0.0001050852604304001\t0.0001008992026308144\t9.687020839871926e-05\t9.299279571844591e-05\t8.926165717713293e-05\t8.567165505618186e-05\t8.2217816536286e-05\t7.889532901429309e-05\t7.569953553016121e-05\t7.262593030225232e-05\t6.967015436921433e-05\t6.682799133669061e-05\t6.40953632271061e-05\t6.146832643076922e-05\t5.894306775653985e-05\t5.651590058030741e-05\t5.418326108954014e-05\t5.194170462215977e-05\t4.978790209801209e-05\t4.771863654120495e-05\t4.573079969160131e-05\t4.382138870375796e-05\t4.198750293161732e-05\t4.022634079726497e-05\t3.853519674208713e-05\t3.691145825866607e-05\t3.53526030017731e-05\t3.385619597682789e-05\t3.241988680421378e-05\t3.104140705785016e-05\t2.97185676764422e-05\t2.84492564458443e-05\t2.723143555099261e-05\t2.606313919587834e-05\t2.494247129005354e-05\t2.38676032001796e-05\t2.283677156514692e-05\t2.184827617331647e-05\t2.090047790045041e-05\t1.999179670692279e-05\t1.912070969281774e-05\t1.828574920954738e-05\t1.748550102663914e-05\t1.671860255236507e-05\t1.598374110690547e-05\t1.527965224676162e-05\t1.460511813915286e-05\t1.395896598515477e-05\t1.334006649035584e-05\t1.274733238183347e-05\t1.217971697026866e-05\t1.163621275604267e-05\t1.111585007817779e-05\t1.061769580500839e-05\t1.014085206548672e-05\t9.684455020051437e-06\t9.247673670005617e-06\t8.829708704374098e-06\t8.429791383228772e-06\t8.047182456492295e-06\t7.681171117250455e-06\t7.331073988623945e-06\t6.996234143270405e-06\t6.676020154607461e-06\t6.36982517886709e-06\t6.077066067111115e-06\t5.797182506357287e-06\t5.529636188984051e-06\t5.273910009601303e-06\t5.029507288592445e-06\t4.795951021552522e-06\t4.572783153864128e-06\t4.359563879671637e-06\t4.155870964531201e-06\t3.961299091032075e-06\t3.775459226701329e-06\t3.597978013521247e-06\t3.428497178405039e-06\t3.266672963993275e-06\t3.112175579148934e-06\t2.964688668545273e-06\t2.823908800755821e-06\t2.689544974271523e-06\t2.561318140884544e-06\t2.438960745893352e-06\t2.322216284597997e-06\t2.210838874568421e-06\t2.104592843183129e-06\t2.003252329948489e-06\t1.906600903122811e-06\t1.81443119018203e-06\t1.726544521677073e-06\t1.642750588045071e-06\t1.56286710894929e-06\t1.486719514734298e-06\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/models/ill100v4_p.gzip |
| b |
| Binary file cnvsim/Wessim/models/ill100v4_p.gzip has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/models/ill100v4_s.gzip |
| b |
| Binary file cnvsim/Wessim/models/ill100v4_s.gzip has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/models/ill100v5_p.gzip |
| b |
| Binary file cnvsim/Wessim/models/ill100v5_p.gzip has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/models/ill100v5_s.gzip |
| b |
| Binary file cnvsim/Wessim/models/ill100v5_s.gzip has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/Wessim/models/r454ti_s.gzip |
| b |
| Binary file cnvsim/Wessim/models/r454ti_s.gzip has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/__init__.pyc |
| b |
| Binary file cnvsim/__init__.pyc has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/exome_simulator.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/exome_simulator.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| b'@@ -0,0 +1,314 @@\n+#!/usr/bin/python\n+\n+__author__ = \'Abdelrahman Hosny\'\n+\n+import random\n+import subprocess\n+import os\n+import sys\n+import datetime\n+import shutil\n+\n+from . import fileio\n+\n+\n+def _log(message):\n+ print \'[CNV SIM {:%Y-%m-%d %H:%M:%S}\'.format(datetime.datetime.now()) + "] " + message\n+\n+\n+def getScriptPath():\n+ return os.path.dirname(os.path.realpath(__file__))\n+\n+\n+def _generateCNVMatrix(targets_list, regions_count):\n+ \'\'\'\n+ A function to randomly divide sequential exonic targets into whole CNV regions\n+ :param targets_list: a list of target exons loaded from the file provided by the user\n+ :return: A matrix where rows represent the region index and the first column as a list of targets in this region\n+ \'\'\'\n+ regions_count -= 1\n+ number_of_targets = len(targets_list)\n+ comine_size = number_of_targets / regions_count\n+ cnv_matrix = []\n+\n+ i = 0\n+ while True:\n+ if i == len(targets_list):\n+ break\n+\n+ first_target = i\n+ i += comine_size\n+ if i > len(targets_list):\n+ i = len(targets_list)\n+ last_target = i\n+ cnv_matrix.append(targets_list[first_target:last_target])\n+\n+ return cnv_matrix\n+\n+\n+def _loadCNVMatrix(targets_list, cnv_regions):\n+ \'\'\'\n+ A function to map targets to their corresponding CNV regions\n+ :param targets_list: a list of target exons loaded from the file provided by the user\n+ :param cnv_regions: a list of CNV regions loaded from the CNV file provided by the user\n+ :return: A matrix where rows represent the region index and the first column as a list of targets in this region\n+ \'\'\'\n+ cnv_matrix = []\n+ for cnv in cnv_regions:\n+ region_chromosome = cnv[0]\n+ region_start = cnv[1]\n+ region_end = cnv[2]\n+\n+ region_targets = []\n+ for target in targets_list:\n+ target_chromosome = target[0]\n+ target_start = target[1]\n+ target_end = target[2]\n+ if target_chromosome == region_chromosome and target_start >= region_start and target_end <= region_end:\n+ region_targets.append(target)\n+\n+ cnv_matrix.append(region_targets)\n+\n+ return cnv_matrix\n+\n+\n+def _generateCNVMask(number_of_exons, p_amplify, p_delete, min_variation, max_variation):\n+ \'\'\'\n+ This function generates random Copy Number Variations mask list\n+ :param number_of_exons: the length of the target regions.\n+ :param p_amplify: percentage of amplifications introduced\n+ :param p_delete: percentage of deletions introduced\n+ :return: a list of the same length as the exons list. each list item\n+ indicates the variation to be added to the exon in the same position.\n+ Positive number: number of amplification\n+ Zero: no change\n+ -1: delete this exon\n+ \'\'\'\n+\n+ number_of_amplifications = int(p_amplify * number_of_exons)\n+ number_of_deletions = int(p_delete * number_of_exons)\n+ cnv_mask = [0] * number_of_exons\n+\n+ # generate CNV mask (a list of amplifications and deletions to be applied to the exons)\n+ while number_of_amplifications > 0:\n+ choice = random.randrange(0, number_of_exons)\n+ while cnv_mask[choice] != 0:\n+ choice = random.randrange(0, number_of_exons)\n+ cnv_mask[choice] = random.randrange(min_variation, max_variation) # random amplifications in the range [min_variation, max_variation)\n+ number_of_amplifications -= 1\n+ random.shuffle(cnv_mask)\n+ while number_of_deletions > 0:\n+ choice = random.randrange(0, number_of_exons)\n+ while cnv_mask[choice] != 0:\n+ choice = random.randrange(0, number_of_exons)\n+ cnv_mask[choice] = -1*random.randrange(min_variation, max_variation) # random deletions in the range [min_variation, max_variation)\n+ number_of_deletions -= 1\n+ random.shuffle(cnv_mask)\n+\n+ return cnv_mask\n+\n+\n+def _simulateCNV(genome, cnv_matrix, mas'..b't_dir\'], "CNVList.bed")\n+\n+ _log("generating CNV list ..")\n+ cnv_matrix = _generateCNVMatrix(targets, cnv_list_parameters[\'regions_count\'])\n+ mask = _generateCNVMask(len(cnv_matrix), cnv_list_parameters[\'amplifications\'], \\\n+ cnv_list_parameters[\'deletions\'], cnv_list_parameters[\'minimum_variations\'], \\\n+ cnv_list_parameters[\'maximum_variations\'])\n+\n+ with open(cnv_list_file, \'w\') as f:\n+ for i, cnv_region in enumerate(cnv_matrix):\n+ line = cnv_region[0][0] + \'\\t\' \\\n+ + `cnv_region[0][1]` + \'\\t\' \\\n+ + `cnv_region[-1][2]` + \'\\t\' \\\n+ + `mask[i]` + \'\\n\'\n+ f.write(line)\n+ _log("randomly generated CNV list saved to " + cnv_list_file)\n+ else:\n+ _log("loading CNV list ..")\n+ with open(simulation_parameters[\'cnv_list_file\'], "r") as f:\n+ cnv_list = []\n+ lines = f.readlines()\n+ for line in lines:\n+ chromosome, region_start, region_end, variation = line.strip().split("\\t")\n+ cnv_list.append((chromosome, int(region_start), int(region_end), int(variation)))\n+\n+ cnv_matrix = _loadCNVMatrix(targets, cnv_list)\n+ mask = map(lambda x: x[3], cnv_list)\n+ _log("successfully loaded CNV list that contains " + `len(lines)` + " regions ..")\n+\n+ # call Wessim to generate reads from the genome file and the target list\n+ _log("generating reads for the target exons ..")\n+ _log("delegating job to Wessim ...")\n+ _callWessim(genome_file, target_file, base_reads_file, simulation_parameters[\'number_of_reads\'], simulation_parameters[\'read_length\'])\n+\n+ # generate genome and target files for the control and the CNV\n+ _log("simulating copy number variations (amplifications/deletions)")\n+ control_genome, control_targets, cnv_genome, cnv_targets = _simulateCNV(genome, cnv_matrix, mask, simulation_parameters[\'read_length\'])\n+\n+ _log("saving to the control genome file ..")\n+ with open(control_genome_file, \'w\') as fw:\n+ fw.write(header + "\\n")\n+ n = 50\n+ l = [control_genome[i:i + n] for i in range(0, len(control_genome), n)]\n+ for line in l:\n+ fw.write(line + "\\n")\n+\n+ _log("saving to the control target file ..")\n+ with open(control_target_file, \'w\') as tw:\n+ for target in control_targets:\n+ line = target[0] + "\\t" + `target[1]` + "\\t" + `target[2]` + "\\n"\n+ tw.write(line)\n+\n+ _log("delegating job to Wessim ...")\n+ _callWessim(control_genome_file, control_target_file, control_reads_file,\n+ int(cnv_list_parameters[\'deletions\'] * simulation_parameters[\'number_of_reads\']), simulation_parameters[\'read_length\'])\n+\n+ _log("saving to the CNV genome file ..")\n+ with open(cnv_genome_file, \'w\') as fw:\n+ fw.write(header + "\\n")\n+ n = 50\n+ l = [cnv_genome[i:i + n] for i in range(0, len(cnv_genome), n)]\n+ for line in l:\n+ fw.write(line + "\\n")\n+\n+ _log("saving to the CNV target file ..")\n+ with open(cnv_target_file, \'w\') as tw:\n+ for target in cnv_targets:\n+ line = target[0] + "\\t" + `target[1]` + "\\t" + `target[2]` + "\\n"\n+ tw.write(line)\n+\n+ _log("delegating job to Wessim ...")\n+ _callWessim(cnv_genome_file, cnv_target_file, cnv_reads_file,\n+ int(cnv_list_parameters[\'amplifications\']* simulation_parameters[\'number_of_reads\']), simulation_parameters[\'read_length\'])\n+\n+ _log("merging results ..")\n+ fileio.mergeWessimReads(simulation_parameters[\'tmp_dir\'], simulation_parameters[\'output_dir\'])\n+\n+ _log("cleaning temporary files ..")\n+ fileio.clean(simulation_parameters[\'tmp_dir\'])\n+\n+ _log("simulation completed. find results in " + simulation_parameters[\'output_dir\'])\n\\ No newline at end of file\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/exome_simulator.pyc |
| b |
| Binary file cnvsim/exome_simulator.pyc has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/fileio.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/fileio.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| @@ -0,0 +1,100 @@ +#!/usr/bin/python + +__author__ = 'Abdelrahman Hosny' + +import subprocess +import os + +def readGenome(filename): + ''' + Read a FASTA genome + :param filename: genome fasta file name + :return: a header string and a genome string + ''' + genome = '' + header = '' + with open(filename, 'r') as f: + for line in f: + # ignore header line with genome information + if not line[0] == '>': + genome += line.rstrip() + else: + header = line.strip() + return header, genome + + +def readTargets(filename): + ''' + Read target file in BED format + :param filename: target bed file name + :return: a list of targets (exons). each target is a dictionary of: + chromosome, start, end, description, score, strand + ''' + exons = [] + with open(filename, 'r') as tf: + for line in tf: + chromosome, start, end = line.strip().split('\t') + exon = (chromosome, int(start), int(end)) + exons.append(exon) + return exons + +def prepareTargetFile(target_file): + ''' + sort and merge targets in the target file and writes it to '.sorted.merged' + :param target_file: target exons in BED format + :return: sorted and merged file name + ''' + sorted_file = target_file + ".sorted" + merged_file = sorted_file + ".merged" + + with open(sorted_file, "w") as f: + subprocess.call(["sort", "-k1,1", "-k2,2n", target_file], stdout=f) + with open(merged_file, "w") as f: + subprocess.call(["bedtools", "merge", "-i", sorted_file], stdout=f) + + return merged_file + +def mergeWessimReads(tmp_dir, output_dir): + ''' + merges the base reads with normal and cnv + :return: null + ''' + base_file_1 = os.path.join(tmp_dir, "base_1.fastq") + base_file_2 = os.path.join(tmp_dir, "base_2.fastq") + normal_file_1 = os.path.join(tmp_dir, "control_1.fastq") + normal_file_2 = os.path.join(tmp_dir, "control_2.fastq") + cnv_file_1 = os.path.join(tmp_dir, "cnv_1.fastq") + cnv_file_2 = os.path.join(tmp_dir, "cnv_2.fastq") + + with open(os.path.join(output_dir, "control_1.fastq"), "w") as f: + subprocess.call(["cat", base_file_1, normal_file_1], stdout=f) + with open(os.path.join(output_dir, "control_2.fastq"), "w") as f: + subprocess.call(["cat", base_file_2, normal_file_2], stdout=f) + with open(os.path.join(output_dir, "cnv_1.fastq"), "w") as f: + subprocess.call(["cat", base_file_1, cnv_file_1], stdout=f) + with open(os.path.join(output_dir, "cnv_2.fastq"), "w") as f: + subprocess.call(["cat", base_file_2, cnv_file_2], stdout=f) + +def mergeARTReads(tmp_dir, output_dir): + ''' + merges the base reads with normal and cnv + :return: null + ''' + base_file_1 = os.path.join(tmp_dir, "base1.fq") + base_file_2 = os.path.join(tmp_dir, "base2.fq") + normal_file_1 = os.path.join(tmp_dir, "control1.fq") + normal_file_2 = os.path.join(tmp_dir, "control2.fq") + cnv_file_1 = os.path.join(tmp_dir, "cnv1.fq") + cnv_file_2 = os.path.join(tmp_dir, "cnv2.fq") + + with open(os.path.join(output_dir, "control_1.fastq"), "w") as f: + subprocess.call(["cat", base_file_1, normal_file_1], stdout=f) + with open(os.path.join(output_dir, "control_2.fastq"), "w") as f: + subprocess.call(["cat", base_file_2, normal_file_2], stdout=f) + with open(os.path.join(output_dir, "cnv_1.fastq"), "w") as f: + subprocess.call(["cat", base_file_1, cnv_file_1], stdout=f) + with open(os.path.join(output_dir, "cnv_2.fastq"), "w") as f: + subprocess.call(["cat", base_file_2, cnv_file_2], stdout=f) + +def clean(tmp_dir): + subprocess.call(["rm", "-rf", tmp_dir]) |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/fileio.pyc |
| b |
| Binary file cnvsim/fileio.pyc has changed |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/genome_simulator.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cnvsim/genome_simulator.py Sat Aug 06 15:27:30 2016 -0400 |
| [ |
| b'@@ -0,0 +1,236 @@\n+__author__ = \'Abdelrahman Hosny\'\n+\n+import random\n+import os\n+import subprocess\n+import datetime\n+import shutil\n+\n+from . import fileio\n+\n+\n+def _log(message):\n+ print \'[CNV SIM {:%Y-%m-%d %H:%M:%S}\'.format(datetime.datetime.now()) + "] " + message\n+\n+\n+def getScriptPath():\n+ return os.path.dirname(os.path.realpath(__file__))\n+\n+\n+def _generateCNVMask(mask_length, p_amplify, p_delete, min_variation, max_variation):\n+ \'\'\'\n+ This function generates random Copy Number Variations mask list\n+ :param exons: list of regions.\n+ :param p_amplify: percentage of amplifications introduced\n+ :param p_delete: percentage of deletions introduced\n+ :return: a list of the same length as the exons list. each list item\n+ indicates the variation to be added to the exon in the same position.\n+ Positive number: number of amplification\n+ Zero: no change\n+ -1: delete this exon\n+ \'\'\'\n+\n+ number_of_amplifications = int(p_amplify * mask_length)\n+ number_of_deletions = int(p_delete * mask_length)\n+ cnv_mask = [0] * mask_length\n+\n+ # generate CNV mask (a list of amplifications and deletions to be applied to the exons)\n+ while number_of_amplifications > 0:\n+ choice = random.randrange(0, mask_length)\n+ while cnv_mask[choice] != 0:\n+ choice = random.randrange(0, mask_length)\n+ cnv_mask[choice] = random.randrange(min_variation, max_variation) # random amplifications in the range [min_variation, max_variation)\n+ number_of_amplifications -= 1\n+ random.shuffle(cnv_mask)\n+ while number_of_deletions > 0:\n+ choice = random.randrange(0, mask_length)\n+ while cnv_mask[choice] != 0:\n+ choice = random.randrange(0, mask_length)\n+ cnv_mask[choice] = -1*random.randrange(min_variation, max_variation) # random deletions in the range [min_variation, max_variation)\n+ number_of_deletions -= 1\n+ random.shuffle(cnv_mask)\n+\n+ return cnv_mask\n+\n+\n+def _generateCNVList(chromosome_length, number_of_regions, p_amplify, p_delete, min_variation, max_variation):\n+ \'\'\'\n+\n+ :param chromosome_length:\n+ :param number_of_regions:\n+ :param p_amplify:\n+ :param p_delete:\n+ :param min_variation:\n+ :param max_variation:\n+ :return:\n+ \'\'\'\n+ region_length = chromosome_length / number_of_regions\n+ start = 0\n+ end = region_length\n+ cnv_list = []\n+ mask = _generateCNVMask(number_of_regions, p_amplify, p_delete, min_variation, max_variation)\n+ for i in range(number_of_regions):\n+ # jump forward start\n+ jump_start = start + int(random.randrange(int(0.40*region_length), int(0.45*region_length)))\n+\n+ # jump backward end\n+ jump_end = end - int(random.randrange(int(0.40*region_length), int(0.45*region_length)))\n+\n+ cnv_list.append([jump_start, jump_end, mask[i]])\n+ start = end\n+ end += region_length\n+\n+ return cnv_list\n+\n+\n+def _simulateCNV(genome, cnv_list, read_length):\n+ \'\'\'\n+ Simulate the control genome and the CNV genome\n+ :param genome: original genome sequence\n+ :param cnv_list: a list of region variations (chromosome, start, end, variation)\n+ :return: control_genome, cnv_genome\n+ \'\'\'\n+ control_genome = []\n+ cnv_genome = []\n+\n+ for cnv in cnv_list:\n+ start = cnv[1]\n+ end = cnv[2]\n+ variation = cnv[3]\n+ sequence = genome[start: end]\n+\n+ if variation > 0:\n+ # amplification\n+ amplification = genome[start-read_length: start] + sequence * variation + genome[end:end+read_length]\n+ cnv_genome.append(amplification)\n+ elif variation < 0:\n+ # deletion\n+ deletion = genome[start-read_length: start] + sequence * variation + genome[end:end+read_length]\n+ control_genome.append(deletion)\n+\n+ return \'\'.join(control_genome), \'\'.join(cnv_genome)\n+\n+\n+def _callART(genome_file, output_file, r'..b' # copy genome to the tmp folder\n+ shutil.copyfile(simulation_parameters[\'genome_file\'], os.path.join(simulation_parameters[\'tmp_dir\'], "reference.fa"))\n+ genome_file = os.path.join(simulation_parameters[\'tmp_dir\'], "reference.fa")\n+\n+ # initialize variables for temporary files\n+ control_genome_file = os.path.join(simulation_parameters[\'tmp_dir\'], "ControlGenome.fa")\n+ cnv_genome_file = os.path.join(simulation_parameters[\'tmp_dir\'], "CNVGenome.fa")\n+ base_reads_file = os.path.join(simulation_parameters[\'tmp_dir\'], "base")\n+ control_reads_file = os.path.join(simulation_parameters[\'tmp_dir\'], "control")\n+ cnv_reads_file = os.path.join(simulation_parameters[\'tmp_dir\'], "cnv")\n+\n+ _log("loading genome file ..")\n+ header, genome = fileio.readGenome(genome_file)\n+ _log("successfully loaded a genome of length " + `len(genome)`)\n+\n+ if simulation_parameters[\'cnv_list_file\'] is None:\n+ # CNV list file\n+ cnv_list_file = os.path.join(simulation_parameters[\'output_dir\'], "CNVList.bed")\n+\n+ _log("generating CNV list ..")\n+ cnv_list = _generateCNVList(len(genome), cnv_list_parameters[\'regions_count\'], \\\n+ cnv_list_parameters[\'amplifications\'], cnv_list_parameters[\'deletions\'], \\\n+ cnv_list_parameters[\'minimum_variations\'], \\\n+ cnv_list_parameters[\'maximum_variations\'])\n+ cnv_list = map(lambda l: [header.replace(\'>\', \'\')] + l, cnv_list)\n+\n+ with open(cnv_list_file, \'w\') as f:\n+ for i, cnv_region in enumerate(cnv_list):\n+ line = cnv_region[0] + \'\\t\' \\\n+ + `cnv_region[1]` + \'\\t\' \\\n+ + `cnv_region[2]` + \'\\t\' \\\n+ + `cnv_region[3]` + \'\\n\'\n+ f.write(line)\n+ _log("randomly generated CNV list saved to " + cnv_list_file)\n+\n+ else:\n+ _log("loading CNV list ..")\n+ with open(simulation_parameters[\'cnv_list_file\'], "r") as f:\n+ cnv_list = []\n+ lines = f.readlines()\n+ for line in lines:\n+ chromosome, region_start, region_end, variation = line.strip().split("\\t")\n+ cnv_list.append((chromosome, int(region_start), int(region_end), int(variation)))\n+\n+ _log("successfully loaded CNV list that contains " + `len(lines)` + " regions ..")\n+\n+ # call ART to generate reads from the genome file\n+ _log("generating reads for the genome ..")\n+ _log("delegating job to ART ...")\n+ _callART(genome_file, base_reads_file, simulation_parameters[\'read_length\'])\n+\n+ _log("simulating copy number variations (amplifications/deletions)")\n+ control_genome, cnv_genome = _simulateCNV(genome, cnv_list, simulation_parameters[\'read_length\'])\n+\n+ _log("saving to the control genome file ..")\n+\n+ with open(control_genome_file, \'w\') as fw:\n+ fw.write(header + "\\n")\n+ n = 50\n+ l = [control_genome[i:i + n] for i in range(0, len(control_genome), n)]\n+ for line in l:\n+ fw.write(line + "\\n")\n+\n+ _log("delegating job to ART ...")\n+ _callART(control_genome_file, control_reads_file, simulation_parameters[\'read_length\'])\n+\n+ _log("saving to the CNV genome file ..")\n+ with open(cnv_genome_file, \'w\') as fw:\n+ fw.write(header + "\\n")\n+ n = 50\n+ l = [cnv_genome[i:i + n] for i in range(0, len(cnv_genome), n)]\n+ for line in l:\n+ fw.write(line + "\\n")\n+\n+ _log("delegating job to ART ...")\n+ _callART(cnv_genome_file, cnv_reads_file, simulation_parameters[\'read_length\'])\n+\n+ _log("merging results ..")\n+ fileio.mergeARTReads(simulation_parameters[\'tmp_dir\'], simulation_parameters[\'output_dir\'])\n+\n+ _log("cleaning temporary files ..")\n+ fileio.clean(simulation_parameters[\'tmp_dir\'])\n+\n+ _log("simulation completed. find results in " + simulation_parameters[\'output_dir\'])\n\\ No newline at end of file\n' |
| b |
| diff -r 4a4d2b78eb55 -r 63955244b2fa cnvsim/genome_simulator.pyc |
| b |
| Binary file cnvsim/genome_simulator.pyc has changed |