[med-svn] [Git][med-team/dwgsim][master] 9 commits: New upstream version 0.1.12

Andreas Tille gitlab at salsa.debian.org
Thu Mar 29 08:45:33 UTC 2018


Andreas Tille pushed to branch master at Debian Med / dwgsim


Commits:
9ed75949 by Andreas Tille at 2018-03-29T09:40:12+02:00
New upstream version 0.1.12
- - - - -
8c7690ab by Andreas Tille at 2018-03-29T09:40:14+02:00
Update upstream source from tag 'upstream/0.1.12'

Update to upstream version '0.1.12'
with Debian dir 5bf86c51d6bec4bf937ab3593f0d0aa4afc1cc62
- - - - -
365a216d by Andreas Tille at 2018-03-29T09:40:40+02:00
New upstream version

- - - - -
8c578d69 by Andreas Tille at 2018-03-29T09:41:12+02:00
cme fix dpkg-control

- - - - -
77e875c7 by Andreas Tille at 2018-03-29T09:41:43+02:00
debhelper 11

- - - - -
3333c667 by Andreas Tille at 2018-03-29T10:10:39+02:00
build time test needs example data from samtools

- - - - -
7c3c2502 by Andreas Tille at 2018-03-29T10:22:13+02:00
Upstream converted README to README.md -> install html

- - - - -
9f96317e by Andreas Tille at 2018-03-29T10:25:03+02:00
hardening=+all

- - - - -
9252fd49 by Andreas Tille at 2018-03-29T10:44:47+02:00
Upload to unstable

- - - - -


23 changed files:

- + .gitignore
- + .travis.yml
- Makefile
- − README
- + README.md
- debian/changelog
- + debian/clean
- debian/compat
- debian/control
- debian/docs
- debian/rules
- scripts/dwgsim_mut_to_vcf.pl
- scripts/dwgsim_pileup_eval.pl
- src/dwgsim.c
- src/dwgsim_eval.c
- src/dwgsim_opt.c
- src/mut.c
- + testdata/ex1.test.bfast.fastq.gz
- + testdata/ex1.test.bwa.read1.fastq.gz
- + testdata/ex1.test.bwa.read2.fastq.gz
- + testdata/ex1.test.mutations.txt.gz
- + testdata/ex1.test.mutations.vcf.gz
- + testdata/test.sh


Changes:

