diff --git a/scripts/test.sh b/scripts/test.sh index 017e91f..66b8d7b 100755 --- 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 index 7467405..466332c 100644 --- 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 index a304663..79cbb45 100755 --- 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 index 3702f5f..ae50539 100755 --- 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 \ No newline at end of file