Merge pull request #21 from hasindu2008/dev

Dev

.github/workflows/c-cpp.yml

@@ -6,6 +6,9 @@ on:
   pull_request:
     branches: [ '*' ]
 
+env:
+  ACTIONS_ALLOW_USE_UNSECURE_NODE_VERSION: true
+
 jobs:
   ubuntu_14:
     name: ubuntu_14
@@ -98,9 +101,38 @@ jobs:
       run: make CC=icc -j8
     - name: test
       run: make test
-  os_x_11:
-    name: OSX 11
-    runs-on: macos-11
+  ubuntu_24:
+    name: Ubuntu 24
+    runs-on: ubuntu-24.04
+    steps:
+    - uses: actions/checkout@v2
+    - name: install packages
+      run: sudo apt-get update && sudo apt-get install zlib1g-dev
+    - name: build
+      run: make -j8
+    - name: test
+      run: make test
+  os_x_12:
+    name: OSX 12
+    runs-on: macos-12
+    steps:
+    - uses: actions/checkout@v2
+    - name: build
+      run: make -j8
+    - name: test
+      run: make test
+  os_x_13:
+    name: OSX 13
+    runs-on: macos-13
+    steps:
+    - uses: actions/checkout@v2
+    - name: build
+      run: make -j8
+    - name: test
+      run: make test
+  os_x_14:
+    name: OSX 14
+    runs-on: macos-14
     steps:
     - uses: actions/checkout@v2
     - name: build

.github/workflows/release-simulation.yml

@@ -7,6 +7,9 @@ on:
   pull_request:
     branches: [ dev ]
 
+env:
+  ACTIONS_ALLOW_USE_UNSECURE_NODE_VERSION: true
+
 jobs:
   ubuntu_14:
     name: ubuntu_14

.github/workflows/release.yml

@@ -6,6 +6,9 @@ on:
     tags:
       - "v*"
 
+env:
+  ACTIONS_ALLOW_USE_UNSECURE_NODE_VERSION: true
+
 jobs:
   ubuntu_14:
     name: ubuntu_14

Makefile

@@ -12,6 +12,7 @@ endif
 BINARY = squigulator
 OBJ = $(BUILD_DIR)/main.o \
       $(BUILD_DIR)/model.o \
+	  $(BUILD_DIR)/methmodel.o \
 	  $(BUILD_DIR)/misc.o \
 	  $(BUILD_DIR)/sim.o \
 	  $(BUILD_DIR)/thread.o \
@@ -43,6 +44,9 @@ $(BUILD_DIR)/main.o: src/main.c $(HEADERS)
 $(BUILD_DIR)/model.o: src/model.c src/model.h  src/misc.h
 	$(CC) $(CFLAGS) $(CPPFLAGS) $< -c -o $@
 
+$(BUILD_DIR)/methmodel.o: src/methmodel.c src/model.h
+	$(CC) $(CFLAGS) $(CPPFLAGS) $< -c -o $@
+
 $(BUILD_DIR)/thread.o: src/thread.c
 	$(CC) $(CFLAGS) $(CPPFLAGS) $< -c -o $@
 

README.md