=====================================
.gitignore
=====================================
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+src/*.o
+dwgsim*


=====================================
.travis.yml
=====================================
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,6 @@
+sudo: required
+language: c
+compiler:
+- clang
+- gcc
+script: make && make test


=====================================
Makefile
=====================================
--- a/Makefile
+++ b/Makefile
@@ -71,3 +71,7 @@ dist:clean
 	tar -vcf dwgsim-${PACKAGE_VERSION}.tar dwgsim-${PACKAGE_VERSION}; \
 	gzip -9 dwgsim-${PACKAGE_VERSION}.tar; \
 	rm -rv dwgsim-${PACKAGE_VERSION};
+
+test:
+	if [ -d tmp ]; then rm -r tmp; fi
+	/bin/bash testdata/test.sh


=====================================
README deleted
=====================================
--- a/README
+++ /dev/null
@@ -1,9 +0,0 @@
-Welcome to DWGSIM.
-
-Please see the file LICENSE for details.
-Please see the file INSTALL for installation instructions;
-
-This software package has limited support since it is in active development.
-
-Please see the following page for more details:
-https://sourceforge.net/apps/mediawiki/dnaa/index.php?title=Whole_Genome_Simulation


=====================================
README.md
=====================================
--- /dev/null
+++ b/README.md
@@ -0,0 +1,15 @@
+[![Build Status](https://travis-ci.org/nh13/DWGSIM.svg?branch=master)](https://travis-ci.org/nh13/DWGSIM)
+
+Welcome to DWGSIM.
+
+Please see the file LICENSE for details.
+Please see the file INSTALL for installation instructions;
+
+This software package has limited support since it is not longer in active development.
+
+Please use the following command when cloning this repository:
+
+```git clone --recursive``` 
+
+Please see the following page for more details:
+https://github.com/nh13/DWGSIM/wiki


=====================================
debian/changelog
=====================================
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,15 @@
+dwgsim (0.1.12-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream version
+  * cme fix dpkg-control
+  * debhelper 11
+  * build time test needs example data from samtools
+  * Upstream converted README to README.md -> install html
+  * hardening=+all
+
+ -- Andreas Tille <tille at debian.org>  Thu, 29 Mar 2018 10:25:11 +0200
+
 dwgsim (0.1.11-3) unstable; urgency=medium
 
   [Andreas Tille]


=====================================
debian/clean
=====================================
--- /dev/null
+++ b/debian/clean
@@ -0,0 +1 @@
+README.html


=====================================
debian/compat
=====================================
--- a/debian/compat
+++ b/debian/compat
@@ -1 +1 @@
-9
+11
\ No newline at end of file


=====================================
debian/control
=====================================
--- a/debian/control
+++ b/debian/control
@@ -1,14 +1,16 @@
 Source: dwgsim
-Section: science
-Priority: optional
 Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
 Uploaders: Kevin Murray <spam at kdmurray.id.au>
-Build-Depends: debhelper (>= 9),
+Section: science
+Priority: optional
+Build-Depends: debhelper (>= 11~),
+               samtools (>= 1.7-2),
+               python-markdown,
                zlib1g-dev
-Standards-Version: 3.9.6
+Standards-Version: 4.1.3
+Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/dwgsim.git
+Vcs-Git: https://anonscm.debian.org/git/debian-med/dwgsim.git
 Homepage: https://github.com/nh13/DWGSIM/
-Vcs-Git: git://anonscm.debian.org/debian-med/dwgsim.git
-Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/dwgsim.git/
 
 Package: dwgsim
 Architecture: any


=====================================
debian/docs
=====================================
--- a/debian/docs
+++ b/debian/docs
@@ -1 +1 @@
-README
+README*


=====================================
debian/rules
=====================================
--- a/debian/rules
+++ b/debian/rules
@@ -4,5 +4,15 @@
 DPKG_EXPORT_BUILDFLAGS = 1
 include /usr/share/dpkg/default.mk
 
+export DEB_BUILD_MAINT_OPTIONS=hardening=+all
+
 %:
-	dh $@ --parallel
+	dh $@
+
+override_dh_auto_build:
+	dh_auto_build
+	markdown_py -f README.html README.md
+
+override_dh_auto_test:
+	ln -s /usr/share/doc/samtools/examples samtools
+	dh_auto_test


=====================================
scripts/dwgsim_mut_to_vcf.pl
=====================================
--- a/scripts/dwgsim_mut_to_vcf.pl
+++ b/scripts/dwgsim_mut_to_vcf.pl
@@ -1,4 +1,4 @@
-#!/bin/perl
+#!/usr/bin/env perl
 
 use strict;
 use warnings;


=====================================
scripts/dwgsim_pileup_eval.pl
=====================================
--- a/scripts/dwgsim_pileup_eval.pl
+++ b/scripts/dwgsim_pileup_eval.pl
@@ -1,9 +1,9 @@
-#!/usr/bin/perl -w
+#!/usr/bin/env perl -w
 
 # Copied from SAMtools (http://samtools.sourceforge.net)
 
-# This perl code is very cryptic.  Evaluates the fraction 
-# of reads correctly aligned at various mapping quality 
+# This perl code is very cryptic.  Evaluates the fraction
+# of reads correctly aligned at various mapping quality
 # thresholds.
 
 use strict;
@@ -79,7 +79,7 @@ sub dwgsim_eval {
 	print STDERR "\r$ctr\n";
 
 	printf(STDERR "%d\t%d\t%d\t%d\t%d\t%d\t%d\n",
-		$m, $found_correct_1, $found_correct_2, 
+		$m, $found_correct_1, $found_correct_2,
 		$is_correct_1, $is_correct_2,
 		$found_pair_correct,
 		$is_pair_correct);
@@ -167,12 +167,12 @@ sub process_reads {
 				if(1 == $end) { # first read
 					if(abs($pos_1 - $left) <= $gap) {
 						$found_correct_1=1;
-						$correct_index_1 = $i; 
+						$correct_index_1 = $i;
 					}
 				} else { # second read
 					if(abs($pos_2 - $left) <= $gap) {
 						$found_correct_2=1;
-						$correct_index_1 = $i; 
+						$correct_index_1 = $i;
 					}
 				}
 			}


=====================================
src/dwgsim.c
=====================================
--- a/src/dwgsim.c
+++ b/src/dwgsim.c
@@ -50,6 +50,8 @@
 #include "dwgsim.h"
 //#include <config.h>
 
+#define QUAL_MAX 40
+
 uint8_t nst_nt4_table[256] = {
     4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
     4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
@@ -95,18 +97,18 @@ uint8_t nst_nt4_table[256] = {
             n_indel_first[x]++;							\
             if(1 == mut_get_ins(currseq, i, &n, &ins)) { \
                 if(0 == strand[x]) { \
+                    if(k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                     while(n > 0 && k < s[x]) { \
                         tmp_seq[x][k++] = ins & 0x3;                \
                         --n, ins >>= 2; \
                     } \
-                    if(k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                 } else { \
-                    tmp_seq[x][k++] = c & 0xf;						\
                     while(n > 0 && k < s[x]) { \
                         ext_coor[x]++; \
                         tmp_seq[x][k++] = (ins >> ((n-1) << 1) & 0x3);                \
                         --n; \
                     } \
+                    if (k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                 } \
             } else { \
                 int32_t byte_index, bit_index; \
@@ -114,7 +116,8 @@ uint8_t nst_nt4_table[256] = {
                 uint8_t *insertion = NULL; \
                 insertion = mut_get_ins_long_n(currseq->ins[ins], &num_ins); \
                 if(0 == strand[x]) { \
-                    byte_index = mut_packed_len(num_ins) - 1; bit_index = 3 - (num_ins & 3); \
+                    byte_index = mut_packed_len(num_ins) - 1; bit_index = (num_ins-1) & 0x3 ; \
+                    if(k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                     while(num_ins > 0 && k < s[x]) { \
                         assert(0 <= byte_index); \
                         tmp_seq[x][k++] = (insertion[byte_index] >> (bit_index << 1)) & 0x3;                \
@@ -125,9 +128,7 @@ uint8_t nst_nt4_table[256] = {
                             byte_index--; \
                         } \
                     } \
-                    if(k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                 } else { \
-                    tmp_seq[x][k++] = c & 0xf;						\
                     byte_index = 0; bit_index = 0; \
                     while(num_ins > 0 && k < s[x]) { \
                         ext_coor[x]++; \
@@ -139,6 +140,7 @@ uint8_t nst_nt4_table[256] = {
                             byte_index++; \
                         } \
                     } \
+                    if (k < s[x]) tmp_seq[x][k++] = c & 0xf;						\
                 } \
             }													\
         } \
@@ -419,7 +421,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
 {
   seq_t seq;
   mutseq_t *mutseq[2]={NULL,NULL};
-  uint64_t tot_len, ii=0, ctr=0;
+  uint64_t tot_len, ii=0, ctr=0, rand_ii=0;
   int i, l, m, n_ref, contig_i;
   char name[1024], *qstr;
   int32_t name_len_max=0;
@@ -486,7 +488,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
           }
       }
   }
-  fprintf(stderr, "[dwgsim_core] %d sequences, total length: %llu\n", n_ref, (long long)tot_len);
+  fprintf(stderr, "[dwgsim_core] %d sequences, total length: %llu\n", n_ref, (unsigned long long)tot_len);
   rewind(opt->fp_fa);
   mut_print_header_post(opt->fp_vcf);
 
@@ -841,12 +843,16 @@ void dwgsim_core(dwgsim_opt_t * opt)
                       }
                       else {
                           for (i = 0; i < s[j]; ++i) {
-                              qstr[i] = (int)(-10.0 * log(e[j]->start + e[j]->by*i) / log(10.0) + 0.499) + 33;
+                              if (e[j]->start+e[j]->by*i>0) {
+                                  qstr[i] = (int)(-10.0 * log(e[j]->start + e[j]->by*i) / log(10.0) + 0.499) + '!';
+                              } else {
+                                  qstr[i] = QUAL_MAX + '!';
+                              }
                               if(0 < opt->quality_std) {
                                   qstr[i] += (int)((ran_normal() * opt->quality_std) + 0.5);
-                                  if(qstr[i] < 33) qstr[i] = 33;
-                                  if(73 < qstr[i]) qstr[i] = 73;
                               }
+                              if(qstr[i] < '!') qstr[i] = '!';
+                              if(QUAL_MAX + '!' < qstr[i]) qstr[i] = QUAL_MAX + '!';
                           }
                       }
                       qstr[i] = 0;
@@ -859,7 +865,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
                                   name, ext_coor[0]+1, ext_coor[1]+1, strand[0], strand[1], 0, 0,
                                   n_err[0], n_sub[0], n_indel[0],
                                   n_err[1], n_sub[1],n_indel[1],
-                                  (long long)ii, j+1);
+                                  (unsigned long long)ii, j+1);
                           for (i = 0; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], fpo);
                           fprintf(fpo, "\n+\n%s\n", qstr);
@@ -877,7 +883,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
                                   name, ext_coor[0]+1, ext_coor[1]+1, strand[0], strand[1], 0, 0,
                                   n_err[0] - n_err_first[0], n_sub[0] - n_sub_first[0], n_indel[0] - n_indel_first[0], 
                                   n_err[1] - n_err_first[1], n_sub[1] - n_sub_first[1], n_indel[1] - n_indel_first[1],
-                                  (long long)ii, 2 - j);
+                                  (unsigned long long)ii, 2 - j);
                           //fputc('A', fpo);
                           for (i = 1; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], fpo);
@@ -893,7 +899,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
                               (NULL == opt->read_prefix) ? "" : "_",
                               name, ext_coor[0]+1, ext_coor[1]+1, strand[0], strand[1], 0, 0,
                               n_err[0], n_sub[0], n_indel[0], n_err[1], n_sub[1], n_indel[1],
-                              (long long)ii);
+                              (unsigned long long)ii);
                       if(ILLUMINA == opt->data_type || IONTORRENT == opt->data_type) {
                           for (i = 0; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], opt->fp_bfast);
@@ -909,7 +915,6 @@ void dwgsim_core(dwgsim_opt_t * opt)
                           fprintf(opt->fp_bfast, "\n");
                       }
                   }
-                  n_sim++;
               }
               else { // random DNA read
                   for(j=0;j<2;j++) {
@@ -931,12 +936,16 @@ void dwgsim_core(dwgsim_opt_t * opt)
                       }
                       else {
                           for (i = 0; i < s[j]; ++i) {
-                              qstr[i] = (int)(-10.0 * log(e[j]->start + e[j]->by*i) / log(10.0) + 0.499) + 33;
+                              if (e[j]->start+e[j]->by*i>0) {
+                                  qstr[i] = (int)(-10.0 * log(e[j]->start + e[j]->by*i) / log(10.0) + 0.499) + '!';
+                              } else {
+                                  qstr[i] = QUAL_MAX + '!';
+                              }
                               if(0 < opt->quality_std) {
                                   qstr[i] += (int)((ran_normal() * opt->quality_std) + 0.5);
-                                  if(qstr[i] < 33) qstr[i] = 33;
-                                  if(73 < qstr[i]) qstr[i] = 73;
                               }
+                              if(qstr[i] < '!') qstr[i] = '!';
+                              if(QUAL_MAX + '!' < qstr[i]) qstr[i] = '!';
                           }
                       }
                       qstr[i] = 0;
@@ -959,7 +968,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
                                   (NULL == opt->read_prefix) ? "" : "_",
                                   "rand", 0, 0, 0, 0, 1, 1,
                                   0, 0, 0, 0, 0, 0,
-                                  (long long)ii,
+                                  (unsigned long long)rand_ii,
                                   j+1);
                           for (i = 0; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], fpo);
@@ -977,7 +986,8 @@ void dwgsim_core(dwgsim_opt_t * opt)
                                   (NULL == opt->read_prefix) ? "" : "_",
                                   "rand", 0, 0, 0, 0, 1, 1,
                                   0, 0, 0, 0, 0, 0,
-                                  (long long)ii, 2 - j);
+                                  (unsigned long long)rand_ii,
+                                  2 - j);
                           //fputc('A', fpo);
                           for (i = 1; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], fpo);
@@ -993,7 +1003,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
                               (NULL == opt->read_prefix) ? "" : "_",
                               "rand", 0, 0, 0, 0, 1, 1,
                               0, 0, 0, 0, 0, 0,
-                              (long long)ii);
+                              (unsigned long long)rand_ii);
                       if(ILLUMINA == opt->data_type || IONTORRENT == opt->data_type) {
                           for (i = 0; i < s[j]; ++i)
                             fputc("ACGTN"[(int)tmp_seq[j][i]], opt->fp_bfast);
@@ -1009,8 +1019,9 @@ void dwgsim_core(dwgsim_opt_t * opt)
                           fprintf(opt->fp_bfast, "\n");
                       }
                   }
-                  n_sim++;
+                  rand_ii++;
               }
+              n_sim++;
           }
           fprintf(stderr, "\r[dwgsim_core] %llu",
                   (unsigned long long int)ctr);


=====================================
src/dwgsim_eval.c
=====================================
--- a/src/dwgsim_eval.c
+++ b/src/dwgsim_eval.c
@@ -614,28 +614,30 @@ dwgsim_eval_counts_print(dwgsim_eval_counts_t *counts, int32_t a, int32_t d, int
 
   // header
   fprintf(stdout, "# thr | the minimum %s threshold\n", (0 == a) ? "mapping quality" : "alignment score");
-  fprintf(stdout, "# mc | the number of reads mapped correctly that should be mapped at the threshold\n");
-  fprintf(stdout, "# mi | the number of reads mapped incorrectly that should be mapped be mapped at the threshold\n");
-          
-  fprintf(stdout, "# mu | the number of reads unmapped that should be mapped be mapped at the threshold\n");
-          
-  fprintf(stdout, "# um | the number of reads mapped that should be unmapped be mapped at the threshold\n");
-  fprintf(stdout, "# uu | the number of reads unmapped that should be unmapped be mapped at the threshold\n");
-  fprintf(stdout, "# mc' + mi' + mu' + um' + uu' | the total number of reads mapped at the threshold\n");
-  fprintf(stdout, "# mc' | the number of reads mapped correctly that should be mapped at or greater than that threshold\n");
-  fprintf(stdout, "# mi' | the number of reads mapped incorrectly that should be mapped be mapped at or greater than that threshold\n");
+
+  fprintf(stdout, "# mc | the number of correctly mapped reads that should be mapped at the threshold\n");
+  fprintf(stdout, "# mi | the number of incorrectly mapped reads that should be mapped at the threshold\n");
+  fprintf(stdout, "# mu | the number of unmapped reads that should be mapped at the threshold\n");
           
-  fprintf(stdout, "# mu' | the number of reads unmapped that should be mapped be mapped at or greater than that threshold\n");
+  fprintf(stdout, "# um | the number of mapped reads that should be unmapped at the threshold\n");
+  fprintf(stdout, "# uu | the number of unmapped reads that should be unmapped at the threshold\n");
+
+  fprintf(stdout, "# mc + mi + mu + um + uu | the total number of reads at the threshold\n");
+
+  fprintf(stdout, "# mc' | the number of correctly mapped reads that should be mapped at or greater than that threshold\n");
+  fprintf(stdout, "# mi' | the number of incorrectly mapped reads that should be mapped at or greater than that threshold\n");
+  fprintf(stdout, "# mu' | the number of unmapped reads that should be mapped at or greater than that threshold\n");
           
-  fprintf(stdout, "# um' | the number of reads mapped that should be unmapped be mapped at or greater than that threshold\n");
-  fprintf(stdout, "# uu' | the number of reads unmapped that should be unmapped be mapped at or greater than that threshold\n");
-  fprintf(stdout, "# mc' + mi' + mu' + um' + uu' | the total number of reads mapped at or greater than the threshold\n");
+  fprintf(stdout, "# um' | the number of mapped reads that should be unmapped at or greater than that threshold\n");
+  fprintf(stdout, "# uu' | the number of unmapped reads that should be unmapped at or greater than that threshold\n");
+
+  fprintf(stdout, "# mc' + mi' + mu' + um' + uu' | the total number of reads at or greater than the threshold\n");
           
-  fprintf(stdout, "# (mc / (mc' + mi' + mu')) | sensitivity: the fraction of reads that should be mapped that are mapped correctly at the threshold\n");
-  fprintf(stdout, "# (mc / mc' + mi') | positive predictive value: the fraction of mapped reads that are mapped correctly at the threshold\n");
+  fprintf(stdout, "# (mc / (mc' + mi' + mu')) | sensitivity: the fraction of mappable reads that are mapped correctly at the threshold\n");
+  fprintf(stdout, "# (mc / (mc' + mi')) | positive predictive value: the fraction of mapped mappable reads that are mapped correctly at the threshold\n");
   fprintf(stdout, "# (um / (um' + uu')) | false discovery rate: the fraction of random reads that are mapped at the threshold\n");
-  fprintf(stdout, "# (mc' / (mc' + mi' + mu')) | sensitivity: the fraction of reads that should be mapped that are mapped correctly at or greater than the threshold\n");
-  fprintf(stdout, "# (mc' / mc' + mi') | positive predictive value: the fraction of mapped reads that are mapped correctly at or greater than the threshold\n");
+  fprintf(stdout, "# (mc' / (mc' + mi' + mu')) | sensitivity: the fraction of mappable reads that are mapped correctly at or greater than the threshold\n");
+  fprintf(stdout, "# (mc' / (mc' + mi')) | positive predictive value: the fraction of mapped mappable reads that are mapped correctly at or greater than the threshold\n");
   fprintf(stdout, "# (um' / (um' + uu')) | false discovery rate: the fraction of random reads that are mapped at or greater than the threshold\n");
 
 


=====================================
src/dwgsim_opt.c
=====================================
--- a/src/dwgsim_opt.c
+++ b/src/dwgsim_opt.c
@@ -163,19 +163,22 @@ static void get_error_rate(const char *str, error_t *e)
 }
 
 int32_t
-dwgsim_opt_is_int(char *optarg)
+dwgsim_opt_is_int(char *optarg, int32_t neg_ok)
 {
   int32_t i;
-  for(i=0;i<strlen(optarg);i++) {
+  int32_t len = strlen(optarg);
+  if (len == 0) return 0;
+  if ('+' != optarg[0] && (neg_ok == 0 || '-' != optarg[0]) && !isdigit(optarg[0])) return 0;
+  for(i=1;i<len;i++) {
       if(!isdigit(optarg[i])) return 0;
   }
   return 1;
 }
 
 int32_t
-dwgsim_atoi(char *optarg, char flag)
+dwgsim_atoi(char *optarg, char flag, int32_t neg_ok)
 {
-  if(0 == dwgsim_opt_is_int(optarg)) {
+  if(0 == dwgsim_opt_is_int(optarg, neg_ok)) {
       fprintf(stderr, "Error: command line option -%c is not a number [%s]\n", flag, optarg);
       exit(1);
   }
@@ -192,22 +195,22 @@ dwgsim_opt_parse(dwgsim_opt_t *opt, int argc, char *argv[])
   while ((c = getopt(argc, argv, "id:s:N:C:1:2:e:E:r:F:R:X:I:c:S:n:y:BHf:z:Mm:b:v:x:P:q:Q:h")) >= 0) {
       switch (c) {
         case 'i': opt->is_inner = 1; break;
-        case 'd': opt->dist = dwgsim_atoi(optarg, 'd'); break;
+        case 'd': opt->dist = dwgsim_atoi(optarg, 'd', 0); break;
         case 's': opt->std_dev = atof(optarg); break;
-        case 'N': opt->N = dwgsim_atoi(optarg, 'N'); opt->C = -1; break;
+        case 'N': opt->N = dwgsim_atoi(optarg, 'N', 1); opt->C = -1; break;
         case 'C': opt->C = atof(optarg); opt->N = -1; break;
-        case '1': opt->length[0] = dwgsim_atoi(optarg, '1'); break;
-        case '2': opt->length[1] = dwgsim_atoi(optarg, '2'); break;
+        case '1': opt->length[0] = dwgsim_atoi(optarg, '1', 0); break;
+        case '2': opt->length[1] = dwgsim_atoi(optarg, '2', 0); break;
         case 'e': get_error_rate(optarg, &opt->e[0]); break;
         case 'E': get_error_rate(optarg, &opt->e[1]); break;
         case 'r': opt->mut_rate = atof(optarg); break;
         case 'F': opt->mut_freq = atof(optarg); break;
         case 'R': opt->indel_frac = atof(optarg); break;
         case 'X': opt->indel_extend = atof(optarg); break;
-        case 'I': opt->indel_min = dwgsim_atoi(optarg, 'I'); break;
-        case 'c': opt->data_type = dwgsim_atoi(optarg, 'c'); break;
-        case 'S': opt->strandedness = dwgsim_atoi(optarg, 'S'); break;
-        case 'n': opt->max_n = dwgsim_atoi(optarg, 'n'); break;
+        case 'I': opt->indel_min = dwgsim_atoi(optarg, 'I', 0); break;
+        case 'c': opt->data_type = dwgsim_atoi(optarg, 'c', 0); break;
+        case 'S': opt->strandedness = dwgsim_atoi(optarg, 'S', 0); break;
+        case 'n': opt->max_n = dwgsim_atoi(optarg, 'n', 0); break;
         case 'y': opt->rand_read = atof(optarg); break;
         case 'f': 
                   if(NULL != opt->flow_order) free(opt->flow_order);
@@ -216,7 +219,7 @@ dwgsim_opt_parse(dwgsim_opt_t *opt, int argc, char *argv[])
         case 'B': opt->use_base_error = 1; break;
         case 'H': opt->is_hap = 1; break;
         case 'h': return 0;
-        case 'z': opt->seed = dwgsim_atoi(optarg, 'n'); break;
+        case 'z': opt->seed = dwgsim_atoi(optarg, 'z', 1); break;
         case 'M': opt->muts_only = 1; break;
         case 'm': free(opt->fn_muts_input); opt->fn_muts_input = strdup(optarg); opt->fn_muts_input_type = MUT_INPUT_TXT; muts_input_type |= 0x1; break;
         case 'b': free(opt->fn_muts_input); opt->fn_muts_input = strdup(optarg); opt->fn_muts_input_type = MUT_INPUT_BED; muts_input_type |= 0x2; break;
@@ -292,7 +295,19 @@ dwgsim_opt_parse(dwgsim_opt_t *opt, int argc, char *argv[])
   if(NULL != opt->read_prefix) {
       fprintf(stderr, "Warning: remember to use the -P option with dwgsim_eval\n");
   }
-  
+
+  if(0 < opt->length[1]){ //paired end / mate pair
+    int s = opt->length[0] + opt->length[1] + opt->is_inner;
+    double r = (s-opt->dist)*1.0/opt->std_dev;
+    if(r>6){
+      fprintf(stderr, "Error: command line option -d is too small for the read length (%d)\n",opt->dist);
+      return 0;
+    }
+    if(r>4){
+      fprintf(stderr, "Warning: command line option -d is small for the read length (%d). Generation speed could be affected.\n",opt->dist);
+    }
+  }
+
   switch(muts_input_type) {
     case 0x0:
     case 0x1:


=====================================
src/mut.c
=====================================
--- a/src/mut.c
+++ b/src/mut.c
@@ -840,7 +840,7 @@ void mut_print(const char *name, const seq_t *seq, mutseq_t *hap1, mutseq_t *hap
                       }
                       if (0 < i) fprintf(fpout_vcf, "\t%c", "ACGTN"[nst_nt4_table[(int)seq->s[i-1]]]);
                       else fprintf(fpout_vcf, "\t.");
-                      fprintf(fpout_vcf, "\t100\tPASS\tAF=1.0;pl=1;mt=DELETE\n"); 
+                      fprintf(fpout_vcf, "\t100\tPASS\tAF=0.5;pl=1;mt=DELETE\n"); 
                   }
                   // NB: convert back 'c'
                   c[0] = nst_nt4_table[(int)seq->s[i]];
@@ -860,7 +860,7 @@ void mut_print(const char *name, const seq_t *seq, mutseq_t *hap1, mutseq_t *hap
                       }
                       if (0 < i) fprintf(fpout_vcf, "\t%c", "ACGTN"[nst_nt4_table[(int)seq->s[i-1]]]);
                       else fprintf(fpout_vcf, "\t.");
-                      fprintf(fpout_vcf, "\t100\tPASS\tAF=1.0;pl=2;mt=DELETE\n"); 
+                      fprintf(fpout_vcf, "\t100\tPASS\tAF=0.5;pl=2;mt=DELETE\n"); 
                   }
                   // NB: convert back 'c'
                   c[0] = nst_nt4_table[(int)seq->s[i]];


=====================================
testdata/ex1.test.bfast.fastq.gz
=====================================
Binary files /dev/null and b/testdata/ex1.test.bfast.fastq.gz differ


=====================================
testdata/ex1.test.bwa.read1.fastq.gz
=====================================
Binary files /dev/null and b/testdata/ex1.test.bwa.read1.fastq.gz differ


=====================================
testdata/ex1.test.bwa.read2.fastq.gz
=====================================
Binary files /dev/null and b/testdata/ex1.test.bwa.read2.fastq.gz differ


=====================================
testdata/ex1.test.mutations.txt.gz
=====================================
Binary files /dev/null and b/testdata/ex1.test.mutations.txt.gz differ


=====================================
testdata/ex1.test.mutations.vcf.gz
=====================================
Binary files /dev/null and b/testdata/ex1.test.mutations.vcf.gz differ


=====================================
testdata/test.sh
=====================================
--- /dev/null
+++ b/testdata/test.sh
@@ -0,0 +1,28 @@
+#!/bin/bash
+
+set -e
+
+TESTDATADIR=$(dirname ${0});
+
+# Decompress the test data
+pushd $TESTDATADIR;
+for FILE in $(ls -1 *gz)
+do 
+	gunzip -c ${FILE} > $(basename ${FILE} .gz);
+done
+popd
+
+mkdir tmp
+
+# Generate the new test data
+./dwgsim -z 13 -N 10000 samtools/examples/ex1.fa ex1.test
+
+# Test the differences
+for FILE in $(ls -1 ex1.test*)
+do 
+	diff -q ${FILE} ${TESTDATADIR}/${FILE}
+done
+
+# Clean up the testdata
+find ${TESTDATADIR} \! -name "*gz" -type f | grep -v sh$ | xargs rm; 
+rm -r tmp;



View it on GitLab: https://salsa.debian.org/med-team/dwgsim/compare/907840c68c965eb9b9a132e93ec32a8770d97b8b...9252fd494a1089f62fea21da74b0e150bbc09fbe

---
View it on GitLab: https://salsa.debian.org/med-team/dwgsim/compare/907840c68c965eb9b9a132e93ec32a8770d97b8b...9252fd494a1089f62fea21da74b0e150bbc09fbe
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.alioth.debian.org/pipermail/debian-med-commit/attachments/20180329/ffb971e8/attachment-0001.html>


More information about the debian-med-commit mailing list