update tests

hasindu2008 · Sep 19, 2024 · 0c773b5 · 0c773b5
1 parent f7bda59
commit 0c773b5
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 11 deletions.
diff --git a/scripts/test.sh b/scripts/test.sh
@@ -143,8 +143,8 @@ diff -q test/r9_mfreq.exp a.slow5  || die "diff failed"
 # 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"
+# 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

diff --git a/src/sim.c b/src/sim.c
@@ -308,6 +308,8 @@ static core_t *init_core(opt_t opt, profile_t p, char *refname, char *output_fil
             }
 
             if(opt.flag & SQ_R10){
+                ERROR("%s","Methylation simulation is not yet supported for R10 data.");
+                exit(EXIT_FAILURE);
                 INFO("%s","builtin DNA R10 cpg model loaded");
                 kmer_size_meth=set_model(core->cpgmodel, MODEL_ID_DNA_R10_CPG);
                 WARNING("%s","Methylation simulation for R10 is still crude. If you have good control data, please share!");

diff --git a/test/test_meth.sh b/test/test_meth.sh
@@ -58,7 +58,7 @@ RUN_TEST(){
     tail -n+2 methcomp.tsv | datamash ppearson 3:5 > a.acc || die "pearson failed"
     # ~/hasindu2008.git/f5c/scripts/plot_methylation.R  -i methcomp.tsv -o methcomp.pdf
     cat a.acc
-    CHECK_ACC 0.92 a.acc
+    CHECK_ACC 0.93 a.acc
     REMOVE_TMP
 }
 

diff --git a/test/test_misc.sh b/test/test_misc.sh
@@ -37,26 +37,26 @@ CHECK_ACC(){
 
 
 
-samtools  view PNXRXX240011_reads_500k.bam | cut -f 3 | sort | uniq -c | awk '{print $2"\t"$1}' | sort -k1,1 > count.tsv || die "failed extracting chr, pos, meth_freq"
+samtools  view -F 4 ${REF_TRANS_BAM} | cut -f 3 | sort | uniq -c | awk '{print $2"\t"$1}' | sort -k1,1 > count.tsv || die "failed extracting chr, pos, meth_freq"
 
 RUN_REST(){
-    minimap2 -cx map-ont -y -Y --secondary=no ${REF_GENCODE} new.fastq > a.paf  2>> a.log
-    cat a.paf | cut -f 6 | sort | uniq -c | awk '{print $2"\t"$1}' -k1,1 > count_new.tsv || die "failed getting counts"
+    minimap2 -cx map-ont --secondary=no  ${REF_GENCODE} new.fastq > a.paf  2>> a.log
+    cat a.paf | cut -f 6 | sort -k1,1 | uniq -c | awk '{print $2"\t"$1}' > count_new.tsv || die "failed getting counts"
     join count.tsv count_new.tsv > count_joined.tsv || die "failed joining counts"
-    cat count_joined.tsv | datamash ppearson 2:3 > a.acc || die "pearson failed"
+    cat count_joined.tsv | awk '{print $1"\t"$2"\t"$3}' | datamash ppearson 2:3 > a.acc || die "pearson failed"
     cat a.acc
     CHECK_ACC 0.95 a.acc
     REMOVE_TMP
 }
 
 PROF=rna004-prom
-MODEL=rna_rp4_130bps_hac_prom.cfg
-./squigulator ${REF_GENCODE} -x ${PROF}  -t 20 -o new.blow5 --trans-count count.tsv 2> a.log || die "squigulator failed"
-eel  -i new.blow5 --config ${MODEL} --device cuda:all -o new.fastq  &>> a.log || die "eel failed"
+MODEL=rna_rp4_130bps_hac.cfg
+./squigulator ${REF_GENCODE} -x ${PROF}  -t 20 -n 500000 -o new.blow5 --trans-count count.tsv 2> a.log || die "squigulator failed"
+/install/buttery-eel-0.4.2+dorado7.2.13/scripts/eel  -i new.blow5 --config ${MODEL} --device cuda:all -o new.fastq  &>> a.log || die "eel failed"
 RUN_REST
 
 PROF=dna-r9-prom
 MODEL=dna_r9.4.1_450bps_hac_prom.cfg
-./squigulator ${REF_GENCODE} -x ${PROF} -t 20 -o new.blow5 --trans-count count.tsv --cdna 2> a.log || die "squigulator failed"
+./squigulator ${REF_GENCODE} -x ${PROF} -t 20 -n 500000 -o new.blow5 --trans-count count.tsv --cdna 2> a.log || die "squigulator failed"
 eel  -i new.blow5 --config ${MODEL} --device cuda:all -o new.fastq  &>> a.log || die "eel failed"
 RUN_REST