@@ -6,6 +6,7 @@
 
 Reads directly extracted from the reference genome are simulated without any mutations/variants. If you want to have variants in your simulated data, you can first apply a set of variants to the reference using [bcftools](http://www.htslib.org/download/) and use that as the input to the *squigulator*.
 
+Publication: [https://doi.org/10.1101/gr.278730.123](https://genome.cshlp.org/content/34/5/778.full?sid=cd2c8aec-be46-4c9e-885c-8452ac069f64)
 Preprint: [https://www.biorxiv.org/content/10.1101/2023.05.09.539953v1](https://www.biorxiv.org/content/10.1101/2023.05.09.539953v1)<br/>
 SLOW5 ecosystem: [https://hasindu2008.github.io/slow5](https://hasindu2008.github.io/slow5)<br/>
 
@@ -18,7 +19,7 @@ SLOW5 ecosystem: [https://hasindu2008.github.io/slow5](https://hasindu2008.githu
 
 Please cite the following in your publications when using *squigulator*:
 
-> Gamaarachchi, H., Ferguson, J. M., Samarakoon, H., Liyanage, K., & Deveson, I. W. (2023). Squigulator: simulation of nanopore sequencing signal data with tunable noise parameters. bioRxiv, 2023-05.
+> Gamaarachchi, H., Ferguson, J. M., Samarakoon, H., Liyanage, K., & Deveson, I. W. (2024). Simulation of nanopore sequencing signal data with tunable parameters. Genome Research, 34(5), 778-783.
 
 ```
 @article{gamaarachchi2023squigulator,
@@ -42,7 +43,7 @@ After getting the basic *ssssim* implemented in ~8 hours and successfully baseca
 For x86-64 Linux, you can use the precompiled binaries under [releases](https://github.com/hasindu2008/squigulator/releases):
 
 ```
-VERSION=0.3.0
+VERSION=0.3.0-dirty
 wget https://github.com/hasindu2008/squigulator/releases/download/v${VERSION}/squigulator-v${VERSION}-x86_64-linux-binaries.tar.gz
 tar xf squigulator-v${VERSION}-x86_64-linux-binaries.tar.gz  && cd squigulator-v${VERSION}
 ./squigulator --help

docs/dna.md

@@ -1,6 +1,6 @@
 # squigulator evaluation - genomic DNA
 
-**This document is not yet complete. The experiments were conducted using commit [2b9f3e886eea16](https://github.com/hasindu2008/squigulator/commit/2b9f3e886eea1601d9f0b2021bcd303ad48005a8). As this tool is under activate development, the results may change in future versions.**
+**This document is based on very early experiments conducted using commit [2b9f3e886eea16](https://github.com/hasindu2008/squigulator/commit/2b9f3e886eea1601d9f0b2021bcd303ad48005a8) when the tool was under activate development. For more up to date information, please refer to the [squigulator publication](https://genome.cshlp.org/content/34/5/778.full?sid=cd2c8aec-be46-4c9e-885c-8452ac069f64).**
 
 For comparing *squigulator* against real data, reads mapping to the chr22 from a NA12878 dataset sequenced on a PromethION sequencer was used. The extracted dataset had ~135,000 reads whose mean read length was ~10,800 bases (will be referred to as *real*). To simulate reads using *squigulator*, first, a heterozygous NA12878 heterozygous chr22 reference was created by applying NA12878 GIAB variants to chr22 in hg38 reference genome. Then, *squigulator* was executed to generate 135,000 reads at mean read length of 10,800 (will be referred to as *simulated*).  See [Methods](#Methods) for more information.
 
@@ -61,7 +61,7 @@ Then, variant calling was performed using Nanopolish and Clair3 and the ROC curv
 ## Different basecalling models
 
 
-All the above experiments were based on high accuracy basecalling. Given below is the comparison of results from variant calling when different basecalling models are used (Clair3 on left and Nanopolish on right). 
+All the above experiments were based on high accuracy basecalling. Given below is the comparison of results from variant calling when different basecalling models are used (Clair3 on left and Nanopolish on right).
 
 <img src="img/model-var-clair3.png"  width=50% height=50%><img src="img/model-var-nanopolish.png"  width=50% height=50%>
 

docs/man.md

@@ -10,7 +10,7 @@ Basic options in *squigulator* are as below:
 - `-t INT`: Number of threads [default: 8]
 - `-h`:  help
 - `--ideal`: To generate perfect signals with no noise. See example [here](img/ideal.svg).
-- `--full-contigs`: generate a complete raw signal per each contig in the input reference genome (incompatible with `-n` and `-r`).
+- `--full-contigs`: generate a complete raw signal per each contig in the input reference genome or each sequence in the input sequences (incompatible with `-n` and `-r`).
 - `--version`: print version
 - `--verbose INT`:  verbosity level [default: 4]
 
@@ -25,13 +25,14 @@ Advanced options are as below:
 -  `--dwell-mean FLOAT`: Mean of number of signal samples per k-mer/base. This is usually the sampling rate (4000Hz for DNA and 3000Hz for RNA) divided by translocation speed in bases per second (450 for R9.4.1 pore for DNA and 70 for RNA).  [default: 9.0]
 -  `--dwell-std FLOAT`: Standard deviation of number of signal samples per k-mer/base. Increasing this will increase time-domain noise. Setting this to 0 is same as `--ideal-time`. See example [here](img/dwell.svg). [default: 4.0]
 - `--bps INT`: translocation speed in bases per second (incompatible with --dwell-mean).
--  `--kmer-model FILE`: custom nucleotide k-mer model file (format similar to [here](https://github.com/hasindu2008/f5c/blob/master/test/r9-models/r9.4_450bps.nucleotide.6mer.template.model))
 -  `--prefix=yes/no`: generate prefixes such as adaptor (and polya for RNA). [default: no]
 -  `--seed INT`: seed for random generators (if 0, will be autogenerated). Giving the same seed will produce same results. [default: 0]
 - `--paf-ref`:  in paf output, use the reference as the target instead of read (needs -c)
 - `--cdna`: generate cDNA reads (only valid with dna profiles and the reference must a transcriptome, experimental)
-- `--trans-count FILE`: simulate relative abundance using specified 2-column tsv with first column containing transcript name and the second containing the count  (only for direct-rna and cDNA, experimental). You may generate this a test fatq dataset using minimap2, for example, `minimap2 -cx map-ont transcripts.fa  reads.fastq --secondary=no -t20 -uf | cut -f 6 | sort | uniq -c | awk '{print$2"\t"$1}'`.
+- `--trans-count FILE`: simulate relative abundance using specified 2-column tsv with first column containing transcript name and the second containing the count. See the example at [test/sequin_count.tsv](test/sequin_count.tsv)  (only for direct-rna and cDNA, experimental). You may generate this using a dataset using minimap2, for example, `minimap2 -cx map-ont transcripts.fa  reads.fastq --secondary=no -t20 -uf | cut -f 6 | sort | uniq -c | awk '{print$2"\t"$1}'`.
 - `--trans-trunc=yes/no`: simulate transcript truncation (only for direct-rna, experimental) [default: no]
+- `--ont-friendly=yes/no`:   generate fake uuid for readids and add a dummy end_reason [default: no]
+-- `--meth-freq FILE`: simulate CpG methylation using a frequency tsv file. The tsv file should have three columns, chr, 0-based pos, and methylation frequency. See the example at [test/mfreq.tsv](test/mfreq.tsv). (for DNA, experimental)
 
 Developer options (which are not much tested and error handling) are as below:
 
@@ -42,3 +43,5 @@ Developer options (which are not much tested and error handling) are as below:
 -  `--offset-std FLOAT`:         ADC offset standard deviation (see [here](https://hasindu2008.github.io/slow5specs/summary))
 -  `--median-before-mean`:       Median before mean (see [here](https://hasindu2008.github.io/slow5specs/summary))
 -  `--median-before-std`:        Median before standard deviation (see [here](https://hasindu2008.github.io/slow5specs/summary))
+-  `--kmer-model FILE`:          custom nucleotide k-mer model file (format similar to [f5c models](https://github.com/hasindu2008/f5c/blob/master/test/r9-models/r9.4_450bps.nucleotide.6mer.template.model))
+-  `--meth-model FILE`:          custom methylation k-mer model file (format similar to [f5c models](https://github.com/hasindu2008/f5c/blob/master/test/r9-models/r9.4_450bps.cpg.6mer.template.model))

docs/profile.md

@@ -43,16 +43,16 @@ These profile sets the following advanced parameters in squigulator.
 
 | Parameter           | dna-r10-min   | dna-r10-prom   | rna004-min   | rna004-prom   |
 |---------------------|---------------|-----------------|--------------|---------------|
-| digitisation        | 8192          | 2048            | 2048         | 2048          |
+| digitisation        | 8192          | 2048            | 8192         | 2048          |
 | sample-rate         | 4000          | 4000            | 4000         | 4000          |
 | bps                 | 400           | 400             | 130          | 130           |
-| range               | 1536.598389   | 281.345551      | TBD  | 299.432068    |
-| offset-mean         | 13.380569389  | -127.5655735    | TBD5 | -259.421128  |
-| offset-std         | 16.311471649  | 19.377283387665 | TBD | 16.010841823643 |
-| median-before-mean  | 202.154074388 | 189.87607393756 | TBD | 189.87607393756 |
-| median-before-std   | 13.406139242  | 15.788097978713 | TBD | 15.788097978713 |
-| dwell-mean          | 10.0          | 10.0            | TBD         | 31.0          |
-| dwell-std           | 4.0           | 4.0             | TBD          | 4.0           |
+| range               | 1536.598389   | 281.345551      | 1437.976685  | 299.432068    |
+| offset-mean         | 13.380569389  | -127.5655735    | 12.47686423863 | -259.421128  |
+| offset-std         | 16.311471649  | 19.377283387665 | 10.442126577137 | 16.010841823643 |
+| median-before-mean  | 202.154074388 | 189.87607393756 | 205.08496731088 | 189.87607393756 |
+| median-before-std   | 13.406139242  | 15.788097978713 | 8.6671292866233 | 15.788097978713 |
+| dwell-mean          | 10.0          | 10.0            | 31.0         | 31.0          |
+| dwell-std           | 4.0           | 4.0             | 0.0          | 0.0           |
 
 ## Determining parameters for a profile
 

scripts/test.sh

@@ -21,15 +21,18 @@ ex() {
     fi
 }
 
+# DNA R9
 echo "Basic DNA"
 ex ./squigulator test/nCoV-2019.reference.fasta -o a.slow5 -q a.fasta -n 10 --seed 1 --dwell-std 1.0 -r 20000 -t1 || die "Running the tool failed"
 diff -q test/fasta.exp a.fasta || die "diff failed"
 diff -q test/slow5.exp a.slow5 || die "diff failed"
 
+# RNA R9
 echo "Basic RNA"
 ex ./squigulator -x rna-r9-prom test/rnasequin_sequences_2.4.fa -o a.slow5 -q a.fastq -n 10 --seed 1 --prefix=yes --dwell-std 3.0 -t1  || die "Running the tool failed"
 diff -q test/rna_slow5.exp a.slow5 || die "diff failed"
 
+# --ideal
 echo "--ideal"
 ex ./squigulator test/nCoV-2019.reference.fasta -o a.slow5 -n 2 --seed 1 --ideal  -r 20000 -t1 || die "Running the tool failed"
 diff -q test/dna_ideal_slow5.exp a.slow5 || die "diff failed"
@@ -44,25 +47,33 @@ diff -q test/dna_ideal_amp_slow5.exp a.slow5 || die "diff failed"
 ex ./squigulator test/nCoV-2019.reference.fasta -o a.slow5 -n 2 --seed 1 --amp-noise 0.0  -r 20000  --dwell-std 5.0 -t1 || die "Running the tool failed"
 diff -q test/dna_ideal_amp_slow5.exp a.slow5 || die "diff failed"
 
-
+# prefixes for dna
 echo "--prefix=yes"
 ex ./squigulator test/nCoV-2019.reference.fasta -o a.slow5 -n 2 --seed 1 --prefix=yes  -r 20000  --dwell-std 5.0 -t1 || die "Running the tool failed"
 diff -q test/dna_prefix_slow5.exp a.slow5 || die "diff failed"
 
+# prefixes for rna
+echo "--prefix=yes"
+ex ./squigulator -x rna-r9-prom test/rnasequin_sequences_2.4.fa -o a.slow5 -n 2 --seed 1 --dwell-std 3.0 -t1 --prefix=yes || die "Running the tool failed"
+diff -q test/rna_prefixyes_slow5.exp a.slow5 || die "diff failed"
+
 echo "--prefix=no"
 ex ./squigulator -x rna-r9-prom test/rnasequin_sequences_2.4.fa -o a.slow5 -n 2 --seed 1 --dwell-std 3.0 -t1 || die "Running the tool failed"
 diff -q test/rna_prefixno_slow5.exp a.slow5 || die "diff failed"
 
+# full contigs
 echo "--full-contigs"
 ex ./squigulator test/nCoV-2019.reference.fasta -o a.slow5 --seed 1 --full-contigs  --dwell-std 5.0 -t1 || die "Running the tool failed"
 diff -q test/dna_full_contig.exp a.slow5 || die "diff failed"
 
+# r10 and paf and fasta
 echo "r10 PAF out"
 ex ./squigulator -x dna-r10-prom -o a.slow5 -n 1 --seed 1 --dwell-std 4.0 -t1 test/nCoV-2019.reference.fasta -c a.paf -q a.fa
 diff -q test/dna_r10_paf.exp a.slow5 || die "diff failed"
 diff -q test/dna_r10_paf.paf.exp a.paf || die "diff failed"
 diff -q test/dna_r10_paf.fa.exp a.fa || die "diff failed"
 
+# r9 rna paf and sam and fasta outputs
 echo "r9 rna paf out and sam out"
 ex ./squigulator -x rna-r9-prom -o a.slow5 -n 1 --seed 1 --dwell-std 3.0 -t1 -t1 test/rnasequin_sequences_2.4.fa -c a.paf -q a.fa -a a.sam
 diff -q test/rna_paf.exp a.slow5 || die "diff failed"
@@ -82,7 +93,63 @@ ex ./squigulator -x dna-r10-prom -o a.slow5 -n 2 --seed 2 --dwell-std 4.0 -t1 te
 diff -q test/dna_r10_paf-ref.exp a.slow5 || die "diff failed"
 diff -q test/dna_r10_paf-ref.sam.exp a.sam || die "diff failed"
 
-
+# rna004
+echo "rna004 test"
+ex ./squigulator -x rna004-prom -o a.slow5 -n 1 --seed 1 --dwell-std 3.0 -t1 test/rnasequin_sequences_2.4.fa
+diff -q test/rna004.slow5.exp a.slow5 || die "diff failed"
+
+# -r and -f and --amp-noise
+echo "read lens and fold coverage"
+ex ./squigulator -x dna-r10-prom -o a.slow5 -r 20000 -f 1 --seed 2 --amp-noise 0.5 -t1 test/nCoV-2019.reference.fasta
+diff -q test/dna_r10_amp_noise.exp a.slow5 || die "diff failed"
+
+# dwell mean and std
+ex ./squigulator -x rna004-min -o a.slow5 -n 1 --seed 1 --dwell-mean 30 --dwell-std 3.0 -t1 test/rnasequin_sequences_2.4.fa
+diff -q test/rna004_dwell.exp a.slow5 || die "diff failed"
+
+# bps
+ex ./squigulator -x dna-r10-prom -o a.slow5 --seed 1 --bps 200 -t1 -n 2 test/nCoV-2019.reference.fasta
+diff -q test/bps.exp a.slow5  || die "diff failed"
+
+# cdna
+ex ./squigulator -x dna-r10-min -o a.slow5 -n 1 --seed 1 --dwell-std 3.0 -t1 test/rnasequin_sequences_2.4.fa --cdna
+diff -q test/cdna.exp a.slow5 || die "diff failed"
+
+# trans count
+ex ./squigulator -x rna004-prom -o a.slow5 -n 3 --seed 3 --trans-count test/sequin_count.tsv -t1 test/rnasequin_sequences_2.4.fa
+diff -q test/trans_count.exp a.slow5 || die "diff failed"
+
+# trans count cdna
+ex ./squigulator -x dna-r10-min -o a.slow5 -n 3 --seed 3 --trans-count test/sequin_count.tsv -t1 test/rnasequin_sequences_2.4.fa --cdna
+diff -q test/trans_count_cdna.exp a.slow5 || die "diff failed"
+
+# trans trunc
+ex ./squigulator -x rna004-prom -o a.slow5 -n 1 --seed 1 --trans-trunc -t1 test/rnasequin_sequences_2.4.fa
+diff -q test/trans_trunc.exp a.slow5 || die "diff failed"
+
+# ont-friendly
+ex ./squigulator -x dna-r10-min -o a.slow5 -n 1 --seed 1 -t1 test/rnasequin_sequences_2.4.fa --ont-friendly=yes
+diff -q test/ont_friendly.exp a.slow5 || die "diff failed"
+
+# digitisation, sample-rate, range, offset-mean, offset-std, median-before-mean, median-before-std
+ex ./squigulator -x dna-r10-min -o a.slow5 -n 1 --seed 1 -t1 test/rnasequin_sequences_2.4.fa --digitisation 4096 --sample-rate 10000 --range 300 --offset-mean -1000 --offset-std 0 --median-before-mean 100 --median-before-std 0
+diff -q test/dev.exp a.slow5 || die "diff failed"
+
+#meth r9
+ex ./squigulator -x dna-r9-prom -o a.slow5 --seed 1 -t1 -n 2 -r 29000 test/nCoV-2019.reference.fasta --meth-freq test/mfreq.tsv
+diff -q test/r9_mfreq.exp a.slow5  || die "diff failed"
+# /install/buttery-eel-0.3.1+6.5.7/scripts/eel -i a.slow5 -o a.fastq --config  dna_r9.4.1_450bps_sup.cfg -x cuda:all
+# minimap2 -ax map-ont /genome/nCoV-2019.reference.fasta  a.fastq --secondary=no | samtools sort - -o a.bam && samtools index a.bam
+# f5c index a.fastq --slow5 a.slow5 && f5c call-methylation -r a.fastq  -g test/nCoV-2019.reference.fasta -b a.bam --slow5 a.slow5 -o meth.tsv && f5c meth-freq -i meth.tsv -s -o meth-freq.tsv
+
+#meth r10
+# ex ./squigulator -x dna-r10-prom -o a.slow5 --seed 1 -t1 -n 2 -r 29000 test/nCoV-2019.reference.fasta --meth-freq test/methfreq.tsv
+# diff -q test/r10_methfreq.exp a.slow5  || die "diff failed"
+# /install/buttery-eel-0.3.1+6.5.7/scripts/eel -i a.slow5 -o a.fastq --config  dna_r10.4.1_e8.2_400bps_5khz_sup.cfg -x cuda:all
+# minimap2 -ax map-ont /genome/nCoV-2019.reference.fasta  a.fastq --secondary=no | samtools sort - -o a.bam && samtools index a.bam
+# f5c index a.fastq --slow5 a.slow5 && f5c call-methylation -r a.fastq  -g test/nCoV-2019.reference.fasta -b a.bam --slow5 a.slow5 -o meth.tsv && f5c meth-freq -i meth.tsv -s -o meth-freq.tsv
+
+# threads and batch size
 redundancy_check () {
     N=$(grep -v ^[@#] a.slow5 | cut -f ${1}  | sort | uniq -c | sort -nr -k1,1 | head -1 | awk '{print $1}')
     [ "$N" != "1" ] && die "failed thread test for column ${1}"

slow5lib/Makefile

@@ -55,10 +55,10 @@ $(SHAREDLIB): $(OBJ) $(SVBLIB)
 $(SVBLIB):
 	make -C $(SVB) no_simd=$(no_simd) libstreamvbyte.a
 
-$(BUILD_DIR)/slow5.o: src/slow5.c src/slow5_extra.h src/slow5_idx.h src/slow5_misc.h src/klib/ksort.h $(SLOW5_H)
+$(BUILD_DIR)/slow5.o: src/slow5.c src/slow5_extra.h src/slow5_idx.h src/slow5_misc.h src/klib/ksort.h src/slow5_byte.h $(SLOW5_H)
 	$(CC) $(CFLAGS) $(CPPFLAGS) $< -c -fpic -o $@
 
-$(BUILD_DIR)/slow5_idx.o: src/slow5_idx.c src/slow5_idx.h src/slow5_extra.h src/slow5_misc.h $(SLOW5_H)
+$(BUILD_DIR)/slow5_idx.o: src/slow5_idx.c src/slow5_idx.h src/slow5_extra.h src/slow5_misc.h src/slow5_byte.h $(SLOW5_H)
 	$(CC) $(CFLAGS) $(CPPFLAGS) $< -c -fpic -o $@
 
 $(BUILD_DIR)/slow5_misc.o: src/slow5_misc.c src/slow5_misc.h include/slow5/slow5_error.h

slow5lib/include/slow5/slow5.h

@@ -624,6 +624,10 @@ void slow5_set_log_level(enum slow5_log_level_opt log_level);
 //sets a global variable, so not thread safe
 void slow5_set_exit_condition(enum slow5_exit_condition_opt exit_condition);
 
+//experimental
+/* no error messages printed and not exit when a requested read ID is not found in index*/
+// being tested, do not use until added to documentation
+void slow5_set_skip_rid();
 
 //get the list of hdr data keys in sorted order (only the returned pointer must be freed, not the ones inside - subjet to change)
 //len is the numberof elements

slow5lib/include/slow5/slow5_defs.h

@@ -44,7 +44,7 @@ The API documentation is available at https://hasindu2008.github.io/slow5tools/
 */
 
 // library version
-#define SLOW5_LIB_VERSION "1.2.0-beta"
+#define SLOW5_LIB_VERSION "1.2.0"
 
 // maximum file version supported by this library - independent of slow5 library version above
 // if updating change all 4 below