From 3ceddb8469ca76a7c13fb3c054146e20850660d4 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Fri, 31 Jul 2026 18:37:44 +0000 Subject: [PATCH 1/6] Update ANGSD: implement all other methods, and update test expectations --- tools/angsd/angsd.xml | 776 +++++++++++++++++++++++++++- tools/angsd/{.shed.yml => shed.yml} | 3 +- tools/angsd/test-data/anc.fa | 2 + tools/angsd/test-data/ref.fa | 2 + tools/angsd/test-data/test.bam | Bin 2989 -> 0 bytes tools/angsd/test-data/test.bam.bai | Bin 2760 -> 0 bytes tools/angsd/test-data/test.phe | 48 ++ tools/angsd/test-data/test.sam | 8 + 8 files changed, 821 insertions(+), 18 deletions(-) rename tools/angsd/{.shed.yml => shed.yml} (72%) create mode 100644 tools/angsd/test-data/anc.fa create mode 100644 tools/angsd/test-data/ref.fa delete mode 100644 tools/angsd/test-data/test.bam delete mode 100644 tools/angsd/test-data/test.bam.bai create mode 100644 tools/angsd/test-data/test.phe create mode 100644 tools/angsd/test-data/test.sam diff --git a/tools/angsd/angsd.xml b/tools/angsd/angsd.xml index 73049ba2942..b958171f5e0 100644 --- a/tools/angsd/angsd.xml +++ b/tools/angsd/angsd.xml @@ -1,5 +1,5 @@ - Extract Internal Counts for Contamination Analysis + Comprehensive Population Genetics Analysis (GL, MAF, SAF, Counts, Geno) macros.xml @@ -10,7 +10,7 @@ - - - + + + + ^[\w\d\._:-]+$ + + +
+ + + + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + +
+ + +
+ + + + + + + + +
+ + +
+ + + + + + + + + + + + +
+ + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ - + + + counts_sec['iCounts'] != "0" + + + + + maf_sec['doMaf'] != "0" or (maf_sec['forceMaf'] is True and asso_sec['doAsso'] != "0") + + + + + gl_sec['doGlf'] == "2" + + + gl_sec['doGlf'] == "1" or gl_sec['doGlf'] == "3" or gl_sec['doGlf'] == "4" or gl_sec['doGlf'] == "5" + + + + + saf_sec['doSaf'] != "0" + + + saf_sec['doSaf'] != "0" + + + saf_sec['doSaf'] != "0" + + + + + geno_sec['geno_cond']['enable_geno'] == "yes" and geno_sec['geno_cond']['doGeno'] != 0 + + + geno_sec['geno_cond']['enable_geno'] == 'yes' and geno_sec['geno_cond']['doGeno'] != 0 and geno_sec['geno_cond']['doPlink'] == '2' + + + geno_sec['geno_cond']['enable_geno'] == 'yes' and geno_sec['geno_cond']['doGeno'] != 0 and geno_sec['geno_cond']['doPlink'] == '2' + + + + + counts_sec['dumpCounts'] != "0" + + + counts_sec['dumpCounts'] != "0" and counts_sec['dumpCounts'] != "1" + + + + + counts_sec['doDepth'] != "0" + + + counts_sec['doDepth'] != "0" + + + + + asso_sec['doAsso'] == "1" or asso_sec['doAsso'] == "2" + + + asso_sec['doAsso'] == "4" or asso_sec['doAsso'] == "5" or asso_sec['doAsso'] == "6" + + + + + hwe_sec['doHWE'] != "0" + + + + + other_sec['doHaploCall'] != "0" + + + + + other_sec['doError'] != "0" + + + other_sec['doAncError'] != "0" + + + other_sec['doHetPlas'] is True + + + + + other_sec['doFasta'] != "0" + + + other_sec['doAbbababa'] != "0" + + - - + +
+ + +
+
+ + +
- + + + + + + + + + +
+ + +
+
+ +
+
+ +
+ + + + + + + + + + + + + +
+ +
+
+ + +
+
+ +
+
+ +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +
+
+ + +
+
+ +
+
+ +
+ + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+
+ + +
+ + + + + + + + + + + + + + + + + + + + +
+ +
+
+ +
+
+ +
+
+ + +
+
+ +
+
+ +
+
+ +
+ + + + + + + + + + + + + + + + + + + +
+ +
+
+ +
+
+ + + +
+ + + + + + + + + + + + + + +
+ +
+
+ +
+ + + + diff --git a/tools/angsd/.shed.yml b/tools/angsd/shed.yml similarity index 72% rename from tools/angsd/.shed.yml rename to tools/angsd/shed.yml index 0e964e8e3bd..7850ccb7dae 100644 --- a/tools/angsd/.shed.yml +++ b/tools/angsd/shed.yml @@ -4,8 +4,7 @@ description: ANGSD and nuclear contamination estimation tools for ancient DNA an homepage_url: https://github.com/angsd/angsd long_description: | This repository contains ANGSD for generating internal counts and an X-chromosome - contamination estimation tool, ported from the nf-core/eager pipeline. - It is designed for paleogenomics research to assess sample quality. + contamination estimation tool. It is designed for paleogenomics research to assess sample quality. remote_repository_url: https://github.com/galaxyproject/tools-iuc/tree/main/tools/angsd type: unrestricted categories: diff --git a/tools/angsd/test-data/anc.fa b/tools/angsd/test-data/anc.fa new file mode 100644 index 00000000000..42e33be4663 --- /dev/null +++ b/tools/angsd/test-data/anc.fa @@ -0,0 +1,2 @@ +>1 +AAAAAACCCCGGGGTTTTTTNNNNNNNNNN \ No newline at end of file diff --git a/tools/angsd/test-data/ref.fa b/tools/angsd/test-data/ref.fa new file mode 100644 index 00000000000..f2ce5212e4b --- /dev/null +++ b/tools/angsd/test-data/ref.fa @@ -0,0 +1,2 @@ +>1 +AAAAACCCCCGGGGGTTTTTNNNNNNNNNN \ No newline at end of file diff --git a/tools/angsd/test-data/test.bam b/tools/angsd/test-data/test.bam deleted file mode 100644 index 2ad644a797dcc51f1bf17d2ecac56c2201df1ef3..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2989 zcmV;e3sUqSiwFb&00000{{{d;LjnNt1g%xga-&8NR(2}~`H*Lr+n$gZ1VX5il(87g z#QYhMEvH;HHNePL76}Egoy{BME%F+9lbrGzxuw#Bgec8=PP?`>HZ!)Vc^UFm@x*dUOR?aTZOyde zKR@zy$}3u|X~t?+@tUfX*5%_mg#%}_8@X=~Byh&IL>j{ zf7EX?C4*6>q)H~~SjUpO&}7LN@kG)_9AlmFg^{K`0JoI^-Rq*rD^Z@8bZII-ngsqTNd@)}PjI%vooMdchsJfKu8JARM=#r5tiKHr{ zp*&I+YL=>pd|g(4v>rG4Qq-an6t8(%EsOPd(8Ilk>oEjLQfyAAwl}4L@7R-DWac>- zZ7-O`$YmS)ysRJS*^PNvaD}pUMzKukifuMf^pdlzqRTZsGhMnU$`z|AMDb^_2=k>{ z4Esty-D@#aQdO7LM3Pf2lT;OBiDil=VJ#i<5j2ISHnoSmE{Q{X647(pAzz#wGWYd; z_B7KO(+5}@B#9xZ1~(*;Q_@F#fVBb7k_F#o4~gqv6Kc%@5Y&uGgy8Go;4q_T<`NqD zVN9&nxQEeb;)kLTQ5?c{dW!%}{d3SXBJ^rzI(BSBFpxEf;SJ{{0}bpsfxR8YrZ*u7 z%QC>ONqFW*c6^IS82TYtu4Xog+r-xp&@yA&_e4fWOgt!m1d@5N(9oRIh+NIUN4oW8&JKSJ9c__>MdRH8MZg2OThPQICLB9Hl$y*bTy}fNx#QnZ*SE;AWvlir z6&2W*qeaZ196*cO9-VoCX@|sWxJE~ay=&swK|mTMUI49A#F>G@VQktSR7!475M9~< zNE${l^{B8!FTN{kNA%M3Z@iQlIX&BEQE!&re)D!N`=I+zdQF% zpCI(7A$-6gfldT|7J7AAw2*1&&jRY1E_sb<>MmRmXiDbj+hS8MSIb{{hEkX!k0nRN z0_Ci#QTOfZH=y_ZLgBn{=oCie)F`w zVfXv$p?JPnW*J`}^}`R=CAG-80|CQV4%<6)h-r1K>gV$euJX?&+J3zqYlEY)_O#^B zNP0k0g0{CGHhN-C&Bt5*vZUeWdGd9q(|Hen1$XZ6|9aOfmCnjl`SAgg?L-^kVQM_#wjOSMHXWW~7@cgNnx0)aA6}lnaBgmPX|RZ9+Nb=p zv-3-Xg_*%Sm)@Jb(C;tsUk3fz;p6_|;N8oEi*tj;{-ulkv$IRHi-V;;>-U!yE??yD zINL|P_8fY3dJgsYJ^{Py&-{vpEIhcuS@*({i)-godW+NmSLf zKu(iT*%#-{C2+ocqZk0IoTX`m#EF8Q%6 zYepc)T402Mq&7COauqxeD!9W{&@7XXzk~|*-r*{^0)!x>7xaQjJc^ry0yiczk;+w| za~CLyRt>bts-bmM74WT#AE565=n+6CxQ%hk++zHrJ-$~7;@aID|2R5a0Kd%zK|D-d z-1IOewvFz=dnT0*tZHY~+O!xZfRXQk5$u~?M$07Rups2t??A|NLqc9I8XO&un$&^d z;5sHl_rW3^#Rz>}23i7H(`DHR4Jk>JC{ELkBFW>8r(QW#BRxEaI90fdnU7i5>Ef+s zY+H}-A!Bdc7Z%!yra?j-8)GMXPq0cmDz}i@c+5o0H?>B+XSQya z9!LKWm=EmAF^{9kVsR8hz#vTPf-FABDp|&C%Os zqSu8-?o2#F|5C9N=LG~dDRn*CbQ}}i29H#eD2Vl{6zh?e#MRIT$7mh+5mgf@{*U%} z1{m7G?a`ubiXBu+tkfnVEb1in^V%x;CuuIWB;9FCYsvhfs0V!5blnrr1 zPFdVN!lhxBDUCk`=r010M*Zrs{0$_dMMJqF7`XW%`< zdNN8ULnTO*K=J+RD!}TR0`U0RyJepJga`{KXL9feIYkoz1qSRUfVY_Kpn0A;*Iq|A zF01~{OM}zW`UYqLag|65612Rwf|eII5eJ&K!yro%#?Yr+6Hi*j`GCITTi$oDQuF0* z!NCl0uwVFs74(Ja0(y8#n-~K4(RGzk=jU8s4h;1br~f2HmE22h0#}H`Ai0#D@bV+H4R9_Buxi0j15P5bB$1XyEpp;T8Q$NYIW)QX1(|6d0it=z5SQO8OK9#+pm?X{hJJglFmy zV=sD4Tk^-2ng zlCKq36!#D@p-I6G0Vg?hs3=A9a66~#N8 zF;G-*>Y(kw{Ly7N6BIiF8aqgA6Jtgl@d*Bj-(eWCSmF_von^L3o0zSw4q&`R8q3E; z`RJs7i`#if{{Es}0_cvq@EF3mV`VLppBlh=tUzzy73*U~*GW5f7z@!qv;7EqSkv6Gq2t02PwbH@w6~6oFI}ZMOcJy|^sTp?Q zLLHz0zISrb8D8iqhNfsR6)`$uDw5HVYGL?vW&={Ti@^bDGX}F~+QlV~QfK}F+%i5Q jni~KBABzYC000000RIL6LPG)o8vp|U0000000000XW*&S diff --git a/tools/angsd/test-data/test.bam.bai b/tools/angsd/test-data/test.bam.bai deleted file mode 100644 index c583bb5db0799eda236f1e8668e1c7eb3da2db60..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2760 zcmZ>A^kigZzyusXp$1_F21X!VDfNj3qGBR9i1wPqz`z6+KoysQD1wt!Q2h)H5Or*f z5V{&F4&(DML&R&K;xN8H8$`SoDh}hva6!cDpyDt-qc}vo9x4vw!`yEH@(k21qckL} sM#Et=93UBXG#p05VYGZ4Ee}V_fzfh+e&qmq8xXyXhu(&R#0@4105C2^0ssI2 diff --git a/tools/angsd/test-data/test.phe b/tools/angsd/test-data/test.phe new file mode 100644 index 00000000000..cfbede7d5bb --- /dev/null +++ b/tools/angsd/test-data/test.phe @@ -0,0 +1,48 @@ +1 +1 +1 +1 +1 +0 +0 +0 +1 +0 +1 +0 +1 +1 +1 +0 +1 +0 +1 +0 +0 +1 +0 +1 +0 +1 +0 +0 +0 +0 +0 +0 +0 +0 +1 +0 +0 +0 +0 +1 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/tools/angsd/test-data/test.sam b/tools/angsd/test-data/test.sam new file mode 100644 index 00000000000..e9169938c47 --- /dev/null +++ b/tools/angsd/test-data/test.sam @@ -0,0 +1,8 @@ +@HD VN:1.0 SO:coordinate +@SQ SN:1 LN:30 +@RG ID:RG0 LB:lb1 SM:sm0 PI:450 CN:SC PL:ILLUMINA +rd1sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN ]asdf]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 +rd2sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN 1234]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 +rd3sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN QWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 +rd4sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN tWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 +rd5sm0 0 1 1 60 30M * 0 0 ATGTNACGTNACGTNACGTNACGTNNNNNN DWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 \ No newline at end of file From 8adc8a77d684366f7168239f33335c77de87044f Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Fri, 31 Jul 2026 19:41:55 +0000 Subject: [PATCH 2/6] Fix FASTA indexing issue --- tools/angsd/{shed.yml => .shed.yml} | 0 tools/angsd/angsd.xml | 9 ++++++ tools/angsd/macros.xml | 2 +- tools/angsd/test-data/test.phe | 47 ----------------------------- 4 files changed, 10 insertions(+), 48 deletions(-) rename tools/angsd/{shed.yml => .shed.yml} (100%) diff --git a/tools/angsd/shed.yml b/tools/angsd/.shed.yml similarity index 100% rename from tools/angsd/shed.yml rename to tools/angsd/.shed.yml diff --git a/tools/angsd/angsd.xml b/tools/angsd/angsd.xml index b958171f5e0..d77db686b82 100644 --- a/tools/angsd/angsd.xml +++ b/tools/angsd/angsd.xml @@ -20,6 +20,15 @@ echo 'sample_${i}.bam' >> bam.filelist && #end for + #if $ref_file: + ln -s '$ref_file' reference.fa && + samtools faidx reference.fa && + #end if + #if $anc_file: + ln -s '$anc_file' ancestral.fa && + samtools faidx ancestral.fa && + #end if + ## 2. Run ANGSD angsd -bam bam.filelist diff --git a/tools/angsd/macros.xml b/tools/angsd/macros.xml index dbb50781719..848054e473c 100644 --- a/tools/angsd/macros.xml +++ b/tools/angsd/macros.xml @@ -1,7 +1,7 @@ 0.940 - 0 + 1 25.1 diff --git a/tools/angsd/test-data/test.phe b/tools/angsd/test-data/test.phe index cfbede7d5bb..d00491fd7e5 100644 --- a/tools/angsd/test-data/test.phe +++ b/tools/angsd/test-data/test.phe @@ -1,48 +1 @@ 1 -1 -1 -1 -1 -0 -0 -0 -1 -0 -1 -0 -1 -1 -1 -0 -1 -0 -1 -0 -0 -1 -0 -1 -0 -1 -0 -0 -0 -0 -0 -0 -0 -0 -1 -0 -0 -0 -0 -1 -0 -0 -0 -0 -0 -0 -0 -0 From 0583eb0e390c82aee63fdaa067fadbd4c82a316c Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 9 Aug 2026 09:58:02 +0000 Subject: [PATCH 3/6] Update ANGSD tool: summaries --- tools/angsd/.shed.yml | 22 +- tools/angsd/angsd.xml | 820 ------------------ tools/angsd/angsd_contamination.xml | 79 -- tools/angsd/angsd_summary.xml | 342 ++++++++ tools/angsd/estError.R | 302 +++++++ tools/angsd/macros.xml | 43 + tools/angsd/test-data/output.icnts.gz | Bin 309 -> 0 bytes tools/angsd/test-data/test.bam | Bin 0 -> 2989 bytes tools/angsd/test-data/test.bam.bai | Bin 0 -> 2760 bytes tools/angsd/test-data/test.phe | 1 - tools/angsd/test-data/test.sam | 10 +- tools/angsd/test-data/{anc.fa => test_anc.fa} | 0 tools/angsd/test-data/{ref.fa => test_ref.fa} | 0 13 files changed, 703 insertions(+), 916 deletions(-) delete mode 100644 tools/angsd/angsd.xml delete mode 100644 tools/angsd/angsd_contamination.xml create mode 100644 tools/angsd/angsd_summary.xml create mode 100644 tools/angsd/estError.R delete mode 100644 tools/angsd/test-data/output.icnts.gz create mode 100644 tools/angsd/test-data/test.bam create mode 100644 tools/angsd/test-data/test.bam.bai delete mode 100644 tools/angsd/test-data/test.phe rename tools/angsd/test-data/{anc.fa => test_anc.fa} (100%) rename tools/angsd/test-data/{ref.fa => test_ref.fa} (100%) diff --git a/tools/angsd/.shed.yml b/tools/angsd/.shed.yml index 7850ccb7dae..cd4a4b74d1e 100644 --- a/tools/angsd/.shed.yml +++ b/tools/angsd/.shed.yml @@ -1,12 +1,12 @@ -name: angsd -owner: iuc -description: ANGSD and nuclear contamination estimation tools for ancient DNA analysis. -homepage_url: https://github.com/angsd/angsd -long_description: | - This repository contains ANGSD for generating internal counts and an X-chromosome - contamination estimation tool. It is designed for paleogenomics research to assess sample quality. -remote_repository_url: https://github.com/galaxyproject/tools-iuc/tree/main/tools/angsd -type: unrestricted -categories: - - Sequence Analysis +name: angsd +owner: iuc +description: ANGSD tools for calculating various genomic statistics, contamination, and error estimation from BAM files. +homepage_url: https://github.com/angsd/angsd +long_description: | + This repository contains ANGSD tools for generating summary statistics from BAM files. + It includes modules for X-chromosome nuclear contamination estimation, sequencing error estimation, allele counting, depth distributions, and mismatch matrix calculations. +remote_repository_url: https://github.com/galaxyproject/tools-iuc/tree/main/tools/angsd +type: unrestricted +categories: + - Sequence Analysis - Statistics \ No newline at end of file diff --git a/tools/angsd/angsd.xml b/tools/angsd/angsd.xml deleted file mode 100644 index d77db686b82..00000000000 --- a/tools/angsd/angsd.xml +++ /dev/null @@ -1,820 +0,0 @@ - - Comprehensive Population Genetics Analysis (GL, MAF, SAF, Counts, Geno) - - - macros.xml - - - - - > bam.filelist && - #end for - - #if $ref_file: - ln -s '$ref_file' reference.fa && - samtools faidx reference.fa && - #end if - #if $anc_file: - ln -s '$anc_file' ancestral.fa && - samtools faidx ancestral.fa && - #end if - - ## 2. Run ANGSD - angsd - -bam bam.filelist - -out output - -nThreads "\${GALAXY_SLOTS:-1}" - - #if $ref_file: - -ref '$ref_file' - #end if - #if $anc_file: - -anc '$anc_file' - #end if - #if $sites_file: - -sites '$sites_file' - #end if - #if $region: - -r '$region' - #end if - - ## Basic BAM Filters - -minMapQ $basic_filters.minMapQ - -minQ $basic_filters.minQ - #if $basic_filters.remove_bads: - -remove_bads 1 - #end if - #if $basic_filters.uniqueOnly: - -uniqueOnly 1 - #end if - #if $basic_filters.only_proper_pairs: - -only_proper_pairs 1 - #end if - #if str($basic_filters.trim) != "0": - -trim $basic_filters.trim - #end if - #if str($basic_filters.c_adj) != "0": - -C $basic_filters.c_adj - #end if - #if str($basic_filters.baq) != "0": - -baq $basic_filters.baq - #end if - #if str($basic_filters.minInd) != "0": - -minInd $basic_filters.minInd - #end if - #if str($basic_filters.downSample) != "0.0": - -downSample $basic_filters.downSample - #end if - - ## Genotype Likelihoods - #if str($gl_sec.gl) != "0": - -GL $gl_sec.gl - #end if - #if str($gl_sec.doGlf) != "0": - -doGlf $gl_sec.doGlf - #end if - - ## Counts and Depths - #if str($counts_sec.doCounts) != "0": - -doCounts $counts_sec.doCounts - #if str($counts_sec.dumpCounts) != "0": - -dumpCounts $counts_sec.dumpCounts - #end if - #if str($counts_sec.iCounts) != "0": - -iCounts $counts_sec.iCounts - #end if - #if str($counts_sec.setMaxDepth) != "-1": - -setMaxDepth $counts_sec.setMaxDepth - #end if - #if str($counts_sec.setMinDepth) != "-1": - -setMinDepth $counts_sec.setMinDepth - #end if - #end if - - #if str($counts_sec.doDepth) != "0": - -doDepth $counts_sec.doDepth - -maxDepth $counts_sec.maxDepth - #end if - - ## Association Studies - #if str($asso_sec.doAsso) != "0": - -doAsso $asso_sec.doAsso - #if $asso_sec.yBin: - -yBin '$asso_sec.yBin' - #end if - #if $asso_sec.yCount: - -yCount '$asso_sec.yCount' - #end if - #if $asso_sec.yQuant: - -yQuant '$asso_sec.yQuant' - #end if - #if $asso_sec.cov: - -cov '$asso_sec.cov' - #end if - -model $asso_sec.model - #if $asso_sec.pvalue: - -Pvalue 1 - #end if - #end if - - ## Allele Frequencies & Posteriors - #if str($maf_sec.doMaf) != "0": - -doMaf $maf_sec.doMaf - #if str($maf_sec.doPost) != "0": - -doPost $maf_sec.doPost - #end if - #if str($maf_sec.minMaf) != "-1.0": - -minMaf $maf_sec.minMaf - #end if - #if str($maf_sec.SNP_pval) != "1.0": - -SNP_pval $maf_sec.SNP_pval - #end if - #if str($maf_sec.rmTriallelic) != "0.0": - -rmTriallelic $maf_sec.rmTriallelic - #end if - #if $maf_sec.skipMissing: - -skipMissing 1 - #end if - #end if - #if $maf_sec.forceMaf and (str($maf_sec.doMaf) != "0" or str($asso_sec.doAsso) != "0"): - -forceMaf 1 - #end if - - ## SAF (Site Frequency Spectrum) - #if str($saf_sec.doSaf) != "0": - -doSaf $saf_sec.doSaf - #if $saf_sec.isHap: - -isHap 1 - #end if - #if $saf_sec.noTrans: - -noTrans 1 - #end if - #if $saf_sec.underFlowProtect: - -underFlowProtect 1 - #end if - -scoreTol $saf_sec.scoreTol - #end if - - ## Genotype Calling - #if str($geno_sec.geno_cond.enable_geno) == 'yes': - #if str($geno_sec.geno_cond.doGeno) != "0": - -doGeno $geno_sec.geno_cond.doGeno - #if str($geno_sec.geno_cond.doPlink) != "0": - -doPlink $geno_sec.geno_cond.doPlink - #end if - #end if - -postCutoff $geno_sec.geno_cond.postCutoff - #end if - - ## Hardy-Weinberg Equilibrium - #if str($hwe_sec.doHWE) != "0": - -doHWE $hwe_sec.doHWE - #if str($hwe_sec.minHWEpval) != "-1.0": - -minHWEpval $hwe_sec.minHWEpval - #end if - #if str($hwe_sec.maxHWEpval) != "-1.0": - -maxHWEpval $hwe_sec.maxHWEpval - #end if - #if str($hwe_sec.maxHetFreq) != "-1.0": - -maxHetFreq $hwe_sec.maxHetFreq - #end if - #if str($hwe_sec.minHetFreq) != "-1.0": - -minHetFreq $hwe_sec.minHetFreq - #end if - #end if - - ## Major/Minor Alleles - #if str($majorminor_sec.doMajorMinor) != "0": - -doMajorMinor $majorminor_sec.doMajorMinor - #if $majorminor_sec.rmTrans: - -rmTrans 1 - #end if - #if $majorminor_sec.skipTriallelic: - -skipTriallelic 1 - #end if - #end if - - ## Other - #if str($other_sec.doError) != "0": - -doError $other_sec.doError - #end if - #if str($other_sec.doAncError) != "0": - -doAncError $other_sec.doAncError - #end if - #if str($other_sec.doFasta) != "0": - -doFasta $other_sec.doFasta - #end if - #if str($other_sec.doAbbababa) != "0": - -doAbbababa $other_sec.doAbbababa - #end if - #if str($other_sec.doHaploCall) != "0": - -doHaploCall $other_sec.doHaploCall - #end if - #if $other_sec.doHetPlas: - -doHetPlas 1 - #end if - - ]]> - - - - - - - - ^[\w\d\._:-]+$ - - - -
- - - - - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - - - -
- - -
- - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - -
- - -
- - - - - - - - -
- - -
- - - - - - - - - - - - -
- - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - counts_sec['iCounts'] != "0" - - - - - maf_sec['doMaf'] != "0" or (maf_sec['forceMaf'] is True and asso_sec['doAsso'] != "0") - - - - - gl_sec['doGlf'] == "2" - - - gl_sec['doGlf'] == "1" or gl_sec['doGlf'] == "3" or gl_sec['doGlf'] == "4" or gl_sec['doGlf'] == "5" - - - - - saf_sec['doSaf'] != "0" - - - saf_sec['doSaf'] != "0" - - - saf_sec['doSaf'] != "0" - - - - - geno_sec['geno_cond']['enable_geno'] == "yes" and geno_sec['geno_cond']['doGeno'] != 0 - - - geno_sec['geno_cond']['enable_geno'] == 'yes' and geno_sec['geno_cond']['doGeno'] != 0 and geno_sec['geno_cond']['doPlink'] == '2' - - - geno_sec['geno_cond']['enable_geno'] == 'yes' and geno_sec['geno_cond']['doGeno'] != 0 and geno_sec['geno_cond']['doPlink'] == '2' - - - - - counts_sec['dumpCounts'] != "0" - - - counts_sec['dumpCounts'] != "0" and counts_sec['dumpCounts'] != "1" - - - - - counts_sec['doDepth'] != "0" - - - counts_sec['doDepth'] != "0" - - - - - asso_sec['doAsso'] == "1" or asso_sec['doAsso'] == "2" - - - asso_sec['doAsso'] == "4" or asso_sec['doAsso'] == "5" or asso_sec['doAsso'] == "6" - - - - - hwe_sec['doHWE'] != "0" - - - - - other_sec['doHaploCall'] != "0" - - - - - other_sec['doError'] != "0" - - - other_sec['doAncError'] != "0" - - - other_sec['doHetPlas'] is True - - - - - other_sec['doFasta'] != "0" - - - other_sec['doAbbababa'] != "0" - - - - - - - -
- - -
-
- - -
- - - - - - - - - - - -
- - -
-
- -
-
- -
- - - - - - - - - - - - - - - - -
- -
-
- - -
-
- -
-
- -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - -
-
- - -
-
- -
-
- -
- - - - - - - - - - - - - - - - - - - - - - - - - - -
- -
-
- - -
- - - - - - - - - - - - - - - - - - - - -
- -
-
- -
-
- -
-
- - -
-
- -
-
- -
-
- -
- - - - - - - - - - - - - - - - - - - -
- -
-
- -
-
- - - -
- - - - - - - - - - - - - - -
- -
-
- -
- - - - - - - - - - diff --git a/tools/angsd/angsd_contamination.xml b/tools/angsd/angsd_contamination.xml deleted file mode 100644 index 5159406791f..00000000000 --- a/tools/angsd/angsd_contamination.xml +++ /dev/null @@ -1,79 +0,0 @@ - - Estimate nuclear contamination on X chromosome - - - macros.xml - - - - - - - - - 'contamination_report.out' && - - ## 3. Run the EAGER python parser - python3 '$__tool_directory__/print_x_contamination.py' 'contamination_report.out' - ]]> - - - - - - - - - - - generate_json - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tools/angsd/angsd_summary.xml b/tools/angsd/angsd_summary.xml new file mode 100644 index 00000000000..64ee6334f34 --- /dev/null +++ b/tools/angsd/angsd_summary.xml @@ -0,0 +1,342 @@ + + Calculate statistics from BAM files + + + macros.xml + + + + + + + + + + > bam.filelist && + #end for + + ## Essential command + angsd -bam bam.filelist -out output -nThreads "\${GALAXY_SLOTS:-1}" + + ## Flags for selected operation + #if $summary_type.summary_type_selector == 'contamination': + -doCounts 1 + -iCounts 1 + @BAM_PARAM@ + @COUNTS_PARAM@ && + + ln -s '$summary_type.hapmap_file' hapmap.gz && + contamination -a output.icnts.gz -h hapmap.gz 2> 'contamination_report.out' && + python3 '$__tool_directory__/print_x_contamination.py' 'contamination_report.out' + + #elif $summary_type.summary_type_selector == 'error_est_poly': + -doCounts 1 -doError 1 -doMajorMinor 2 + @BAM_PARAM@ + @COUNTS_PARAM@ + -minSites $summary_type.minSites + -minPhat $summary_type.minPhat + -eps $summary_type.eps + #if $summary_type.errors: + -errors '$summary_type.errors' + #end if + + #elif $summary_type.summary_type_selector == 'error_est_out': + @BAM_PARAM@ + -doAncError $summary_type.doAncError + -ref '$summary_type.ref' + -anc '$summary_type.anc' + && Rscript '$__tool_directory__/estError.R' file=output.ancError > output.estError.txt + + #elif $summary_type.summary_type_selector == 'allele_counts': + -doCounts 1 + @BAM_PARAM@ + @COUNTS_PARAM@ + #if $summary_type.minQfile: + -minQfile '$summary_type.minQfile' + #end if + -dumpCounts $summary_type.dumpCounts + #if $summary_type.doDepth: + -doDepth 1 + -maxDepth $summary_type.maxDepth + #end if + #if $summary_type.doQsDist: + -doQsDist 1 + #end if + + #elif $summary_type.summary_type_selector == 'mismatch': + -doMisMatch 1 + @BAM_PARAM@ + -ref '$summary_type.ref' + #end if + ]]> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + summary_type['summary_type_selector'] == 'contamination' + + + summary_type['summary_type_selector'] == 'contamination' + + + summary_type['summary_type_selector'] == 'contamination' and summary_type['generate_json'] + + + + summary_type['summary_type_selector'] == 'error_est_poly' + + + summary_type['summary_type_selector'] == 'error_est_out' + + + summary_type['summary_type_selector'] == 'error_est_out' + + + + summary_type['summary_type_selector'] == 'allele_counts' and summary_type['dumpCounts'] != '0' + + + summary_type['summary_type_selector'] == 'allele_counts' and summary_type['dumpCounts'] in ['2','3','4'] + + + summary_type['summary_type_selector'] == 'allele_counts' and summary_type['doDepth'] + + + summary_type['summary_type_selector'] == 'allele_counts' and summary_type['doDepth'] + + + summary_type['summary_type_selector'] == 'allele_counts' and summary_type['doQsDist'] + + + + summary_type['summary_type_selector'] == 'mismatch' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tools/angsd/estError.R b/tools/angsd/estError.R new file mode 100644 index 00000000000..09f43253cf1 --- /dev/null +++ b/tools/angsd/estError.R @@ -0,0 +1,302 @@ +bases<-c("A","C","G","T") +b <- c(bases,"N") +########### do not change ################3 +l<-commandArgs(TRUE) +getArgs<-function(x,l) + unlist(strsplit(grep(paste("^",x,"=",sep=""),l,val=T),"="))[2] +Args<-function(l,args){ + if(! all(sapply(strsplit(l,"="),function(x)x[1])%in%names(args))){ + cat("Error -> ",l[!sapply(strsplit(l,"="),function(x)x[1])%in%names(args)]," is not a valid argument") + q("no") + } + arguments<-list() + for(a in names(args)) + arguments[[a]]<-getArgs(a,l) + + if(any(!names(args)%in%names(arguments)&sapply(args,is.null))){ + cat("Error -> ",names(args)[!names(args)%in%names(arguments)&sapply(args,is.null)]," is not optional!\n") + q("no") + } + for(a in names(args)) + if(is.null(arguments[[a]])) + arguments[[a]]<-args[[match(a,names(args))]] + + + arguments +} + +print.args<-function(args,des){ + if(missing(des)){ + des<-as.list(rep("",length(args))) + names(des)<-names(args) + } + cat("-> Needed arguments:\n") + mapply(function(x)cat("\t",x,":",des[[x]],"\n"),cbind(names(args)[sapply(args,is.null)])) + cat("-> Optional arguments (defaults):\n") + mapply(function(x)cat("\t",x," (",args[[x]],")",":",des[[x]],"\n"),cbind(names(args)[!sapply(args,is.null)])) + q("no") +} + +## choose your parameters and defaults +## NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments +args<-list(file=NULL, + out="errorEst", + indNames="ind", + nIter=100, + subset=NA, + main="Error rate using an outgroup and a high quality genome", + height=7, + width=12, + srt=90, + cex=1, + doPng=FALSE, + ylim1=NA, + ylim2=NA + ) +#if no argument aree given prints the need arguments and the optional ones with default +des<-list(file="the ancError File", + out="Name of the out files", + indNames="postFix, file with names or comma seperated names of individuals", + nIter="Numer of optimazation attemps", + subset="comma seperated numbers of the individuals to include", + width="width of the pdf", + height="height of the pdf", + srt="angle of ind names", + cex="scale of names", + doPng="Make png instead of pdf", + ylim1="Set max ylim for barplot1", + ylim2="Set max ylim for barplot2" + ) + +#######get arguments and add to workspace +### do not change +if(length(l)==0) print.args(args,des) +attach(Args(l,args)) +args <- commandArgs(TRUE) +if(length(args)==0){ + cat(" Arguments: output prefix\n") + q("no") +} + +################################### +tol21rainbow= c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#117744", "#44AA77", "#88CCAA", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455", "#DD7788") + + + +width=as.numeric(width) +height=as.numeric(height) + +if(!is.na(ylim1)) + ylim1=as.numeric(ylim1) +if(!is.na(ylim2)) + ylim2=as.numeric(ylim2) + + +palette(tol21rainbow[c(1:8,10:21)]) +srt<-as.numeric(srt) + + +nIter<-as.integer(nIter) +cex=as.numeric(cex) +cat("----------\nfile: ",file," out: ",out," nIter: ",nIter," subset: ",subset,"\n width: ",width," height: ",height," srt: ",srt," cex: ",cex," doPng: ",doPng,"\n-----------\n") + +r<-as.matrix(read.table(file)) + +if(!is.na(subset)){ + subset<- as.integer(unlist(strsplit(subset,","))) + print(subset) + r<-r[subset,] +} + +nInd<-nrow(r) +cat("Number of individuals read:",nInd,"\n") +{ + if(length(grep(",",indNames))>0) + indNames<-unlist(strsplit(indNames,",")) + else { + options("warn"=-1) + try(indNames<-basename(scan(indNames,what="theFuck")),silent=TRUE) + options("warn"=0) + } +} + +if(length(indNames)==1&nInd>1){ + indNames<-paste(indNames,1:nInd,sep="") +} +cat("Ind names:\n") +print(indNames) + +if(length(indNames)!=nInd){ + cat("Error: Wrong number of ind Names\n") + q("no") +} + +getMat<-function(x){ + m<-array(0,dim=c(5,5,5),dimnames=list(b,b,b)) + for(s in 0:4) + for(p in 0:4) + for(a in 0:4) + m[a+1,p+1,s+1]<-x[a*25+p*5+s+1] + m +} + +logLike<-function(x,Xch,Pch){ + eMat<-matrix(0,4,4) + eMat[-c(1,6,11,16)]<-x + diag(eMat)<-1-rowSums(eMat) + P <- Pch %*% eMat + ll <- -sum(log(P)*Xch) + return(ll) +} + +res<-NULL +for(j in 1:nInd){ + m<-getMat(r[j,]) + print(m) + + ##remove if missing + m<-m[-5,-5,-5] + + Pch<-matrix(0,4,4) + for(i in 1:4) + Pch<-Pch+m[,,i] + print(Pch) + Pch<-Pch/rowSums(Pch) + cat("Matrix of allele change probabilities.\nA diagonal close to 1 denotes a low probability of errors.\n\n") + print(Pch) + + Xch<-matrix(0,4,4) + for(i in 1:4) + Xch<-Xch+m[,i,] + cat("\n\nObserved allele changes:\n\n") + print(Xch) + maxErr=0.02; + + conv<-list(par=1,objective=Inf) + + for(i in 1:nIter){ + Tempconv <- nlminb(runif(12)/100,logLike,upper=rep(maxErr,12),lower=rep(1e-10,12),Xch=Xch,Pch=Pch) + if( any(Tempconv$par > maxErr*0.99) ) + maxErr<-min(1,maxErr*10); + if(Tempconv$objective",rep(bases,each=4))[-c(1,6,11,16)] + +{ + if(doPng==FALSE){ + pdf(paste(out,".pdf",sep=""),width=width,height=height) + } else + bitmap(paste(out,".png",sep=""),w=width,h=height,res=300) +} + +{ +if(is.na(ylim1)) + barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate") +else + barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate",ylim=c(0,ylim1)) +} +legend("top",paste(indNames,round(over*100,2),"%"),fill=1:nInd,bty="n") +dev.off() + +colnames(res)<-nam +rownames(res)<-indNames + +write.table(res,file=paste(out,".txt",sep=""),sep="\t") +{ +if(doPng==FALSE) + pdf(paste(out,"Overall.pdf",sep=""),w=width,h=height) +else + bitmap(paste(out,"Overall.png",sep=""),w=width,h=height,res=300) +} + +par(mar=par()$mar+c(3,0,0,0)) +{ +if(is.na(ylim2)){ + h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate") +}else{ + h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate",ylim=c(0,ylim2)) +} +} +text(h,rep(-max(over)/50,length(h)),indNames,xpd=T,srt=srt,adj=1,cex=cex) +dev.off() +{ + if(doPng==FALSE){ + cat("figure:",paste(out,".pdf",sep=""),"\n") + cat("figure2:",paste(out,"Overall.pdf",sep=""),"\n") + }else{ + cat("figure:",paste(out,".png",sep=""),"\n") + cat("figure2:",paste(out,"Overall.png",sep=""),"\n") + } + +} +cat("table:",paste(out,".txt",sep=""),"\n") + +write.table(cbind(paste(indNames,round(over*100,4),"%")),file=paste(out,".txt",sep=""),sep="\t",append=T,col=F,row=F) + +chrFile<-paste(file,"Chr",sep="") + +if(file.info(chrFile)$size!=0 & !is.na(file.info(chrFile)$size)){ + + overList<-list() + con<-file(chrFile,"r") + while(length(cchr<-scan(con,nlines=1,what="theFuck"))>0){ + r<-read.table(con,nrow=nInd) + r<-as.matrix(r) + overList[[cchr[2]]]<-getover(r,nInd) + } + close(con) + + rbind.list<-function (l) + { + n <- length(l) + results <- NULL + for (i in 1:n) { + results <- rbind(results, l[[i]]) + } + results + } + res<-rbind.list(overList) + chrNames<-names(overList) + rownames(res) <- chrNames + colnames(res) <- indNames + pdf(paste(out,"OverChr.pdf",sep=""),w=width,h=height) + for(i in 1:nInd){ + dotchart(res[,i],chrNames,xlab="Error rates",col=3,main=indNames[i]) + dotchart(res[res[,i]>0,i],chrNames,xlab="Error rates (positive)",col=3,main=indNames[i]) + } + dev.off() + write.table(res,file=paste(out,"OverChr.txt",sep=""),qu=F,row=T,col=T) +} + +#warnings() diff --git a/tools/angsd/macros.xml b/tools/angsd/macros.xml index 848054e473c..6fd26f9c9f8 100644 --- a/tools/angsd/macros.xml +++ b/tools/angsd/macros.xml @@ -9,6 +9,7 @@ angsd samtools python + r-base @@ -19,4 +20,46 @@ 10.7717/peerj.10947 + + + + + + + ^([\w\d\._:-]+)?$ + + + + + + + + + + + + + -minMapQ ${summary_type.minMapQ} + -minQ ${summary_type.minQ} + -trim ${summary_type.trim} + #if str($summary_type.r): + -r '$summary_type.r' + #end if + + + + -minInd ${summary_type.minInd} + #if str($summary_type.setMinDepth) != '-1': + -setMinDepth ${summary_type.setMinDepth} + #end if + #if str($summary_type.setMaxDepth) != '-1': + -setMaxDepth ${summary_type.setMaxDepth} + #end if + #if str($summary_type.setMinDepthInd) != '-1': + -setMinDepthInd ${summary_type.setMinDepthInd} + #end if + #if str($summary_type.setMaxDepthInd) != '-1': + -setMaxDepthInd ${summary_type.setMaxDepthInd} + #end if + \ No newline at end of file diff --git a/tools/angsd/test-data/output.icnts.gz b/tools/angsd/test-data/output.icnts.gz deleted file mode 100644 index 7bff3041ba4334074824f2dd2ccb7da2397aad73..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 309 zcmb2|=3rp}f&Xj_PR>jW5{$W*T>B0?2)I3*xjJ=jO#Y0Sif7LFl%F({wl12o-hA4; z_|U(xuO|Ec(Xcsw@Vm(&r@C#O@Bg(wm)R=L*u87|i}Kc6*&C!cM5pXLe5ZN?`xU!G z)lJ`fbo3L}?|3i0r@TV`gWbgU+&}j{{M+z%(Kr5k_7Or`H7hbJBr75}$b6f0JdP{D z^^5kxc*&ny387c68EofU&v=a`dttnE9Q)5*AF6IR%HLzi;V@&8R{qR#Ph@8Lin)iQ z8*~qb?RX#c-ej7^jM;sveD^aSyxkjIVfp3sVeJFjZV@NEXGFX0=ZX+r!xqM#-C%k! xDj|0RUxfG?_G`!2u_ZWP={s=j#!H}aIou38+MQ(B85rcz0!EsF85~X^0sz(EgX#bP diff --git a/tools/angsd/test-data/test.bam b/tools/angsd/test-data/test.bam new file mode 100644 index 0000000000000000000000000000000000000000..2ad644a797dcc51f1bf17d2ecac56c2201df1ef3 GIT binary patch literal 2989 zcmV;e3sUqSiwFb&00000{{{d;LjnNt1g%xga-&8NR(2}~`H*Lr+n$gZ1VX5il(87g z#QYhMEvH;HHNePL76}Egoy{BME%F+9lbrGzxuw#Bgec8=PP?`>HZ!)Vc^UFm@x*dUOR?aTZOyde zKR@zy$}3u|X~t?+@tUfX*5%_mg#%}_8@X=~Byh&IL>j{ zf7EX?C4*6>q)H~~SjUpO&}7LN@kG)_9AlmFg^{K`0JoI^-Rq*rD^Z@8bZII-ngsqTNd@)}PjI%vooMdchsJfKu8JARM=#r5tiKHr{ zp*&I+YL=>pd|g(4v>rG4Qq-an6t8(%EsOPd(8Ilk>oEjLQfyAAwl}4L@7R-DWac>- zZ7-O`$YmS)ysRJS*^PNvaD}pUMzKukifuMf^pdlzqRTZsGhMnU$`z|AMDb^_2=k>{ z4Esty-D@#aQdO7LM3Pf2lT;OBiDil=VJ#i<5j2ISHnoSmE{Q{X647(pAzz#wGWYd; z_B7KO(+5}@B#9xZ1~(*;Q_@F#fVBb7k_F#o4~gqv6Kc%@5Y&uGgy8Go;4q_T<`NqD zVN9&nxQEeb;)kLTQ5?c{dW!%}{d3SXBJ^rzI(BSBFpxEf;SJ{{0}bpsfxR8YrZ*u7 z%QC>ONqFW*c6^IS82TYtu4Xog+r-xp&@yA&_e4fWOgt!m1d@5N(9oRIh+NIUN4oW8&JKSJ9c__>MdRH8MZg2OThPQICLB9Hl$y*bTy}fNx#QnZ*SE;AWvlir z6&2W*qeaZ196*cO9-VoCX@|sWxJE~ay=&swK|mTMUI49A#F>G@VQktSR7!475M9~< zNE${l^{B8!FTN{kNA%M3Z@iQlIX&BEQE!&re)D!N`=I+zdQF% zpCI(7A$-6gfldT|7J7AAw2*1&&jRY1E_sb<>MmRmXiDbj+hS8MSIb{{hEkX!k0nRN z0_Ci#QTOfZH=y_ZLgBn{=oCie)F`w zVfXv$p?JPnW*J`}^}`R=CAG-80|CQV4%<6)h-r1K>gV$euJX?&+J3zqYlEY)_O#^B zNP0k0g0{CGHhN-C&Bt5*vZUeWdGd9q(|Hen1$XZ6|9aOfmCnjl`SAgg?L-^kVQM_#wjOSMHXWW~7@cgNnx0)aA6}lnaBgmPX|RZ9+Nb=p zv-3-Xg_*%Sm)@Jb(C;tsUk3fz;p6_|;N8oEi*tj;{-ulkv$IRHi-V;;>-U!yE??yD zINL|P_8fY3dJgsYJ^{Py&-{vpEIhcuS@*({i)-godW+NmSLf zKu(iT*%#-{C2+ocqZk0IoTX`m#EF8Q%6 zYepc)T402Mq&7COauqxeD!9W{&@7XXzk~|*-r*{^0)!x>7xaQjJc^ry0yiczk;+w| za~CLyRt>bts-bmM74WT#AE565=n+6CxQ%hk++zHrJ-$~7;@aID|2R5a0Kd%zK|D-d z-1IOewvFz=dnT0*tZHY~+O!xZfRXQk5$u~?M$07Rups2t??A|NLqc9I8XO&un$&^d z;5sHl_rW3^#Rz>}23i7H(`DHR4Jk>JC{ELkBFW>8r(QW#BRxEaI90fdnU7i5>Ef+s zY+H}-A!Bdc7Z%!yra?j-8)GMXPq0cmDz}i@c+5o0H?>B+XSQya z9!LKWm=EmAF^{9kVsR8hz#vTPf-FABDp|&C%Os zqSu8-?o2#F|5C9N=LG~dDRn*CbQ}}i29H#eD2Vl{6zh?e#MRIT$7mh+5mgf@{*U%} z1{m7G?a`ubiXBu+tkfnVEb1in^V%x;CuuIWB;9FCYsvhfs0V!5blnrr1 zPFdVN!lhxBDUCk`=r010M*Zrs{0$_dMMJqF7`XW%`< zdNN8ULnTO*K=J+RD!}TR0`U0RyJepJga`{KXL9feIYkoz1qSRUfVY_Kpn0A;*Iq|A zF01~{OM}zW`UYqLag|65612Rwf|eII5eJ&K!yro%#?Yr+6Hi*j`GCITTi$oDQuF0* z!NCl0uwVFs74(Ja0(y8#n-~K4(RGzk=jU8s4h;1br~f2HmE22h0#}H`Ai0#D@bV+H4R9_Buxi0j15P5bB$1XyEpp;T8Q$NYIW)QX1(|6d0it=z5SQO8OK9#+pm?X{hJJglFmy zV=sD4Tk^-2ng zlCKq36!#D@p-I6G0Vg?hs3=A9a66~#N8 zF;G-*>Y(kw{Ly7N6BIiF8aqgA6Jtgl@d*Bj-(eWCSmF_von^L3o0zSw4q&`R8q3E; z`RJs7i`#if{{Es}0_cvq@EF3mV`VLppBlh=tUzzy73*U~*GW5f7z@!qv;7EqSkv6Gq2t02PwbH@w6~6oFI}ZMOcJy|^sTp?Q zLLHz0zISrb8D8iqhNfsR6)`$uDw5HVYGL?vW&={Ti@^bDGX}F~+QlV~QfK}F+%i5Q jni~KBABzYC000000RIL6LPG)o8vp|U0000000000XW*&S literal 0 HcmV?d00001 diff --git a/tools/angsd/test-data/test.bam.bai b/tools/angsd/test-data/test.bam.bai new file mode 100644 index 0000000000000000000000000000000000000000..c583bb5db0799eda236f1e8668e1c7eb3da2db60 GIT binary patch literal 2760 zcmZ>A^kigZzyusXp$1_F21X!VDfNj3qGBR9i1wPqz`z6+KoysQD1wt!Q2h)H5Or*f z5V{&F4&(DML&R&K;xN8H8$`SoDh}hva6!cDpyDt-qc}vo9x4vw!`yEH@(k21qckL} sM#Et=93UBXG#p05VYGZ4Ee}V_fzfh+e&qmq8xXyXhu(&R#0@4105C2^0ssI2 literal 0 HcmV?d00001 diff --git a/tools/angsd/test-data/test.phe b/tools/angsd/test-data/test.phe deleted file mode 100644 index d00491fd7e5..00000000000 --- a/tools/angsd/test-data/test.phe +++ /dev/null @@ -1 +0,0 @@ -1 diff --git a/tools/angsd/test-data/test.sam b/tools/angsd/test-data/test.sam index e9169938c47..b603f0e4921 100644 --- a/tools/angsd/test-data/test.sam +++ b/tools/angsd/test-data/test.sam @@ -1,8 +1,8 @@ @HD VN:1.0 SO:coordinate @SQ SN:1 LN:30 @RG ID:RG0 LB:lb1 SM:sm0 PI:450 CN:SC PL:ILLUMINA -rd1sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN ]asdf]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 -rd2sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN 1234]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 -rd3sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN QWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 -rd4sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN tWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 -rd5sm0 0 1 1 60 30M * 0 0 ATGTNACGTNACGTNACGTNACGTNNNNNN DWER]]]]]]]]]]]]]]]]]]]]]]]]]] RG:Z:RG0 \ No newline at end of file +rd1sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 +rd2sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 +rd3sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 +rd4sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 +rd5sm0 0 1 1 60 30M * 0 0 ATGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 \ No newline at end of file diff --git a/tools/angsd/test-data/anc.fa b/tools/angsd/test-data/test_anc.fa similarity index 100% rename from tools/angsd/test-data/anc.fa rename to tools/angsd/test-data/test_anc.fa diff --git a/tools/angsd/test-data/ref.fa b/tools/angsd/test-data/test_ref.fa similarity index 100% rename from tools/angsd/test-data/ref.fa rename to tools/angsd/test-data/test_ref.fa From e3ae69da421e5956bf242f65e6c239ee04bab40f Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 9 Aug 2026 11:03:22 +0000 Subject: [PATCH 4/6] fix: resolve bam inputs --- tools/angsd/{angsd_summary.xml => angsd.xml} | 12 +- tools/angsd/estError.R | 618 +++++++++--------- tools/angsd/test-data/test.bam.bai | Bin 2760 -> 0 bytes tools/angsd/test-data/test.sam | 8 - tools/angsd/test-data/test_chr1.bam | Bin 0 -> 323 bytes .../test-data/{test.bam => test_chrX.bam} | Bin 6 files changed, 322 insertions(+), 316 deletions(-) rename tools/angsd/{angsd_summary.xml => angsd.xml} (97%) delete mode 100644 tools/angsd/test-data/test.bam.bai delete mode 100644 tools/angsd/test-data/test.sam create mode 100644 tools/angsd/test-data/test_chr1.bam rename tools/angsd/test-data/{test.bam => test_chrX.bam} (100%) diff --git a/tools/angsd/angsd_summary.xml b/tools/angsd/angsd.xml similarity index 97% rename from tools/angsd/angsd_summary.xml rename to tools/angsd/angsd.xml index 64ee6334f34..c80e88359c3 100644 --- a/tools/angsd/angsd_summary.xml +++ b/tools/angsd/angsd.xml @@ -1,4 +1,4 @@ - + Calculate statistics from BAM files @@ -181,7 +181,7 @@ - + @@ -207,7 +207,7 @@ - + @@ -224,7 +224,7 @@ - + @@ -246,7 +246,7 @@ - + @@ -283,7 +283,7 @@ - + diff --git a/tools/angsd/estError.R b/tools/angsd/estError.R index 09f43253cf1..57430a24d22 100644 --- a/tools/angsd/estError.R +++ b/tools/angsd/estError.R @@ -1,302 +1,316 @@ -bases<-c("A","C","G","T") -b <- c(bases,"N") -########### do not change ################3 -l<-commandArgs(TRUE) -getArgs<-function(x,l) - unlist(strsplit(grep(paste("^",x,"=",sep=""),l,val=T),"="))[2] -Args<-function(l,args){ - if(! all(sapply(strsplit(l,"="),function(x)x[1])%in%names(args))){ - cat("Error -> ",l[!sapply(strsplit(l,"="),function(x)x[1])%in%names(args)]," is not a valid argument") - q("no") - } - arguments<-list() - for(a in names(args)) - arguments[[a]]<-getArgs(a,l) - - if(any(!names(args)%in%names(arguments)&sapply(args,is.null))){ - cat("Error -> ",names(args)[!names(args)%in%names(arguments)&sapply(args,is.null)]," is not optional!\n") - q("no") - } - for(a in names(args)) - if(is.null(arguments[[a]])) - arguments[[a]]<-args[[match(a,names(args))]] - - - arguments -} - -print.args<-function(args,des){ - if(missing(des)){ - des<-as.list(rep("",length(args))) - names(des)<-names(args) - } - cat("-> Needed arguments:\n") - mapply(function(x)cat("\t",x,":",des[[x]],"\n"),cbind(names(args)[sapply(args,is.null)])) - cat("-> Optional arguments (defaults):\n") - mapply(function(x)cat("\t",x," (",args[[x]],")",":",des[[x]],"\n"),cbind(names(args)[!sapply(args,is.null)])) - q("no") -} - -## choose your parameters and defaults -## NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments -args<-list(file=NULL, - out="errorEst", - indNames="ind", - nIter=100, - subset=NA, - main="Error rate using an outgroup and a high quality genome", - height=7, - width=12, - srt=90, - cex=1, - doPng=FALSE, - ylim1=NA, - ylim2=NA - ) -#if no argument aree given prints the need arguments and the optional ones with default -des<-list(file="the ancError File", - out="Name of the out files", - indNames="postFix, file with names or comma seperated names of individuals", - nIter="Numer of optimazation attemps", - subset="comma seperated numbers of the individuals to include", - width="width of the pdf", - height="height of the pdf", - srt="angle of ind names", - cex="scale of names", - doPng="Make png instead of pdf", - ylim1="Set max ylim for barplot1", - ylim2="Set max ylim for barplot2" - ) - -#######get arguments and add to workspace -### do not change -if(length(l)==0) print.args(args,des) -attach(Args(l,args)) -args <- commandArgs(TRUE) -if(length(args)==0){ - cat(" Arguments: output prefix\n") - q("no") -} - -################################### -tol21rainbow= c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#117744", "#44AA77", "#88CCAA", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455", "#DD7788") - - - -width=as.numeric(width) -height=as.numeric(height) - -if(!is.na(ylim1)) - ylim1=as.numeric(ylim1) -if(!is.na(ylim2)) - ylim2=as.numeric(ylim2) - - -palette(tol21rainbow[c(1:8,10:21)]) -srt<-as.numeric(srt) - - -nIter<-as.integer(nIter) -cex=as.numeric(cex) -cat("----------\nfile: ",file," out: ",out," nIter: ",nIter," subset: ",subset,"\n width: ",width," height: ",height," srt: ",srt," cex: ",cex," doPng: ",doPng,"\n-----------\n") - -r<-as.matrix(read.table(file)) - -if(!is.na(subset)){ - subset<- as.integer(unlist(strsplit(subset,","))) - print(subset) - r<-r[subset,] -} - -nInd<-nrow(r) -cat("Number of individuals read:",nInd,"\n") -{ - if(length(grep(",",indNames))>0) - indNames<-unlist(strsplit(indNames,",")) - else { - options("warn"=-1) - try(indNames<-basename(scan(indNames,what="theFuck")),silent=TRUE) - options("warn"=0) - } -} - -if(length(indNames)==1&nInd>1){ - indNames<-paste(indNames,1:nInd,sep="") -} -cat("Ind names:\n") -print(indNames) - -if(length(indNames)!=nInd){ - cat("Error: Wrong number of ind Names\n") - q("no") -} - -getMat<-function(x){ - m<-array(0,dim=c(5,5,5),dimnames=list(b,b,b)) - for(s in 0:4) - for(p in 0:4) - for(a in 0:4) - m[a+1,p+1,s+1]<-x[a*25+p*5+s+1] - m -} - -logLike<-function(x,Xch,Pch){ - eMat<-matrix(0,4,4) - eMat[-c(1,6,11,16)]<-x - diag(eMat)<-1-rowSums(eMat) - P <- Pch %*% eMat - ll <- -sum(log(P)*Xch) - return(ll) -} - -res<-NULL -for(j in 1:nInd){ - m<-getMat(r[j,]) - print(m) - - ##remove if missing - m<-m[-5,-5,-5] - - Pch<-matrix(0,4,4) - for(i in 1:4) - Pch<-Pch+m[,,i] - print(Pch) - Pch<-Pch/rowSums(Pch) - cat("Matrix of allele change probabilities.\nA diagonal close to 1 denotes a low probability of errors.\n\n") - print(Pch) - - Xch<-matrix(0,4,4) - for(i in 1:4) - Xch<-Xch+m[,i,] - cat("\n\nObserved allele changes:\n\n") - print(Xch) - maxErr=0.02; - - conv<-list(par=1,objective=Inf) - - for(i in 1:nIter){ - Tempconv <- nlminb(runif(12)/100,logLike,upper=rep(maxErr,12),lower=rep(1e-10,12),Xch=Xch,Pch=Pch) - if( any(Tempconv$par > maxErr*0.99) ) - maxErr<-min(1,maxErr*10); - if(Tempconv$objective",rep(bases,each=4))[-c(1,6,11,16)] - -{ - if(doPng==FALSE){ - pdf(paste(out,".pdf",sep=""),width=width,height=height) - } else - bitmap(paste(out,".png",sep=""),w=width,h=height,res=300) -} - -{ -if(is.na(ylim1)) - barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate") -else - barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate",ylim=c(0,ylim1)) -} -legend("top",paste(indNames,round(over*100,2),"%"),fill=1:nInd,bty="n") -dev.off() - -colnames(res)<-nam -rownames(res)<-indNames - -write.table(res,file=paste(out,".txt",sep=""),sep="\t") -{ -if(doPng==FALSE) - pdf(paste(out,"Overall.pdf",sep=""),w=width,h=height) -else - bitmap(paste(out,"Overall.png",sep=""),w=width,h=height,res=300) -} - -par(mar=par()$mar+c(3,0,0,0)) -{ -if(is.na(ylim2)){ - h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate") -}else{ - h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate",ylim=c(0,ylim2)) -} -} -text(h,rep(-max(over)/50,length(h)),indNames,xpd=T,srt=srt,adj=1,cex=cex) -dev.off() -{ - if(doPng==FALSE){ - cat("figure:",paste(out,".pdf",sep=""),"\n") - cat("figure2:",paste(out,"Overall.pdf",sep=""),"\n") - }else{ - cat("figure:",paste(out,".png",sep=""),"\n") - cat("figure2:",paste(out,"Overall.png",sep=""),"\n") - } - -} -cat("table:",paste(out,".txt",sep=""),"\n") - -write.table(cbind(paste(indNames,round(over*100,4),"%")),file=paste(out,".txt",sep=""),sep="\t",append=T,col=F,row=F) - -chrFile<-paste(file,"Chr",sep="") - -if(file.info(chrFile)$size!=0 & !is.na(file.info(chrFile)$size)){ - - overList<-list() - con<-file(chrFile,"r") - while(length(cchr<-scan(con,nlines=1,what="theFuck"))>0){ - r<-read.table(con,nrow=nInd) - r<-as.matrix(r) - overList[[cchr[2]]]<-getover(r,nInd) - } - close(con) - - rbind.list<-function (l) - { - n <- length(l) - results <- NULL - for (i in 1:n) { - results <- rbind(results, l[[i]]) - } - results - } - res<-rbind.list(overList) - chrNames<-names(overList) - rownames(res) <- chrNames - colnames(res) <- indNames - pdf(paste(out,"OverChr.pdf",sep=""),w=width,h=height) - for(i in 1:nInd){ - dotchart(res[,i],chrNames,xlab="Error rates",col=3,main=indNames[i]) - dotchart(res[res[,i]>0,i],chrNames,xlab="Error rates (positive)",col=3,main=indNames[i]) - } - dev.off() - write.table(res,file=paste(out,"OverChr.txt",sep=""),qu=F,row=T,col=T) -} - -#warnings() +bases <- c("A", "C", "G", "T") +b <- c(bases, "N") +########### do not change ################3 +l <- commandArgs(TRUE) +getArgs <- function(x, l) { + unlist(strsplit(grep(paste("^", x, "=", sep = ""), l, val = T), "="))[2] +} +Args <- function(l, args) { + if (!all(sapply(strsplit(l, "="), function(x) x[1]) %in% names(args))) { + cat("Error -> ", l[!sapply(strsplit(l, "="), function(x) x[1]) %in% names(args)], " is not a valid argument") + q("no") + } + arguments <- list() + for (a in names(args)) { + arguments[[a]] <- getArgs(a, l) + } + + if (any(!names(args) %in% names(arguments) & sapply(args, is.null))) { + cat("Error -> ", names(args)[!names(args) %in% names(arguments) & sapply(args, is.null)], " is not optional!\n") + q("no") + } + for (a in names(args)) { + if (is.null(arguments[[a]])) { + arguments[[a]] <- args[[match(a, names(args))]] + } + } + + + arguments +} + +print.args <- function(args, des) { + if (missing(des)) { + des <- as.list(rep("", length(args))) + names(des) <- names(args) + } + cat("-> Needed arguments:\n") + mapply(function(x) cat("\t", x, ":", des[[x]], "\n"), cbind(names(args)[sapply(args, is.null)])) + cat("-> Optional arguments (defaults):\n") + mapply(function(x) cat("\t", x, " (", args[[x]], ")", ":", des[[x]], "\n"), cbind(names(args)[!sapply(args, is.null)])) + q("no") +} + +## choose your parameters and defaults +## NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments +args <- list( + file = NULL, + out = "errorEst", + indNames = "ind", + nIter = 100, + subset = NA, + main = "Error rate using an outgroup and a high quality genome", + height = 7, + width = 12, + srt = 90, + cex = 1, + doPng = FALSE, + ylim1 = NA, + ylim2 = NA +) +# if no argument aree given prints the need arguments and the optional ones with default +des <- list( + file = "the ancError File", + out = "Name of the out files", + indNames = "postFix, file with names or comma seperated names of individuals", + nIter = "Numer of optimazation attemps", + subset = "comma seperated numbers of the individuals to include", + width = "width of the pdf", + height = "height of the pdf", + srt = "angle of ind names", + cex = "scale of names", + doPng = "Make png instead of pdf", + ylim1 = "Set max ylim for barplot1", + ylim2 = "Set max ylim for barplot2" +) + +####### get arguments and add to workspace +### do not change +if (length(l) == 0) print.args(args, des) +attach(Args(l, args)) +args <- commandArgs(TRUE) +if (length(args) == 0) { + cat(" Arguments: output prefix\n") + q("no") +} + +################################### +tol21rainbow <- c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#117744", "#44AA77", "#88CCAA", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455", "#DD7788") + + +width <- as.numeric(width) +height <- as.numeric(height) + +if (!is.na(ylim1)) { + ylim1 <- as.numeric(ylim1) +} +if (!is.na(ylim2)) { + ylim2 <- as.numeric(ylim2) +} + + +palette(tol21rainbow[c(1:8, 10:21)]) +srt <- as.numeric(srt) + + +nIter <- as.integer(nIter) +cex <- as.numeric(cex) +cat("----------\nfile: ", file, " out: ", out, " nIter: ", nIter, " subset: ", subset, "\n width: ", width, " height: ", height, " srt: ", srt, " cex: ", cex, " doPng: ", doPng, "\n-----------\n") + +r <- as.matrix(read.table(file)) + +if (!is.na(subset)) { + subset <- as.integer(unlist(strsplit(subset, ","))) + print(subset) + r <- r[subset, ] +} + +nInd <- nrow(r) +cat("Number of individuals read:", nInd, "\n") +{ + if (length(grep(",", indNames)) > 0) { + indNames <- unlist(strsplit(indNames, ",")) + } else { + options("warn" = -1) + try(indNames <- basename(scan(indNames, what = "theFuck")), silent = TRUE) + options("warn" = 0) + } +} + +if (length(indNames) == 1 & nInd > 1) { + indNames <- paste(indNames, 1:nInd, sep = "") +} +cat("Ind names:\n") +print(indNames) + +if (length(indNames) != nInd) { + cat("Error: Wrong number of ind Names\n") + q("no") +} + +getMat <- function(x) { + m <- array(0, dim = c(5, 5, 5), dimnames = list(b, b, b)) + for (s in 0:4) { + for (p in 0:4) { + for (a in 0:4) { + m[a + 1, p + 1, s + 1] <- x[a * 25 + p * 5 + s + 1] + } + } + } + m +} + +logLike <- function(x, Xch, Pch) { + eMat <- matrix(0, 4, 4) + eMat[-c(1, 6, 11, 16)] <- x + diag(eMat) <- 1 - rowSums(eMat) + P <- Pch %*% eMat + ll <- -sum(log(P) * Xch) + return(ll) +} + +res <- NULL +for (j in 1:nInd) { + m <- getMat(r[j, ]) + print(m) + + ## remove if missing + m <- m[-5, -5, -5] + + Pch <- matrix(0, 4, 4) + for (i in 1:4) { + Pch <- Pch + m[, , i] + } + print(Pch) + Pch <- Pch / rowSums(Pch) + cat("Matrix of allele change probabilities.\nA diagonal close to 1 denotes a low probability of errors.\n\n") + print(Pch) + + Xch <- matrix(0, 4, 4) + for (i in 1:4) { + Xch <- Xch + m[, i, ] + } + cat("\n\nObserved allele changes:\n\n") + print(Xch) + maxErr <- 0.02 + + conv <- list(par = 1, objective = Inf) + + for (i in 1:nIter) { + Tempconv <- nlminb(runif(12) / 100, logLike, upper = rep(maxErr, 12), lower = rep(1e-10, 12), Xch = Xch, Pch = Pch) + if (any(Tempconv$par > maxErr * 0.99)) { + maxErr <- min(1, maxErr * 10) + } + if (Tempconv$objective < conv$objective) { + conv <- Tempconv + } + } + + res <- rbind(res, conv$par) +} + +getover <- function(r, nInd) { + over <- NULL + for (j in 1:nInd) { + m <- getMat(r[j, ]) + + ## remove if missing + m <- m[-5, -5, -5] + Xch <- matrix(0, 4, 4) + for (i in 1:4) { + Xch <- Xch + m[, i, ] + } + + Pch <- matrix(0, 4, 4) + for (i in 1:4) { + Pch <- Pch + m[, , i] + } + + N1dot <- sum(Xch) - sum(diag(Xch)) + Ndot1 <- sum(Pch) - sum(diag(Pch)) + Ndot0 <- sum(diag(Pch)) + err <- (N1dot - Ndot1) / (Ndot0 - Ndot1) + over <- c(over, err) + } + over +} +over <- getover(r, nInd) +nam <- paste(rep(bases, 4), "->", rep(bases, each = 4))[-c(1, 6, 11, 16)] + +{ + if (doPng == FALSE) { + pdf(paste(out, ".pdf", sep = ""), width = width, height = height) + } else { + bitmap(paste(out, ".png", sep = ""), w = width, h = height, res = 300) + } +} + +{ + if (is.na(ylim1)) { + barplot(res, beside = T, col = 1:nInd, names = nam, main = main, ylab = "error rate") + } else { + barplot(res, beside = T, col = 1:nInd, names = nam, main = main, ylab = "error rate", ylim = c(0, ylim1)) + } +} +legend("top", paste(indNames, round(over * 100, 2), "%"), fill = 1:nInd, bty = "n") +dev.off() + +colnames(res) <- nam +rownames(res) <- indNames + +write.table(res, file = paste(out, ".txt", sep = ""), sep = "\t") +{ + if (doPng == FALSE) { + pdf(paste(out, "Overall.pdf", sep = ""), w = width, h = height) + } else { + bitmap(paste(out, "Overall.png", sep = ""), w = width, h = height, res = 300) + } +} + +par(mar = par()$mar + c(3, 0, 0, 0)) +{ + if (is.na(ylim2)) { + h <- barplot(over, col = 1:nInd, names = NULL, main = main, ylab = "error rate") + } else { + h <- barplot(over, col = 1:nInd, names = NULL, main = main, ylab = "error rate", ylim = c(0, ylim2)) + } +} +text(h, rep(-max(over) / 50, length(h)), indNames, xpd = T, srt = srt, adj = 1, cex = cex) +dev.off() +{ + if (doPng == FALSE) { + cat("figure:", paste(out, ".pdf", sep = ""), "\n") + cat("figure2:", paste(out, "Overall.pdf", sep = ""), "\n") + } else { + cat("figure:", paste(out, ".png", sep = ""), "\n") + cat("figure2:", paste(out, "Overall.png", sep = ""), "\n") + } +} +cat("table:", paste(out, ".txt", sep = ""), "\n") + +write.table(cbind(paste(indNames, round(over * 100, 4), "%")), file = paste(out, ".txt", sep = ""), sep = "\t", append = T, col = F, row = F) + +chrFile <- paste(file, "Chr", sep = "") + +if (file.info(chrFile)$size != 0 & !is.na(file.info(chrFile)$size)) { + overList <- list() + con <- file(chrFile, "r") + while (length(cchr <- scan(con, nlines = 1, what = "theFuck")) > 0) { + r <- read.table(con, nrow = nInd) + r <- as.matrix(r) + overList[[cchr[2]]] <- getover(r, nInd) + } + close(con) + + rbind.list <- function(l) { + n <- length(l) + results <- NULL + for (i in 1:n) { + results <- rbind(results, l[[i]]) + } + results + } + res <- rbind.list(overList) + chrNames <- names(overList) + rownames(res) <- chrNames + colnames(res) <- indNames + pdf(paste(out, "OverChr.pdf", sep = ""), w = width, h = height) + for (i in 1:nInd) { + dotchart(res[, i], chrNames, xlab = "Error rates", col = 3, main = indNames[i]) + dotchart(res[res[, i] > 0, i], chrNames, xlab = "Error rates (positive)", col = 3, main = indNames[i]) + } + dev.off() + write.table(res, file = paste(out, "OverChr.txt", sep = ""), qu = F, row = T, col = T) +} + +# warnings() diff --git a/tools/angsd/test-data/test.bam.bai b/tools/angsd/test-data/test.bam.bai deleted file mode 100644 index c583bb5db0799eda236f1e8668e1c7eb3da2db60..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2760 zcmZ>A^kigZzyusXp$1_F21X!VDfNj3qGBR9i1wPqz`z6+KoysQD1wt!Q2h)H5Or*f z5V{&F4&(DML&R&K;xN8H8$`SoDh}hva6!cDpyDt-qc}vo9x4vw!`yEH@(k21qckL} sM#Et=93UBXG#p05VYGZ4Ee}V_fzfh+e&qmq8xXyXhu(&R#0@4105C2^0ssI2 diff --git a/tools/angsd/test-data/test.sam b/tools/angsd/test-data/test.sam deleted file mode 100644 index b603f0e4921..00000000000 --- a/tools/angsd/test-data/test.sam +++ /dev/null @@ -1,8 +0,0 @@ -@HD VN:1.0 SO:coordinate -@SQ SN:1 LN:30 -@RG ID:RG0 LB:lb1 SM:sm0 PI:450 CN:SC PL:ILLUMINA -rd1sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 -rd2sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 -rd3sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 -rd4sm0 0 1 1 60 30M * 0 0 ACGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 -rd5sm0 0 1 1 60 30M * 0 0 ATGTNACGTNACGTNACGTNACGTNNNNNN IIIIIIIIIIIIIIIIIIIIIIIIIIIIII RG:Z:RG0 \ No newline at end of file diff --git a/tools/angsd/test-data/test_chr1.bam b/tools/angsd/test-data/test_chr1.bam new file mode 100644 index 0000000000000000000000000000000000000000..11c3aaf011656420ef93c5814d6af56527841027 GIT binary patch literal 323 zcmb2|=3rp}f&Xj_PR>jWCm34$efgRUL|RtwKCZ|k&JbNMwokO4bJCQk%o6@XO)NQ1 z(tEeG3q)p=Jbt7)*?eZcQQ(a8axD`Q<{yYu_^i2tIc|&Gt8Y=aj&WWQypq8nzGInR zz>*275lciWT4IknJk6-$GCF*Anv+hzQ3cDiIrBgHwD>73Z_%t?U7GID)jv<_=DB#s zGxy`aUo$JezB=Z%#Otd|b*Jk;?f<8LCOzIl^Z&wLxi!Yp-_IRSPI2U@WLngFlZlCu zfk7V4Be4wK=?)BK(*BBW3cBn^{x4)uiJOuXF(>=L2L=uGCm#cs%6l=Snto_f_`X7+ toVl%H(HetAdqt0J*?y?yTiA?4YbGch=3?-i##8kki|Nt~%-~Q15db<6dUpT- literal 0 HcmV?d00001 diff --git a/tools/angsd/test-data/test.bam b/tools/angsd/test-data/test_chrX.bam similarity index 100% rename from tools/angsd/test-data/test.bam rename to tools/angsd/test-data/test_chrX.bam From ed8cc608273ef88ea9c9c1d697a916dacf809f21 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 9 Aug 2026 11:40:31 +0000 Subject: [PATCH 5/6] fix: symlink fasta files --- tools/angsd/angsd.xml | 13 +- tools/angsd/estError.R | 620 ++++++++++++++++++++--------------------- 2 files changed, 314 insertions(+), 319 deletions(-) diff --git a/tools/angsd/angsd.xml b/tools/angsd/angsd.xml index c80e88359c3..a59a1951402 100644 --- a/tools/angsd/angsd.xml +++ b/tools/angsd/angsd.xml @@ -25,6 +25,13 @@ echo 'sample_${i}.bam' >> bam.filelist && #end for + #if str($summary_type.summary_type_selector) == 'error_est_out' + ln -f -s '$summary_type.ref' ref.fa && + ln -f -s '$summary_type.anc' anc.fa && + #elif str($summary_type.summary_type_selector) == 'mismatch' + ln -f -s '$summary_type.ref' ref.fa && + #end if + ## Essential command angsd -bam bam.filelist -out output -nThreads "\${GALAXY_SLOTS:-1}" @@ -53,8 +60,8 @@ #elif $summary_type.summary_type_selector == 'error_est_out': @BAM_PARAM@ -doAncError $summary_type.doAncError - -ref '$summary_type.ref' - -anc '$summary_type.anc' + -ref ref.fa + -anc anc.fa && Rscript '$__tool_directory__/estError.R' file=output.ancError > output.estError.txt #elif $summary_type.summary_type_selector == 'allele_counts': @@ -76,7 +83,7 @@ #elif $summary_type.summary_type_selector == 'mismatch': -doMisMatch 1 @BAM_PARAM@ - -ref '$summary_type.ref' + -ref ref.fa #end if ]]> diff --git a/tools/angsd/estError.R b/tools/angsd/estError.R index 57430a24d22..90615ab50ea 100644 --- a/tools/angsd/estError.R +++ b/tools/angsd/estError.R @@ -1,316 +1,304 @@ -bases <- c("A", "C", "G", "T") -b <- c(bases, "N") -########### do not change ################3 -l <- commandArgs(TRUE) -getArgs <- function(x, l) { - unlist(strsplit(grep(paste("^", x, "=", sep = ""), l, val = T), "="))[2] -} -Args <- function(l, args) { - if (!all(sapply(strsplit(l, "="), function(x) x[1]) %in% names(args))) { - cat("Error -> ", l[!sapply(strsplit(l, "="), function(x) x[1]) %in% names(args)], " is not a valid argument") - q("no") - } - arguments <- list() - for (a in names(args)) { - arguments[[a]] <- getArgs(a, l) - } - - if (any(!names(args) %in% names(arguments) & sapply(args, is.null))) { - cat("Error -> ", names(args)[!names(args) %in% names(arguments) & sapply(args, is.null)], " is not optional!\n") - q("no") - } - for (a in names(args)) { - if (is.null(arguments[[a]])) { - arguments[[a]] <- args[[match(a, names(args))]] - } - } - - - arguments -} - -print.args <- function(args, des) { - if (missing(des)) { - des <- as.list(rep("", length(args))) - names(des) <- names(args) - } - cat("-> Needed arguments:\n") - mapply(function(x) cat("\t", x, ":", des[[x]], "\n"), cbind(names(args)[sapply(args, is.null)])) - cat("-> Optional arguments (defaults):\n") - mapply(function(x) cat("\t", x, " (", args[[x]], ")", ":", des[[x]], "\n"), cbind(names(args)[!sapply(args, is.null)])) - q("no") -} - -## choose your parameters and defaults -## NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments -args <- list( - file = NULL, - out = "errorEst", - indNames = "ind", - nIter = 100, - subset = NA, - main = "Error rate using an outgroup and a high quality genome", - height = 7, - width = 12, - srt = 90, - cex = 1, - doPng = FALSE, - ylim1 = NA, - ylim2 = NA -) -# if no argument aree given prints the need arguments and the optional ones with default -des <- list( - file = "the ancError File", - out = "Name of the out files", - indNames = "postFix, file with names or comma seperated names of individuals", - nIter = "Numer of optimazation attemps", - subset = "comma seperated numbers of the individuals to include", - width = "width of the pdf", - height = "height of the pdf", - srt = "angle of ind names", - cex = "scale of names", - doPng = "Make png instead of pdf", - ylim1 = "Set max ylim for barplot1", - ylim2 = "Set max ylim for barplot2" -) - -####### get arguments and add to workspace -### do not change -if (length(l) == 0) print.args(args, des) -attach(Args(l, args)) -args <- commandArgs(TRUE) -if (length(args) == 0) { - cat(" Arguments: output prefix\n") - q("no") -} - -################################### -tol21rainbow <- c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#117744", "#44AA77", "#88CCAA", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455", "#DD7788") - - -width <- as.numeric(width) -height <- as.numeric(height) - -if (!is.na(ylim1)) { - ylim1 <- as.numeric(ylim1) -} -if (!is.na(ylim2)) { - ylim2 <- as.numeric(ylim2) -} - - -palette(tol21rainbow[c(1:8, 10:21)]) -srt <- as.numeric(srt) - - -nIter <- as.integer(nIter) -cex <- as.numeric(cex) -cat("----------\nfile: ", file, " out: ", out, " nIter: ", nIter, " subset: ", subset, "\n width: ", width, " height: ", height, " srt: ", srt, " cex: ", cex, " doPng: ", doPng, "\n-----------\n") - -r <- as.matrix(read.table(file)) - -if (!is.na(subset)) { - subset <- as.integer(unlist(strsplit(subset, ","))) - print(subset) - r <- r[subset, ] -} - -nInd <- nrow(r) -cat("Number of individuals read:", nInd, "\n") -{ - if (length(grep(",", indNames)) > 0) { - indNames <- unlist(strsplit(indNames, ",")) - } else { - options("warn" = -1) - try(indNames <- basename(scan(indNames, what = "theFuck")), silent = TRUE) - options("warn" = 0) - } -} - -if (length(indNames) == 1 & nInd > 1) { - indNames <- paste(indNames, 1:nInd, sep = "") -} -cat("Ind names:\n") -print(indNames) - -if (length(indNames) != nInd) { - cat("Error: Wrong number of ind Names\n") - q("no") -} - -getMat <- function(x) { - m <- array(0, dim = c(5, 5, 5), dimnames = list(b, b, b)) - for (s in 0:4) { - for (p in 0:4) { - for (a in 0:4) { - m[a + 1, p + 1, s + 1] <- x[a * 25 + p * 5 + s + 1] - } - } - } - m -} - -logLike <- function(x, Xch, Pch) { - eMat <- matrix(0, 4, 4) - eMat[-c(1, 6, 11, 16)] <- x - diag(eMat) <- 1 - rowSums(eMat) - P <- Pch %*% eMat - ll <- -sum(log(P) * Xch) - return(ll) -} - -res <- NULL -for (j in 1:nInd) { - m <- getMat(r[j, ]) - print(m) - - ## remove if missing - m <- m[-5, -5, -5] - - Pch <- matrix(0, 4, 4) - for (i in 1:4) { - Pch <- Pch + m[, , i] - } - print(Pch) - Pch <- Pch / rowSums(Pch) - cat("Matrix of allele change probabilities.\nA diagonal close to 1 denotes a low probability of errors.\n\n") - print(Pch) - - Xch <- matrix(0, 4, 4) - for (i in 1:4) { - Xch <- Xch + m[, i, ] - } - cat("\n\nObserved allele changes:\n\n") - print(Xch) - maxErr <- 0.02 - - conv <- list(par = 1, objective = Inf) - - for (i in 1:nIter) { - Tempconv <- nlminb(runif(12) / 100, logLike, upper = rep(maxErr, 12), lower = rep(1e-10, 12), Xch = Xch, Pch = Pch) - if (any(Tempconv$par > maxErr * 0.99)) { - maxErr <- min(1, maxErr * 10) - } - if (Tempconv$objective < conv$objective) { - conv <- Tempconv - } - } - - res <- rbind(res, conv$par) -} - -getover <- function(r, nInd) { - over <- NULL - for (j in 1:nInd) { - m <- getMat(r[j, ]) - - ## remove if missing - m <- m[-5, -5, -5] - Xch <- matrix(0, 4, 4) - for (i in 1:4) { - Xch <- Xch + m[, i, ] - } - - Pch <- matrix(0, 4, 4) - for (i in 1:4) { - Pch <- Pch + m[, , i] - } - - N1dot <- sum(Xch) - sum(diag(Xch)) - Ndot1 <- sum(Pch) - sum(diag(Pch)) - Ndot0 <- sum(diag(Pch)) - err <- (N1dot - Ndot1) / (Ndot0 - Ndot1) - over <- c(over, err) - } - over -} -over <- getover(r, nInd) -nam <- paste(rep(bases, 4), "->", rep(bases, each = 4))[-c(1, 6, 11, 16)] - -{ - if (doPng == FALSE) { - pdf(paste(out, ".pdf", sep = ""), width = width, height = height) - } else { - bitmap(paste(out, ".png", sep = ""), w = width, h = height, res = 300) - } -} - -{ - if (is.na(ylim1)) { - barplot(res, beside = T, col = 1:nInd, names = nam, main = main, ylab = "error rate") - } else { - barplot(res, beside = T, col = 1:nInd, names = nam, main = main, ylab = "error rate", ylim = c(0, ylim1)) - } -} -legend("top", paste(indNames, round(over * 100, 2), "%"), fill = 1:nInd, bty = "n") -dev.off() - -colnames(res) <- nam -rownames(res) <- indNames - -write.table(res, file = paste(out, ".txt", sep = ""), sep = "\t") -{ - if (doPng == FALSE) { - pdf(paste(out, "Overall.pdf", sep = ""), w = width, h = height) - } else { - bitmap(paste(out, "Overall.png", sep = ""), w = width, h = height, res = 300) - } -} - -par(mar = par()$mar + c(3, 0, 0, 0)) -{ - if (is.na(ylim2)) { - h <- barplot(over, col = 1:nInd, names = NULL, main = main, ylab = "error rate") - } else { - h <- barplot(over, col = 1:nInd, names = NULL, main = main, ylab = "error rate", ylim = c(0, ylim2)) - } -} -text(h, rep(-max(over) / 50, length(h)), indNames, xpd = T, srt = srt, adj = 1, cex = cex) -dev.off() -{ - if (doPng == FALSE) { - cat("figure:", paste(out, ".pdf", sep = ""), "\n") - cat("figure2:", paste(out, "Overall.pdf", sep = ""), "\n") - } else { - cat("figure:", paste(out, ".png", sep = ""), "\n") - cat("figure2:", paste(out, "Overall.png", sep = ""), "\n") - } -} -cat("table:", paste(out, ".txt", sep = ""), "\n") - -write.table(cbind(paste(indNames, round(over * 100, 4), "%")), file = paste(out, ".txt", sep = ""), sep = "\t", append = T, col = F, row = F) - -chrFile <- paste(file, "Chr", sep = "") - -if (file.info(chrFile)$size != 0 & !is.na(file.info(chrFile)$size)) { - overList <- list() - con <- file(chrFile, "r") - while (length(cchr <- scan(con, nlines = 1, what = "theFuck")) > 0) { - r <- read.table(con, nrow = nInd) - r <- as.matrix(r) - overList[[cchr[2]]] <- getover(r, nInd) - } - close(con) - - rbind.list <- function(l) { - n <- length(l) - results <- NULL - for (i in 1:n) { - results <- rbind(results, l[[i]]) - } - results - } - res <- rbind.list(overList) - chrNames <- names(overList) - rownames(res) <- chrNames - colnames(res) <- indNames - pdf(paste(out, "OverChr.pdf", sep = ""), w = width, h = height) - for (i in 1:nInd) { - dotchart(res[, i], chrNames, xlab = "Error rates", col = 3, main = indNames[i]) - dotchart(res[res[, i] > 0, i], chrNames, xlab = "Error rates (positive)", col = 3, main = indNames[i]) - } - dev.off() - write.table(res, file = paste(out, "OverChr.txt", sep = ""), qu = F, row = T, col = T) -} - -# warnings() +# nolint start +bases<-c("A","C","G","T") +b <- c(bases,"N") +########### do not change ################3 +l<-commandArgs(TRUE) +getArgs<-function(x,l) + unlist(strsplit(grep(paste("^",x,"=",sep=""),l,val=T),"="))[2] +Args<-function(l,args){ + if(! all(sapply(strsplit(l,"="),function(x)x[1])%in%names(args))){ + cat("Error -> ",l[!sapply(strsplit(l,"="),function(x)x[1])%in%names(args)]," is not a valid argument") + q("no") + } + arguments<-list() + for(a in names(args)) + arguments[[a]]<-getArgs(a,l) + + if(any(!names(args)%in%names(arguments)&sapply(args,is.null))){ + cat("Error -> ",names(args)[!names(args)%in%names(arguments)&sapply(args,is.null)]," is not optional!\n") + q("no") + } + for(a in names(args)) + if(is.null(arguments[[a]])) + arguments[[a]]<-args[[match(a,names(args))]] + + + arguments +} + +print.args<-function(args,des){ + if(missing(des)){ + des<-as.list(rep("",length(args))) + names(des)<-names(args) + } + cat("-> Needed arguments:\n") + mapply(function(x)cat("\t",x,":",des[[x]],"\n"),cbind(names(args)[sapply(args,is.null)])) + cat("-> Optional arguments (defaults):\n") + mapply(function(x)cat("\t",x," (",args[[x]],")",":",des[[x]],"\n"),cbind(names(args)[!sapply(args,is.null)])) + q("no") +} + +## choose your parameters and defaults +## NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments +args<-list(file=NULL, + out="errorEst", + indNames="ind", + nIter=100, + subset=NA, + main="Error rate using an outgroup and a high quality genome", + height=7, + width=12, + srt=90, + cex=1, + doPng=FALSE, + ylim1=NA, + ylim2=NA + ) +#if no argument aree given prints the need arguments and the optional ones with default +des<-list(file="the ancError File", + out="Name of the out files", + indNames="postFix, file with names or comma seperated names of individuals", + nIter="Numer of optimazation attemps", + subset="comma seperated numbers of the individuals to include", + width="width of the pdf", + height="height of the pdf", + srt="angle of ind names", + cex="scale of names", + doPng="Make png instead of pdf", + ylim1="Set max ylim for barplot1", + ylim2="Set max ylim for barplot2" + ) + +#######get arguments and add to workspace +### do not change +if(length(l)==0) print.args(args,des) +attach(Args(l,args)) +args <- commandArgs(TRUE) +if(length(args)==0){ + cat(" Arguments: output prefix\n") + q("no") +} + +################################### +tol21rainbow= c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#117744", "#44AA77", "#88CCAA", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455", "#DD7788") + + + +width=as.numeric(width) +height=as.numeric(height) + +if(!is.na(ylim1)) + ylim1=as.numeric(ylim1) +if(!is.na(ylim2)) + ylim2=as.numeric(ylim2) + + +palette(tol21rainbow[c(1:8,10:21)]) +srt<-as.numeric(srt) + + +nIter<-as.integer(nIter) +cex=as.numeric(cex) +cat("----------\nfile: ",file," out: ",out," nIter: ",nIter," subset: ",subset,"\n width: ",width," height: ",height," srt: ",srt," cex: ",cex," doPng: ",doPng,"\n-----------\n") + +r<-as.matrix(read.table(file)) + +if(!is.na(subset)){ + subset<- as.integer(unlist(strsplit(subset,","))) + print(subset) + r<-r[subset,] +} + +nInd<-nrow(r) +cat("Number of individuals read:",nInd,"\n") +{ + if(length(grep(",",indNames))>0) + indNames<-unlist(strsplit(indNames,",")) + else { + options("warn"=-1) + try(indNames<-basename(scan(indNames,what="theFuck")),silent=TRUE) + options("warn"=0) + } +} + +if(length(indNames)==1&nInd>1){ + indNames<-paste(indNames,1:nInd,sep="") +} +cat("Ind names:\n") +print(indNames) + +if(length(indNames)!=nInd){ + cat("Error: Wrong number of ind Names\n") + q("no") +} + +getMat<-function(x){ + m<-array(0,dim=c(5,5,5),dimnames=list(b,b,b)) + for(s in 0:4) + for(p in 0:4) + for(a in 0:4) + m[a+1,p+1,s+1]<-x[a*25+p*5+s+1] + m +} + +logLike<-function(x,Xch,Pch){ + eMat<-matrix(0,4,4) + eMat[-c(1,6,11,16)]<-x + diag(eMat)<-1-rowSums(eMat) + P <- Pch %*% eMat + ll <- -sum(log(P)*Xch) + return(ll) +} + +res<-NULL +for(j in 1:nInd){ + m<-getMat(r[j,]) + print(m) + + ##remove if missing + m<-m[-5,-5,-5] + + Pch<-matrix(0,4,4) + for(i in 1:4) + Pch<-Pch+m[,,i] + print(Pch) + Pch<-Pch/rowSums(Pch) + cat("Matrix of allele change probabilities.\nA diagonal close to 1 denotes a low probability of errors.\n\n") + print(Pch) + + Xch<-matrix(0,4,4) + for(i in 1:4) + Xch<-Xch+m[,i,] + cat("\n\nObserved allele changes:\n\n") + print(Xch) + maxErr=0.02; + + conv<-list(par=1,objective=Inf) + + for(i in 1:nIter){ + Tempconv <- nlminb(runif(12)/100,logLike,upper=rep(maxErr,12),lower=rep(1e-10,12),Xch=Xch,Pch=Pch) + if( any(Tempconv$par > maxErr*0.99) ) + maxErr<-min(1,maxErr*10); + if(Tempconv$objective",rep(bases,each=4))[-c(1,6,11,16)] + +{ + if(doPng==FALSE){ + pdf(paste(out,".pdf",sep=""),width=width,height=height) + } else + bitmap(paste(out,".png",sep=""),w=width,h=height,res=300) +} + +{ +if(is.na(ylim1)) + barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate") +else + barplot(res,beside=T,col=1:nInd,names=nam,main=main,ylab="error rate",ylim=c(0,ylim1)) +} +legend("top",paste(indNames,round(over*100,2),"%"),fill=1:nInd,bty="n") +dev.off() + +colnames(res)<-nam +rownames(res)<-indNames + +write.table(res,file=paste(out,".txt",sep=""),sep="\t") +{ +if(doPng==FALSE) + pdf(paste(out,"Overall.pdf",sep=""),w=width,h=height) +else + bitmap(paste(out,"Overall.png",sep=""),w=width,h=height,res=300) +} + +par(mar=par()$mar+c(3,0,0,0)) +{ +if(is.na(ylim2)){ + h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate") +}else{ + h<-barplot(over,col=1:nInd,names=NULL,main=main,ylab="error rate",ylim=c(0,ylim2)) +} +} +text(h,rep(-max(over)/50,length(h)),indNames,xpd=T,srt=srt,adj=1,cex=cex) +dev.off() +{ + if(doPng==FALSE){ + cat("figure:",paste(out,".pdf",sep=""),"\n") + cat("figure2:",paste(out,"Overall.pdf",sep=""),"\n") + }else{ + cat("figure:",paste(out,".png",sep=""),"\n") + cat("figure2:",paste(out,"Overall.png",sep=""),"\n") + } + +} +cat("table:",paste(out,".txt",sep=""),"\n") + +write.table(cbind(paste(indNames,round(over*100,4),"%")),file=paste(out,".txt",sep=""),sep="\t",append=T,col=F,row=F) + +chrFile<-paste(file,"Chr",sep="") + +if(file.info(chrFile)$size!=0 & !is.na(file.info(chrFile)$size)){ + + overList<-list() + con<-file(chrFile,"r") + while(length(cchr<-scan(con,nlines=1,what="theFuck"))>0){ + r<-read.table(con,nrow=nInd) + r<-as.matrix(r) + overList[[cchr[2]]]<-getover(r,nInd) + } + close(con) + + rbind.list<-function (l) + { + n <- length(l) + results <- NULL + for (i in 1:n) { + results <- rbind(results, l[[i]]) + } + results + } + res<-rbind.list(overList) + chrNames<-names(overList) + rownames(res) <- chrNames + colnames(res) <- indNames + pdf(paste(out,"OverChr.pdf",sep=""),w=width,h=height) + for(i in 1:nInd){ + dotchart(res[,i],chrNames,xlab="Error rates",col=3,main=indNames[i]) + dotchart(res[res[,i]>0,i],chrNames,xlab="Error rates (positive)",col=3,main=indNames[i]) + } + dev.off() + write.table(res,file=paste(out,"OverChr.txt",sep=""),qu=F,row=T,col=T) +} + +#warnings() +# nolint end \ No newline at end of file From 229fc6c8ecb5c6a26c14990000766ec46f76e1f7 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 9 Aug 2026 12:10:35 +0000 Subject: [PATCH 6/6] fix: force faidx --- tools/angsd/angsd.xml | 3 +++ tools/angsd/estError.R | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/tools/angsd/angsd.xml b/tools/angsd/angsd.xml index a59a1951402..6752fa1d07e 100644 --- a/tools/angsd/angsd.xml +++ b/tools/angsd/angsd.xml @@ -28,8 +28,11 @@ #if str($summary_type.summary_type_selector) == 'error_est_out' ln -f -s '$summary_type.ref' ref.fa && ln -f -s '$summary_type.anc' anc.fa && + samtools faidx ref.fa && + samtools faidx anc.fa && #elif str($summary_type.summary_type_selector) == 'mismatch' ln -f -s '$summary_type.ref' ref.fa && + samtools faidx ref.fa && #end if ## Essential command diff --git a/tools/angsd/estError.R b/tools/angsd/estError.R index 90615ab50ea..ddc50c065b4 100644 --- a/tools/angsd/estError.R +++ b/tools/angsd/estError.R @@ -301,4 +301,4 @@ if(file.info(chrFile)$size!=0 & !is.na(file.info(chrFile)$size)){ } #warnings() -# nolint end \ No newline at end of file +# nolint end