[med-svn] [Git][med-team/dwgsim][upstream] New upstream version 0.1.14
Mohd Bilal (@rmb)
gitlab at salsa.debian.org
Fri Aug 12 16:06:12 BST 2022
Mohd Bilal pushed to branch upstream at Debian Med / dwgsim
Commits:
ec016fec by Mohammed Bilal at 2022-08-12T14:12:49+00:00
New upstream version 0.1.14
- - - - -
11 changed files:
- + .github/workflows/testing.yml
- − .travis.yml
- Makefile
- README.md
- + docs/01_Introduction.md
- + docs/02_Installation.md
- + docs/03_Simulating_Reads.md
- + docs/04_Evaluating_Mappings.md
- src/dwgsim.c
- src/dwgsim_eval.c
- src/regions_bed.c
Changes:
=====================================
.github/workflows/testing.yml
=====================================
@@ -0,0 +1,22 @@
+name: Checks
+
+on:
+ push:
+ branches: [ main ]
+ pull_request:
+ branches: [ main ]
+
+jobs:
+ test:
+ name: Checks
+ runs-on: ubuntu-latest
+ steps:
+ - name: Checkout sources
+ uses: actions/checkout at v2
+ with:
+ submodules: 'recursive'
+ - name: Build
+ run: make
+ - name: Test
+ run: make test
+
=====================================
.travis.yml deleted
=====================================
@@ -1,6 +0,0 @@
-sudo: required
-language: c
-compiler:
-- clang
-- gcc
-script: make && make test
=====================================
Makefile
=====================================
@@ -1,4 +1,4 @@
-PACKAGE_VERSION="0.1.13"
+PACKAGE_VERSION="0.1.14"
CC= gcc
CFLAGS= -g -Wall -O3 #-m64 #-arch ppc
DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_USE_KNETFILE -DPACKAGE_VERSION=\\\"${PACKAGE_VERSION}\\\"
=====================================
README.md
=====================================
@@ -1,15 +1,10 @@
-[](https://travis-ci.org/nh13/DWGSIM)
+[](https://github.com/nh13/DWGSIM/actions/workflows/testing.yml)
+[](http://bioconda.github.io/recipes/dwgsim/README.html)
+[](https://github.com/nh13/dwgsim/blob/main/LICENSE)
+[](https://en.wikipedia.org/wiki/C_(programming_language))
-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```
+Welcome to DWGSIM.
+----
-Please see the following page for more details:
-https://github.com/nh13/DWGSIM/wiki
+Documentation can be found in the [docs folder](docs/01_Introduction.md)
=====================================
docs/01_Introduction.md
=====================================
@@ -0,0 +1,5 @@
+# Introduction
+
+* [Installation](02_Installation.md)
+* [Simulating Reads](03_Simulating_Reads.md)
+* [Evaluating Mappings](04_Evaluating_Mappings.md)
=====================================
docs/02_Installation.md
=====================================
@@ -0,0 +1,19 @@
+# Installation
+
+Please use the following command when cloning this repository:
+
+```console
+git clone --recursive git at github.com:nh13/DWGSIM.git
+```
+
+Next, build the `dwgsim` and `dwgsim_eval` binaries:
+
+```console
+make
+```
+
+You can run a simple integration test with:
+
+```console
+make test
+```
=====================================
docs/03_Simulating_Reads.md
=====================================
@@ -0,0 +1,158 @@
+# Simulating Reads
+
+<!---toc start-->
+ * [Overview](#overview)
+ * [Error rates explained](#error-rates-explained)
+ * [Read names explained](#read-names-explained)
+ * [Mate pair or paired end modes](#mate-pair-or-paired-end-modes)
+ * [Output mutations file](#output-mutations-file)
+ * [Output FASTQ files](#output-fastq-files)
+
+<!---toc end-->
+
+## Overview
+
+The `dwgsim` tool simulates reads from a reference genome FASTA.
+It can be used to evaluate both mapping and variant calling.
+See the `dwgsim_eval` tool for [Evaluating Mappings](04_Evaluating_Mappings.md).
+
+Use `dwgsim -h` to see the full set of command line options (usage).
+
+An example command is as follows:
+
+```console
+dwgsim -N 10000 -1 100 -2 100 -y 0 phix.fasta output
+```
+
+This will simulate 10,000 reads (`-N 10000`), which are paired end 2x100bp (`-1 100` for R1 and `-2 100` for R2), with no random reads (`-y 0`),
+from the given genome FASTA (`phix.fasta`), producing output with prefix `output`.
+
+The following output will be created:
+
+| Name | Description |
+| --- | --- |
+| output.bfast.fastq.gz | Interleaved FASTQ containing both read one and read two |
+| output.bwa.read1.fastq.gz | FASTQ containing only read one |
+| output.bwa.read2.fastq.gz | FASTQ containing only read two |
+| output.mutations.vcf | VCF containing simulated mutations |
+| output.mutations.txt | TXT in a custom format containing simulated mutations (see [below](#output-mutations-file)) |
+
+
+Notes:
+
+- the longest supported insertion is 255bp
+- The `-H` mode will simulate a haploid genome, whereas the default is to simulate a diploid genome.
+
+## Error rates explained
+
+The `-e` and `-E` options accept a uniform error rate (i.e. `-e 0.01` for 1%), or a uniformly increasing/decreasing error rate (i.e. `-e 0.01-0.1` for an error rate of 1% at the start of the read increasing to 10% at the end of the read).
+
+## Read names explained
+
+Read names are of the form:
+
+```
+ @<#1>_<#2>_<#3>_<#4>_<#5>_<#6>_<#7>_<#8>:<#9>:<#10>_<#11>:<#12>:<#13>_<#14>
+```
+
+| Field | Description |
+| --- | --- |
+| 1 | contig name (chromsome name) |
+| 2 | start read 1 (one-based) |
+| 3 | start read 2 (one-based) |
+| 4 | strand read 1 (0 - forward, 1 - reverse) |
+| 5 | strand read 2 (0 - forward, 1 - reverse) |
+| 6 | random read 1 (0 - from the mutated reference, 1 - random) |
+| 7 | random read 2 (0 - from the mutated reference, 1 - random) |
+| 8 | number of sequencing errors read 1 (color errors for colorspace) |
+| 9 | number of SNPs read 1 |
+| 10 | number of indels read 1 |
+| 11 | number of sequencing errors read 2 (color errors for colorspace) |
+| 12 | number of SNPs read 2 |
+| 13 | number of indels read 2 |
+| 14 | read number (unique within a given contig/chromsome) |
+
+Read 1 and read 2 correspond to the first and second reads from a paired-end/mate-pair read respectively.
+
+## Mate pair or paired end modes
+
+This utility can generate mate pair or paired end reads using the `-S` option.
+By default, Illumina (nucleotide) data are paired end, and SOLiD (color space) data are mate pair.
+For clarity, lets call the first end sequence E1 and the second end E2.
+
+Paired end reads have the following orientation:
+
+```
+ 5' E1 -----> .... 3'
+ 3' .... <------- E2 5'
+```
+
+Above, the start co-ordinate of E1 is less than E2, with E1 and E2 reported on opposite strands.
+
+Mate pair reads have following orientation
+
+```
+ 5' E2 -----> .... E1 -------> 3'
+ 3' .... 5'
+```
+
+Above, the start co-ordinate of E1 is greater than E2, with E1 and E2 reported on the same strand.
+
+So for SOLiD mate pair reads, the R3 tag (E2) is listed before the F3 tag (E1).
+For SOLiD paired end reads, the F3 tag (E1) is listed before the F5 tag (E2).
+
+## Output mutations file
+
+The locations of introduced mutations are given in a `<prefix>.mutations.txt` text file.
+There are file columns:
+
+1. the chromosome/contig name
+2. the one-based position
+3. the original reference base
+4. the new reference base(s)
+5. the variant strand(s)
+
+SNPs are represented on one line, and in the case of heterozygous mutations, the new reference base is an [IUPAC](http://www.bioinformatics.org/sms/iupac.html) code.
+
+```
+ contig4 4 T K 1
+```
+
+The above shows a heterozygous mutation at position 4 of contig4 on the first strand, mutating the T base to a heterozygous K (G or T) SNP.
+
+Insertions are represented on one line, where the reference base is missing (indicated by a '-' in the third column).
+
+```
+ contig5 13 - TAC 3
+```
+
+The above shows a homozygous insertion of TAC prior to position 13.
+
+Each base of a deletion is represented on one line, where the new reference base is missing and represented by a '-'.
+
+```
+ contig6 22 A - 2
+```
+
+The above shows a heterozygous deletion of T at position 22 on the second strand.
+Multi-base deletions are show on consecutive lines.
+
+```
+ contig6 22 A - 2
+ contig6 23 C - 2
+```
+
+The above shows a two base homozygous deletion of positions 22 and 23 on the second strand.
+
+## Output FASTQ files
+
+Three FASTQ files are produced, for use with BFAST (interleaved FASTQ) and BWA (one FASTQ per read end).
+
+The FASTQ for BFAST is formatted so that the multi-end reads (paired end or mate pair) occur consecutively in the FASTQ (interleaved), with the read that is 5' of the other listed first.
+For paired end reads, this means that E1 is always listed before E2.
+For mate pair reads, this means that E2 is always listed before E1.
+
+The FASTQs for BWA are split into two files, the first file for one end, the second file for the other, with the read that is 5' of the other in the first file.
+For paired end reads, this means that E1 is in the first file and E2 is in the second file.
+For mate pair reads, this means that E2 is in the first file and E1 is in the second file.
+
=====================================
docs/04_Evaluating_Mappings.md
=====================================
@@ -0,0 +1,62 @@
+# Evaluating Mappings
+
+<!---toc start-->
+ * [Overview](#overview)
+ * [Key Options](#key-options)
+ * [Output](#output)
+
+<!---toc end-->
+
+## Overview
+
+The `dwgsim_eval` tool evaluates the mappings from reads produced by `dwgsim`.
+See the `dwgsim` tool for [Simulating Reads](03_Simulating_Reads.md).
+
+Use `dwgsim_eval -h` to see the full set of command line options (usage)..
+
+An example command is as follows:
+
+```console
+dwgsim_eval mapped.bam
+```
+
+A read is correctly mapped if the mapped start position is within `INT` bases (`-g`) of the simulated read's start position.
+Otherwise, it is counted as incorrectly mapped.
+Furthermore, random reads created by `dwgsim` are expected to be unmapped (counted as false mappings otherwise).
+
+## Key Options
+
+- Use the `-c` option if the reads are from SOLID
+- Use the `-b` option if the alignments are from BWA when mapping SOLiD reads
+- Use the `-z` option if the input reads are expected to be single end (not paired end)
+- Use the `-s` option if the input is a SAM file (default: BAM)
+- Use the `-n` option to specify the number of reads if only a subset of the input reads were used (affects the denominator)
+
+
+## Output
+
+The `dwgsim_eval` tool outputs a table with the following columns:
+
+| Column | Description |
+| --- | --- |
+| `thr` | The alignment threshold (see `-a`) |
+| `mc` | The number of reads mapped correctly that should be mapped at the threshold |
+| `mi` | The number of reads mapped incorrectly that should be mapped be mapped at the threshold |
+| `mu` | The number of reads unmapped that should be mapped be mapped at the threshold |
+| `um` | The number of reads mapped that should be unmapped be mapped at the threshold |
+| `uu` | The number of reads unmapped that should be unmapped be mapped at the threshold |
+| `mc + mi + mu + um + uu` | The total number of reads that should be unmapped be mapped at the threshold |
+| `mc'` | The number of reads mapped correctly that should be mapped at or greater than the threshold |
+| `mi'` | The number of reads mapped incorrectly that should be mapped be mapped at or greater than the threshold |
+| `mu'` | The number of reads unmapped that should be mapped be mapped at or greater than the threshold |
+| `um'` | The number of reads mapped that should be unmapped be mapped at or greater than the threshold |
+| `uu'` | The number of reads unmapped that should be unmapped be mapped at or greater than the threshold |
+| `mc' + mi' + mu' + um' + uu'` | The total number of reads that should be unmapped be mapped at or greater than the threshold |
+| `mc / (mc + mi + mu)` | Sensitivity at the threshold. I.e. the fraction of reads that should be mapped that are mapped correctly. |
+| `mc / (mc + mi)` | Positive predictive value at the threshold. I.e. The fraction of mapped reads that are mapped correctly. |
+| `um / (um + uu)` | False discovery rate at the threshold. I.e. The fraction of random reads that are mapped. |
+| `mc' / (mc' + mi' + mu')` | Sensitivity at or greater than the threshold. I.e. the fraction of reads that should be mapped that are mapped correctly. |
+| `mc' / (mc' + mi')` | Positive predictive value at or greater than the threshold. I.e. The fraction of mapped reads that are mapped correctly. |
+| `um' / (um' + uu')` | False discovery rate at or greater than the threshold. I.e. The fraction of random reads that are mapped. |
+
+"At or greater than the threshold" tells us what our sensitivity, PPV, and FDR would be if we filtered based on that threshold.
=====================================
src/dwgsim.c
=====================================
@@ -571,10 +571,14 @@ void dwgsim_core(dwgsim_opt_t * opt)
case 'G':
case 't':
case 'T':
- break; default:
- num_n++; break;
- } }
- } }
+ break;
+ default:
+ num_n++;
+ break;
+ }
+ }
+ }
+ }
if(0.95 < num_n / (double)l) { // TODO: arbitrary cutoff
fprintf(stderr, "[dwgsim_core] #1 skip sequence '%s' as %d out of %d bases are non-ACGT\n", name, num_n, l);
contig_i++;
@@ -682,7 +686,7 @@ void dwgsim_core(dwgsim_opt_t * opt)
|| pos >= seq.l
|| pos + d - 1 >= seq.l
|| (0 < s[1] && 0 == opt->is_inner && ((0 < s[0] && d <= s[1]) || (d <= s[0] && 0 < s[1])))
- || 0 == regions_bed_query(regions_bed, contig_i, pos, pos + s[0] + s[1] + d - 1));
+ || 0 == regions_bed_query(regions_bed, contig_i, pos, pos + d));
}
// generate the read sequences
=====================================
src/dwgsim_eval.c
=====================================
@@ -27,7 +27,7 @@ print_usage(dwgsim_eval_args_t *args)
fprintf(stderr, "\t\t\t\t\t1: by alignment score\n");
fprintf(stderr, "\t\t\t\t\t2: by suboptimal alignment score\n");
fprintf(stderr, "\t\t\t\t\t3: by alignment score - suboptimal alignment score\n");
- fprintf(stderr, "\t-b\t\talignments are from BWA (SOLiD only) [%s]\n", __IS_TRUE(args->b));
+ fprintf(stderr, "\t-b\t\talignments are from BWA (for SOLiD data only) [%s]\n", __IS_TRUE(args->b));
fprintf(stderr, "\t-c\t\tcolor space alignments [%s]\n", __IS_TRUE(args->c));
fprintf(stderr, "\t-d\tINT\tdivide quality/alignment score by this factor [%d]\n", args->d);
fprintf(stderr, "\t-g\t\tgap \"wiggle\" [%d]\n", args->g);
=====================================
src/regions_bed.c
=====================================
@@ -58,7 +58,7 @@ regions_bed_txt *regions_bed_init(FILE *fp, contigs_t *c)
i++;
}
if(c->n == i) {
- fprintf(stderr, "Error: contig not found [%s]\n", name);
+ fprintf(stderr, "Error: contig not found [%s]. Are you sure your BED is coordinate sorted?\n", name);
exit(1);
}
else if(c->contigs[i].len < start) {
View it on GitLab: https://salsa.debian.org/med-team/dwgsim/-/commit/ec016fec92b831d53f249bc607d90d609f3928f9
--
View it on GitLab: https://salsa.debian.org/med-team/dwgsim/-/commit/ec016fec92b831d53f249bc607d90d609f3928f9
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20220812/af841096/attachment-0001.htm>
More information about the debian-med-commit
mailing list