import%20marimo%0A%0A__generated_with%20%3D%20%220.20.2%22%0Aapp%20%3D%20marimo.App(width%3D%22full%22%2C%20app_title%3D%22Sheep%20Population%20Genetics%20Analysis%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Saudi%20Sheep%20Population%20Genetics%20Analysis%0A%20%20%20%20%23%23%23%20OvineSNP50%20whole-genome%20genotyping%20Array%20%E2%80%94%20*Ovis%20aries*%20(Rambouillet%20v1.0%20Reference)%0A%0A%20%20%20%20This%20notebook%20documents%20a%20complete%20population%20genetics%20pipeline%20for%20the%0A%20%20%20%20**OvineSNP50%20whole-genome%20genotyping%20array%20dataset%2C%20following%20the%20methodology%20of%20%5BComparative%20genome-wide%20analysis%20of%20Ovis%20aries%20in%20Saudi%20Arabia%2C%20highlighting%20inbreeding%20and%20genetic%20isolation%20of%20the%20Najdi%20sheep%20breed%5D(https%3A%2F%2Fwww.frontiersin.org%2Fjournals%2Fgenetics%2Farticles%2F10.3389%2Ffgene.2025.1646127%2Ffull%23s2)%0A%20%20%20%20study.%20PLINK%201.9%20performs%20all%20genotype-level%20computations.%0A%20%20%20%20Python%20handles%20post-processing%20and%20visualisation.%0A%0A%20%20%20%20%7C%20Attribute%20%7C%20Value%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20VCF%20%7C%20%60OvineSNPIII_for_EVA.submit.vcf.gz%60%20%7C%0A%20%20%20%20%7C%20Raw%20SNPs%20%7C%2053%2C370%20%7C%0A%20%20%20%20%7C%20Samples%20%7C%2096%20%7C%0A%20%20%20%20%7C%20Breeds%20%7C%20Harri%20(n%20%3D%2039)%20%C2%B7%20Naemi%20(n%20%3D%2024)%20%C2%B7%20Najdi%20(n%20%3D%2033)%20%7C%0A%20%20%20%20%7C%20Reference%20genome%20%7C%20*Ovis%20aries*%20Rambouillet%20v1.0%20(GCF_002742125.1)%20%7C%0A%0A%20%20%20%20**Pipeline%20%E2%80%94%20aligned%20with%20reference%20study**%0A%0A%20%20%20%20%7C%20%20%7C%20Step%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%201b%20%7C%20VCF%20inspection%20%7C%0A%20%20%20%20%7C%202%20%7C%20Chromosome%20renaming%20%7C%0A%20%20%20%20%7C%203%20%7C%20Breed%20assignment%20%7C%0A%20%20%20%20%7C%204%20%7C%20VCF%20%E2%86%92%20PLINK%20binary%20%7C%0A%20%20%20%20%7C%205%20%7C%20QC%20filtering%20%2B%20heterozygosity%20outlier%20detection%20%7C%0A%20%20%20%20%7C%206%20%7C%20LD%20pruning%20(r%C2%B2%20%3D%200.3)%20%7C%0A%20%20%20%20%7C%207%20%7C%20MAF%20distribution%20per%20breed%20%7C%0A%20%20%20%20%7C%208%20%7C%20PCA%20%7C%0A%20%20%20%20%7C%209%20%7C%20Population%20structure%20%E2%80%94%20NMF%20K%20%3D%201%E2%80%936%20with%20cross-entropy%20%7C%0A%20%20%20%20%7C%209b%20%7C%20Population%20structure%20%E2%80%94%20ADMIXTURE%20K%20%3D%201%E2%80%936%20(alternative)%20%7C%0A%20%20%20%20%7C%2010%20%7C%20Genetic%20diversity%20%E2%80%94%20Ho%2C%20He%2C%20F_IS%2C%20Ne%20trajectory%20%7C%0A%20%20%20%20%7C%2011%20%7C%20Pairwise%20F_ST%20%2B%20Manhattan%20plot%20%7C%0A%20%20%20%20%7C%2012%20%7C%20Runs%20of%20Homozygosity%20%E2%80%94%20F_ROH%20%2B%20robustness%20tests%20%7C%0A%20%20%20%20%7C%2013%20%7C%20Per-chromosome%20SNP%20density%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%201%20%C2%B7%20Setup%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20import%20re%0A%20%20%20%20import%20subprocess%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20matplotlib.patches%20as%20mpatches%0A%20%20%20%20import%20seaborn%20as%20sns%0A%0A%20%20%20%20PROJECT%20%3D%20Path(%22%2FUsers%2Ftalal.alyazeedi%2FDesktop%2FSheep_project%22)%0A%20%20%20%20VCF%20%20%20%20%20%3D%20PROJECT%20%2F%20%22OvineSNPIII_for_EVA.submit.vcf.gz%22%0A%20%20%20%20OUT%20%20%20%20%20%3D%20PROJECT%20%2F%20%22plink_analysis%22%0A%20%20%20%20PLINK%20%20%20%3D%20%22%2FUsers%2Ftalal.alyazeedi%2Fminiconda3%2Fbin%2Fplink%22%0A%20%20%20%20OUT.mkdir(exist_ok%3DTrue)%0A%0A%20%20%20%20mo.callout(mo.md(%0A%20%20%20%20%20%20%20%20f%22**Project%3A**%20%60%7BPROJECT%7D%60%20%20%5Cn%22%0A%20%20%20%20%20%20%20%20f%22**Output%3A**%20%60%7BOUT%7D%60%20%20%5Cn%22%0A%20%20%20%20%20%20%20%20f%22**PLINK%3A**%20%60%7BPLINK%7D%60%22%0A%20%20%20%20)%2C%20kind%3D%22info%22)%0A%20%20%20%20return%20OUT%2C%20PLINK%2C%20Path%2C%20VCF%2C%20mpatches%2C%20np%2C%20pd%2C%20plt%2C%20re%2C%20sns%2C%20subprocess%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%201b%20%C2%B7%20VCF%20Inspection%0A%0A%20%20%20%20Before%20any%20analysis%2C%20we%20inspect%20the%20raw%20VCF%20to%20confirm%20sample%20count%2C%0A%20%20%20%20variant%20count%2C%20variant%20types%2C%20and%20chromosome%20naming%20convention.%0A%0A%20%20%20%20%7C%20Command%20%7C%20Purpose%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20%60bcftools%20stats%60%20%7C%20Summary%20statistics%20%E2%80%94%20N%20variants%2C%20N%20samples%2C%20ts%2Ftv%20ratio%20%7C%0A%20%20%20%20%7C%20%60bcftools%20query%20-l%60%20%7C%20List%20all%20sample%20IDs%20present%20in%20the%20file%20%7C%0A%20%20%20%20%7C%20%60bcftools%20view%20--header-only%60%20%7C%20Inspect%20contig%20names%20and%20array%20metadata%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(VCF%2C%20mo%2C%20pd%2C%20subprocess)%3A%0A%20%20%20%20_cmd_stats%20%3D%20%22bcftools%20stats%20OvineSNPIII_for_EVA.submit.vcf.gz%20%7C%20grep%20'%5ESN'%22%0A%20%20%20%20_cmd_samps%20%3D%20%22bcftools%20query%20-l%20OvineSNPIII_for_EVA.submit.vcf.gz%22%0A%20%20%20%20_cmd_ctg%20%20%20%3D%20%22bcftools%20view%20--header-only%20OvineSNPIII_for_EVA.submit.vcf.gz%20%7C%20grep%20'%5E%23%23contig'%22%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20bcftools%20stats%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_r_st%20%20%20%20%3D%20subprocess.run(%5B%22bcftools%22%2C%20%22stats%22%2C%20str(VCF)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2C%20text%3DTrue)%0A%20%20%20%20_sn_rows%20%3D%20%5Bl.split(%22%5Ct%22)%20for%20l%20in%20_r_st.stdout.splitlines()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20l.startswith(%22SN%22)%5D%0A%20%20%20%20_stats_df%20%3D%20pd.DataFrame(%0A%20%20%20%20%20%20%20%20%5B(r%5B2%5D.rstrip(%22%3A%22)%2C%20r%5B3%5D.strip())%20for%20r%20in%20_sn_rows%5D%2C%0A%20%20%20%20%20%20%20%20columns%3D%5B%22Statistic%22%2C%20%22Value%22%5D%2C%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20sample%20list%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_r_ql%20%20%3D%20subprocess.run(%5B%22bcftools%22%2C%20%22query%22%2C%20%22-l%22%2C%20str(VCF)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2C%20text%3DTrue)%0A%20%20%20%20_samps%20%3D%20_r_ql.stdout.strip().splitlines()%0A%20%20%20%20_samp_df%20%3D%20pd.DataFrame(%7B%22Sample%20ID%22%3A%20_samps%7D)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20contig%20names%20from%20header%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_r_hdr%20%3D%20subprocess.run(%5B%22bcftools%22%2C%20%22view%22%2C%20%22--header-only%22%2C%20str(VCF)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2C%20text%3DTrue)%0A%20%20%20%20_ctgs%20%20%3D%20%5Bl.split(%22ID%3D%22)%5B1%5D.rstrip(%22%3E%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20l%20in%20_r_hdr.stdout.splitlines()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20l.startswith(%22%23%23contig%22)%5D%0A%20%20%20%20_ctg_df%20%3D%20pd.DataFrame(%7B%22Contig%20(NCBI%20accession)%22%3A%20_ctgs%7D)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20summary%20callout%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_n_snp%20%3D%20next((r%5B3%5D.strip()%20for%20r%20in%20_sn_rows%20if%20%22number%20of%20SNPs%22%20in%20r%5B2%5D)%2C%20%22%3F%22)%0A%20%20%20%20_n_sam%20%3D%20next((r%5B3%5D.strip()%20for%20r%20in%20_sn_rows%20if%20%22number%20of%20samples%22%20in%20r%5B2%5D)%2C%20%22%3F%22)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_stats%7D%5Cn%5Cn%7B_cmd_samps%7D%5Cn%5Cn%7B_cmd_ctg%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22**%7B_n_sam%7D%20samples**%20%C2%B7%20**%7B_n_snp%7D%20SNPs**%20%C2%B7%20all%20biallelic%20%C2%B7%20no%20indels%22%0A%20%20%20%20%20%20%20%20)%2C%20kind%3D%22success%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20%60bcftools%20stats%60%20summary%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_stats_df)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Sample%20IDs%20(%60bcftools%20query%20-l%60)%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22Prefix%20**H**%20%E2%86%92%20Harri%20(n%20%3D%2039)%20%20%C2%B7%20%20**M**%20%E2%86%92%20Naemi%20(n%20%3D%2024)%20%20%C2%B7%20%20**N**%20%E2%86%92%20Najdi%20(n%20%3D%2033)%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_samp_df)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Chromosome%20names%20in%20raw%20VCF%20(%60%23%23contig%60%20header%20lines)%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22NCBI%20RefSeq%20accessions%20%E2%80%94%20these%20must%20be%20renamed%20to%20integers%201%E2%80%9323%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22before%20PLINK%20import%20(see%202).%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_ctg_df)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%202%20%C2%B7%20Chromosome%20Renaming%20Map%0A%0A%20%20%20%20The%20VCF%20carries%20NCBI%20RefSeq%20accessions%20(e.g.%20%60NC_040252.1%60)%20as%20chromosome%0A%20%20%20%20identifiers.%20%20%60bcftools%20annotate%20--rename-chrs%60%20replaces%20them%20with%20the%0A%20%20%20%20integers%201%E2%80%9323%20before%20PLINK%20import.%0A%0A%20%20%20%20%7C%20Accession%20%7C%20Chr%20%7C%20Notes%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20NC_040252.1%20%7C%201%20%7C%20%E2%80%A6%20%7C%0A%20%20%20%20%7C%20NC_040273.1%20%7C%2022%20%7C%20last%20autosome%20present%20on%20this%20array%20%7C%0A%20%20%20%20%7C%20NC_040278.1%20%7C%2023%20%7C%20Chr%20X%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20mo%2C%20pd)%3A%0A%20%20%20%20_map%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22NC_040252.1%22%3A%221%22%2C%20%20%22NC_040253.1%22%3A%222%22%2C%20%20%22NC_040254.1%22%3A%223%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040255.1%22%3A%224%22%2C%20%20%22NC_040256.1%22%3A%225%22%2C%20%20%22NC_040257.1%22%3A%226%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040258.1%22%3A%227%22%2C%20%20%22NC_040259.1%22%3A%228%22%2C%20%20%22NC_040260.1%22%3A%229%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040261.1%22%3A%2210%22%2C%20%22NC_040262.1%22%3A%2211%22%2C%20%22NC_040263.1%22%3A%2212%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040264.1%22%3A%2213%22%2C%20%22NC_040265.1%22%3A%2214%22%2C%20%22NC_040266.1%22%3A%2215%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040267.1%22%3A%2216%22%2C%20%22NC_040268.1%22%3A%2217%22%2C%20%22NC_040269.1%22%3A%2218%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040270.1%22%3A%2219%22%2C%20%22NC_040271.1%22%3A%2220%22%2C%20%22NC_040272.1%22%3A%2221%22%2C%0A%20%20%20%20%20%20%20%20%22NC_040273.1%22%3A%2222%22%2C%20%22NC_040278.1%22%3A%2223%22%2C%0A%20%20%20%20%7D%0A%20%20%20%20chr_map_file%20%3D%20OUT%20%2F%20%22chr_map.txt%22%0A%20%20%20%20with%20open(chr_map_file%2C%20%22w%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20for%20k%2C%20v%20in%20_map.items()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_f.write(f%22%7Bk%7D%5Ct%7Bv%7D%5Cn%22)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(f%22Wrote%20%60%7Bchr_map_file.name%7D%60%20%E2%80%94%20%7Blen(_map)%7D%20entries%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(pd.DataFrame(list(_map.items())%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20columns%3D%5B%22NCBI%20accession%22%2C%20%22Numeric%22%5D))%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(chr_map_file%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%203%20%C2%B7%20Breed%20Assignment%0A%0A%20%20%20%20%7C%20Prefix%20%7C%20Breed%20%7C%20n%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20H%20%7C%20**Harri**%20%7C%2039%20%7C%0A%20%20%20%20%7C%20M%20%7C%20**Naemi**%20%7C%2024%20%7C%0A%20%20%20%20%7C%20N%20%7C%20**Najdi**%20%7C%2033%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20VCF%2C%20mo%2C%20pd%2C%20subprocess)%3A%0A%20%20%20%20_r%20%3D%20subprocess.run(%5B%22bcftools%22%2C%22view%22%2C%22--header-only%22%2Cstr(VCF)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2C%20text%3DTrue)%0A%20%20%20%20_hdr%20%3D%20%5Bl%20for%20l%20in%20_r.stdout.splitlines()%20if%20l.startswith(%22%23CHROM%22)%5D%5B-1%5D%0A%20%20%20%20_samples%20%3D%20_hdr.strip().split(%22%5Ct%22)%5B9%3A%5D%0A%0A%20%20%20%20pop_df%20%3D%20pd.DataFrame(%7B%0A%20%20%20%20%20%20%20%20%22FID%22%3A%20_samples%2C%20%22IID%22%3A%20_samples%2C%0A%20%20%20%20%20%20%20%20%22POP%22%3A%20%5B%7B%22H%22%3A%22Harri%22%2C%22M%22%3A%22Naemi%22%2C%22N%22%3A%22Najdi%22%7D.get(s%5B0%5D%2C%22%3F%22)%20for%20s%20in%20_samples%5D%2C%0A%20%20%20%20%7D)%0A%20%20%20%20_clust%20%3D%20OUT%20%2F%20%22sheep_clusters.txt%22%0A%20%20%20%20pop_df%5B%5B%22FID%22%2C%22IID%22%2C%22POP%22%5D%5D.to_csv(_clust%2C%20sep%3D%22%5Ct%22%2C%20index%3DFalse%2C%20header%3DFalse)%0A%20%20%20%20for%20_b%20in%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%3A%0A%20%20%20%20%20%20%20%20pop_df%5Bpop_df%5B%22POP%22%5D%3D%3D_b%5D%5B%5B%22FID%22%2C%22IID%22%5D%5D.to_csv(%0A%20%20%20%20%20%20%20%20%20%20%20%20OUT%2Ff%22keep_%7B_b%7D.txt%22%2C%20sep%3D%22%5Ct%22%2C%20index%3DFalse%2C%20header%3DFalse)%0A%0A%20%20%20%20_cnt%20%3D%20pop_df%5B%22POP%22%5D.value_counts().reset_index()%0A%20%20%20%20_cnt.columns%20%3D%20%5B%22Breed%22%2C%22n%22%5D%0A%20%20%20%20mo.vstack(%5Bmo.md(f%22**%7Blen(_samples)%7D%20samples**%22)%2C%20mo.ui.table(_cnt)%5D)%0A%20%20%20%20return%20(pop_df%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%204%20%C2%B7%20VCF%20%E2%86%92%20PLINK%20Binary%20Conversion%0A%0A%20%20%20%20Two%20commands%20convert%20the%20VCF%20to%20PLINK's%20efficient%20binary%20format%3A%0A%0A%20%20%20%201.%20%60bcftools%20annotate%20--rename-chrs%60%20%E2%80%94%20swap%20NCBI%20accessions%20for%20integers%0A%20%20%20%202.%20%60plink%20--vcf%20--double-id%20--chr-set%2026%20--make-bed%60%20%E2%80%94%20write%20%60.bed%2F.bim%2F.fam%60%0A%0A%20%20%20%20%60--chr-set%2026%60%20declares%2026%20autosomes%20(sheep%202n%20%3D%2054).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20Path%2C%20VCF%2C%20chr_map_file%2C%20mo%2C%20subprocess)%3A%0A%20%20%20%20_vcf2%20%20%3D%20OUT%20%2F%20%22sheep_renamed.vcf.gz%22%0A%20%20%20%20raw_prefix%20%3D%20str(OUT%20%2F%20%22sheep_raw%22)%0A%20%20%20%20_cmd%20%3D%20(%0A%20%20%20%20%20%20%20%20%22bcftools%20annotate%20--rename-chrs%20chr_map.txt%20-Oz%20-o%20sheep_renamed.vcf.gz%20input.vcf.gz%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--vcf%20sheep_renamed.vcf.gz%20%5C%5C%5Cn%20%20%20%20--double-id%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--make-bed%20%5C%5C%5Cn%20%20%20%20--out%20sheep_raw%22%0A%20%20%20%20)%0A%20%20%20%20if%20not%20_vcf2.exists()%3A%0A%20%20%20%20%20%20%20%20_r%20%3D%20subprocess.run(%5B%22bcftools%22%2C%22annotate%22%2C%22--rename-chrs%22%2Cstr(chr_map_file)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22-Oz%22%2C%22-o%22%2Cstr(_vcf2)%2Cstr(VCF)%5D%2Ccapture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20%20%20%20%20subprocess.run(%5B%22bcftools%22%2C%22index%22%2Cstr(_vcf2)%5D%2Ccapture_output%3DTrue)%0A%20%20%20%20if%20not%20Path(raw_prefix%2B%22.bed%22).exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--vcf%22%2Cstr(_vcf2)%2C%22--double-id%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--make-bed%22%2C%22--out%22%2Craw_prefix%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20_log%20%3D%20Path(raw_prefix%2B%22.log%22).read_text()%20if%20Path(raw_prefix%2B%22.log%22).exists()%20else%20%22%22%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20PLINK%20log%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60%5Cn%7B_log%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(raw_prefix%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%205%20%C2%B7%20Quality%20Control%20%2B%20Heterozygosity%20Outlier%20Detection%0A%0A%20%20%20%20%23%23%23%205a%20%E2%80%94%20SNP%20and%20individual%20filters%0A%0A%20%20%20%20%7C%20Flag%20%7C%20Threshold%20%7C%20Removes%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20%60--geno%60%20%7C%20%3E%205%20%25%20missing%20%7C%20Low-quality%20SNPs%20%7C%0A%20%20%20%20%7C%20%60--maf%60%20%7C%20%3C%201%20%25%20%7C%20Near-monomorphic%20SNPs%20%7C%0A%20%20%20%20%7C%20%60--mind%60%20%7C%20%3E%2010%20%25%20missing%20%7C%20Low-quality%20individuals%20%7C%0A%20%20%20%20%7C%20%60--hwe%60%20%7C%20p%20%3C%201%20%C3%97%2010%E2%81%BB%E2%81%B6%20%7C%20Genotyping-error%20SNPs%20(lenient%20threshold)%20%7C%0A%0A%20%20%20%20%23%23%23%205b%20%E2%80%94%20Heterozygosity%20outlier%20removal%0A%20%20%20%20Following%20the%20reference%20study%2C%20individuals%20whose%20observed%20heterozygosity%0A%20%20%20%20falls%20**%3E%203%20SD**%20from%20the%20breed-mean%20are%20flagged%20and%20removed.%20%20The%0A%20%20%20%20reference%20detected%20and%20removed%20**H16**%20from%20their%20dataset.%20%20We%20inspect%20our%0A%20%20%20%2096%20samples%20for%20the%20same%20pattern.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20Path%2C%20mo%2C%20pd%2C%20plt%2C%20raw_prefix%2C%20re%2C%20subprocess)%3A%0A%20%20%20%20qc_prefix%20%3D%20str(OUT%20%2F%20%22sheep_qc%22)%0A%20%20%20%20_cmd_qc%20%3D%20(%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_raw%20%5C%5C%5Cn%20%20%20%20--geno%200.05%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--maf%20%200.01%20%5C%5C%5Cn%20%20%20%20--mind%200.10%20%5C%5C%5Cn%20%20%20%20--hwe%20%201e-6%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--make-bed%20%5C%5C%5Cn%20%20%20%20--out%20sheep_qc%22%0A%20%20%20%20)%0A%20%20%20%20_cmd_het%20%3D%20(%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--het%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20sheep_het%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20step%201%3A%20QC%20filters%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20if%20not%20Path(qc_prefix%2B%22.bed%22).exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Craw_prefix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--geno%22%2C%220.05%22%2C%22--maf%22%2C%220.01%22%2C%22--mind%22%2C%220.10%22%2C%22--hwe%22%2C%221e-6%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--make-bed%22%2C%22--out%22%2Cqc_prefix%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20_log_qc%20%3D%20Path(qc_prefix%2B%22.log%22).read_text()%20if%20Path(qc_prefix%2B%22.log%22).exists()%20else%20%22%22%0A%20%20%20%20_pass%20%3D%20re.search(r%22(%5Cd%2B)%20variants%3F%20and%20(%5Cd%2B)%20people%20pass%22%2C%20_log_qc)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20step%202%3A%20heterozygosity%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20het_file%20%3D%20OUT%20%2F%20%22sheep_het.het%22%0A%20%20%20%20if%20not%20het_file.exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%22--het%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2Cstr(OUT%2F%22sheep_het%22)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20_het%20%3D%20pd.read_csv(het_file%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20_het.columns%20%3D%20%5Bc.strip()%20for%20c%20in%20_het.columns%5D%0A%20%20%20%20_het%5B%22Ho%22%5D%20%3D%20(_het%5B%22N(NM)%22%5D%20-%20_het%5B%22O(HOM)%22%5D)%20%2F%20_het%5B%22N(NM)%22%5D%0A%20%20%20%20_mn%2C%20_sd%20%3D%20_het%5B%22Ho%22%5D.mean()%2C%20_het%5B%22Ho%22%5D.std()%0A%20%20%20%20_outliers%20%3D%20_het%5B(_het%5B%22Ho%22%5D%20%3C%20_mn%20-%203*_sd)%20%7C%20(_het%5B%22Ho%22%5D%20%3E%20_mn%20%2B%203*_sd)%5D%0A%0A%20%20%20%20%23%20write%20remove%20file%20(may%20be%20empty)%0A%20%20%20%20_outliers%5B%5B%22FID%22%2C%22IID%22%5D%5D.to_csv(OUT%2F%22het_outliers.txt%22%2Csep%3D%22%5Ct%22%2Cindex%3DFalse%2Cheader%3DFalse)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20plot%20het%20distribution%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig0%2C%20_ax0%20%3D%20plt.subplots(figsize%3D(8%2C3))%0A%20%20%20%20_ax0.hist(_het%5B%22Ho%22%5D%2C%20bins%3D30%2C%20color%3D%22%23457B9D%22%2C%20edgecolor%3D%22white%22)%0A%20%20%20%20for%20_b%20in%20%5B_mn-3*_sd%2C%20_mn%2B3*_sd%5D%3A%0A%20%20%20%20%20%20%20%20_ax0.axvline(_b%2C%20color%3D%22%23E63946%22%2C%20linestyle%3D%22--%22%2C%20linewidth%3D1.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22%C2%B13%20SD%20(%7B_b%3A.3f%7D)%22)%0A%20%20%20%20_ax0.set_xlabel(%22Observed%20heterozygosity%20(Ho)%22)%0A%20%20%20%20_ax0.set_ylabel(%22Count%22)%0A%20%20%20%20_ax0.set_title(%22Heterozygosity%20distribution%20%E2%80%94%20outlier%20detection%22)%0A%20%20%20%20_ax0.legend(fontsize%3D9)%0A%20%20%20%20_ax0.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20_msg%20%3D%20(f%22No%20heterozygosity%20outliers%20detected%20(all%20samples%20within%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%5B%7B_mn-3*_sd%3A.4f%7D%2C%20%7B_mn%2B3*_sd%3A.4f%7D%5D).%22)%0A%20%20%20%20if%20len(_outliers)%3A%0A%20%20%20%20%20%20%20%20_msg%20%3D%20(f%22%7Blen(_outliers)%7D%20outlier(s)%20flagged%20for%20removal%3A%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%2B%20%22%2C%20%22.join(_outliers%5B%22IID%22%5D.tolist()))%0A%0A%20%20%20%20_summary%20%3D%20(f%22**After%20QC%3A**%20%7B_pass.group(1)%7D%20SNPs%20%C2%B7%20%7B_pass.group(2)%7D%20individuals%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20_pass%20else%20%22See%20log%22)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_qc%7D%5Cn%5Cn%7B_cmd_het%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(_summary)%2C%20kind%3D%22success%22)%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(f%22**Outlier%20check%3A**%20%7B_msg%7D%22)%2C%20kind%3D%22info%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig0)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60%5Cn%7B_log_qc%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20het_file%2C%20qc_prefix%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%206%20%C2%B7%20LD%20Pruning%20(r%C2%B2%20%3D%200.3)%0A%0A%20%20%20%20The%20reference%20study%20retained%20**34%2C079%20SNPs**%20after%20LD%20pruning%20of%20a%0A%20%20%20%2058%2C116-SNP%20dataset%20(58.6%20%25%20retention).%20%20Using%20our%2051%2C457-SNP%20dataset%0A%20%20%20%20with%20r%C2%B2%20%3D%200.3%20yields%20a%20comparable%20retention%20rate.%0A%0A%20%20%20%20%3E%20**Why%20r%C2%B2%20%3D%200.3%3F**%20%20This%20is%20the%20standard%20threshold%20in%20livestock%0A%20%20%20%20%3E%20population-structure%20analyses.%20%20Our%20previous%20r%C2%B2%20%3D%200.1%20was%20too%20strict%0A%20%20%20%20%3E%20(retaining%20only%2012%2C404%20SNPs%2C%20~24%20%25)%2C%20which%20over-corrects%20for%20LD%20and%0A%20%20%20%20%3E%20discards%20informative%20variation.%0A%0A%20%20%20%20%7C%20Parameter%20%7C%20Value%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20Window%20%7C%2050%20SNPs%20%7C%0A%20%20%20%20%7C%20Step%20%7C%2010%20SNPs%20%7C%0A%20%20%20%20%7C%20r%C2%B2%20threshold%20%7C%20**0.3**%20%7C%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20Path%2C%20mo%2C%20qc_prefix%2C%20re%2C%20subprocess)%3A%0A%20%20%20%20prune_prefix%20%3D%20str(OUT%20%2F%20%22sheep_pruned_r03%22)%0A%20%20%20%20_cmd1%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--indep-pairwise%2050%2010%200.3%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--out%20sheep_pruned_r03%22)%0A%20%20%20%20_cmd2%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--extract%20sheep_pruned_r03.prune.in%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--make-bed%20%5C%5C%5Cn%20%20%20%20--out%20sheep_pruned_r03%22)%0A%0A%20%20%20%20if%20not%20Path(prune_prefix%2B%22.prune.in%22).exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--indep-pairwise%22%2C%2250%22%2C%2210%22%2C%220.3%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2Cprune_prefix%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20if%20not%20Path(prune_prefix%2B%22.bed%22).exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--extract%22%2Cprune_prefix%2B%22.prune.in%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--make-bed%22%2C%22--out%22%2Cprune_prefix%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20_log%20%3D%20Path(prune_prefix%2B%22.log%22).read_text()%20if%20Path(prune_prefix%2B%22.log%22).exists()%20else%20%22%22%0A%20%20%20%20_rm%20%20%3D%20re.search(r%22Pruning%20complete%5C.%5Cs%2B(%5Cd%2B)%20of%20(%5Cd%2B)%20variants%3F%20removed%22%2C%20_log)%0A%20%20%20%20_kp%20%20%3D%20re.search(r%22(%5Cd%2B)%20variants%3F%20and%20(%5Cd%2B)%20people%20pass%22%2C%20_log)%0A%20%20%20%20_parts%20%3D%20%5B%5D%0A%20%20%20%20if%20_rm%3A%0A%20%20%20%20%20%20%20%20_parts.append(f%22Removed%3A%20**%7B_rm.group(1)%7D**%20of%20%7B_rm.group(2)%7D%20SNPs%20%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22(retention%3A%20%7B100*(1-int(_rm.group(1))%2Fint(_rm.group(2)))%3A.1f%7D%20%25)%22)%0A%20%20%20%20if%20_kp%3A%0A%20%20%20%20%20%20%20%20_parts.append(f%22Retained%3A%20**%7B_kp.group(1)%7D%20SNPs**%20across%20%7B_kp.group(2)%7D%20individuals%22)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd1%7D%5Cn%5Cn%7B_cmd2%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(%22%20%20%5Cn%22.join(_parts))%2C%20kind%3D%22success%22)%20if%20_parts%20else%20mo.md(%22%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60%5Cn%7B_log%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(prune_prefix%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%207%20%C2%B7%20MAF%20Distribution%20per%20Breed%0A%0A%20%20%20%20Following%20the%20reference%20study%2C%20SNPs%20are%20binned%20into%20four%20frequency%20classes%3A%0A%0A%20%20%20%20%7C%20Class%20%7C%20Range%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20Fixed%20%7C%20MAF%20%3D%200%20%7C%0A%20%20%20%20%7C%20Rare%20%7C%200%20%3C%20MAF%20%3C%200.05%20%7C%0A%20%20%20%20%7C%20Intermediate%20%7C%200.05%20%E2%89%A4%20MAF%20%3C%200.10%20%7C%0A%20%20%20%20%7C%20Common%20%7C%200.10%20%E2%89%A4%20MAF%20%E2%89%A4%200.50%20%7C%0A%0A%20%20%20%20Per-breed%20allele%20frequencies%20are%20computed%20with%20%60plink%20--freq%20--keep%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20mo%2C%20pd%2C%20plt%2C%20qc_prefix%2C%20subprocess)%3A%0A%20%20%20%20_cmd_frq%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--keep%20keep_%3CBreed%3E.txt%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--freq%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20%3CBreed%3E_freq%22)%0A%0A%20%20%20%20for%20_b%20in%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%3A%0A%20%20%20%20%20%20%20%20_frq_f%20%3D%20OUT%20%2F%20f%22%7B_b%7D_freq.frq%22%0A%20%20%20%20%20%20%20%20if%20not%20_frq_f.exists()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--keep%22%2Cstr(OUT%2Ff%22keep_%7B_b%7D.txt%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--freq%22%2C%22--chr-set%22%2C%2226%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--out%22%2Cstr(OUT%2Ff%22%7B_b%7D_freq%22)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20_bins%20%20%3D%20%5B(%22Fixed%22%2C%22%232A9D8F%22)%2C%20(%22Rare%22%2C%22%23457B9D%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20(%22Intermediate%22%2C%22%23E9C46A%22)%2C%20(%22Common%22%2C%22%23E63946%22)%5D%0A%20%20%20%20_fig_m%2C%20_axes_m%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(13%2C%204)%2C%20sharey%3DFalse)%0A%20%20%20%20_fig_m.suptitle(%22MAF%20Distribution%20per%20Breed%22%2C%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%0A%20%20%20%20_maf_rows%20%3D%20%5B%5D%0A%20%20%20%20for%20_ax%2C%20_b%20in%20zip(_axes_m%2C%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D)%3A%0A%20%20%20%20%20%20%20%20_frq%20%3D%20pd.read_csv(OUT%2Ff%22%7B_b%7D_freq.frq%22%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20%20%20%20%20_m%20%20%20%3D%20_frq%5B%22MAF%22%5D.values%0A%20%20%20%20%20%20%20%20_counts_maf%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Fixed%22%3A%20%20%20%20%20%20%20%20((_m%20%3D%3D%200)).sum()%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Rare%22%3A%20%20%20%20%20%20%20%20%20((_m%20%3E%200)%20%26%20(_m%20%3C%200.05)).sum()%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Intermediate%22%3A%20((_m%20%3E%3D%200.05)%20%26%20(_m%20%3C%200.10)).sum()%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Common%22%3A%20%20%20%20%20%20%20((_m%20%3E%3D%200.10)%20%26%20(_m%20%3C%3D%200.50)).sum()%2C%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20_labels%20%3D%20list(_counts_maf.keys())%0A%20%20%20%20%20%20%20%20_vals%20%20%20%3D%20list(_counts_maf.values())%0A%20%20%20%20%20%20%20%20_colors%20%3D%20%5Bc%20for%20_%2C%20c%20in%20_bins%5D%0A%20%20%20%20%20%20%20%20_ax.bar(_labels%2C%20_vals%2C%20color%3D_colors%2C%20edgecolor%3D%22white%22%2C%20width%3D0.6)%0A%20%20%20%20%20%20%20%20_ax.set_title(_b%2C%20fontsize%3D11%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22MAF%20class%22)%3B%20_ax.set_ylabel(%22Number%20of%20SNPs%22)%0A%20%20%20%20%20%20%20%20_ax.grid(axis%3D%22y%22%2C%20alpha%3D0.3%2C%20linestyle%3D%22--%22)%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20%20%20%20%20for%20_i%2C%20(_l%2C%20_v)%20in%20enumerate(zip(_labels%2C%20_vals))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.text(_i%2C%20_v%20%2B%2020%2C%20str(_v)%2C%20ha%3D%22center%22%2C%20fontsize%3D8)%0A%20%20%20%20%20%20%20%20_total%20%3D%20sum(_vals)%0A%20%20%20%20%20%20%20%20_maf_rows.append(%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Breed%22%3A%20_b%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20**%7Bk%3A%20f%22%7Bv%3A%2C%7D%20(%7B100*v%2F_total%3A.1f%7D%25)%22%20for%20k%2C%20v%20in%20_counts_maf.items()%7D%2C%0A%20%20%20%20%20%20%20%20%7D)%0A%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20_maf_df%20%3D%20pd.DataFrame(_maf_rows)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Command%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_frq%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_m)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_maf_df)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%208%20%C2%B7%20Principal%20Component%20Analysis%20(PCA)%0A%0A%20%20%20%20PCA%20is%20computed%20on%20the%20**LD-pruned**%20dataset%20(r%C2%B2%20%3D%200.3%2C%20%E2%89%88%2040%2C930%20SNPs)%2C%0A%20%20%20%20matching%20the%20reference%20study's%20approach%20of%20building%20a%20GRM%20from%20pruned%20SNPs.%0A%20%20%20%20PLINK%20computes%20the%20top%2010%20eigenvectors%3B%20we%20plot%20PC1%20vs%20PC2%20and%20PC1%20vs%20PC3.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20Path%2C%20mo%2C%20prune_prefix%2C%20subprocess)%3A%0A%20%20%20%20pca_prefix%20%20%20%3D%20str(OUT%20%2F%20%22sheep_pca_r03%22)%0A%20%20%20%20pca_eigenvec%20%3D%20Path(pca_prefix%20%2B%20%22.eigenvec%22)%0A%20%20%20%20_cmd%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_pruned_r03%20%5C%5C%5Cn%20%20%20%20--pca%2010%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20sheep_pca_r03%22)%0A%20%20%20%20if%20not%20pca_eigenvec.exists()%3A%0A%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cprune_prefix%2C%22--pca%22%2C%2210%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2Cpca_prefix%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Command%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22Output%3A%20%60%7Bpca_eigenvec.name%7D%60%20%C2%B7%20%60sheep_pca_r03.eigenval%60%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(pca_eigenvec%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20mpatches%2C%20np%2C%20pca_eigenvec%2C%20pd%2C%20plt%2C%20pop_df)%3A%0A%20%20%20%20_ev%20%3D%20pd.read_csv(pca_eigenvec%2C%20sep%3Dr%22%5Cs%2B%22%2C%20header%3DNone)%0A%20%20%20%20_ev.columns%20%3D%20%5B%22FID%22%2C%22IID%22%5D%20%2B%20%5Bf%22PC%7Bi%7D%22%20for%20i%20in%20range(1%2C%20_ev.shape%5B1%5D-1)%5D%0A%20%20%20%20_ev%20%3D%20_ev.merge(pop_df%5B%5B%22IID%22%2C%22POP%22%5D%5D%2C%20on%3D%22IID%22%2C%20how%3D%22left%22)%0A%20%20%20%20_evals%20%3D%20np.loadtxt(str(pca_eigenvec).replace(%22.eigenvec%22%2C%22.eigenval%22))%0A%20%20%20%20pct_var%20%20%20%20%20%3D%20_evals%20%2F%20_evals.sum()%20*%20100%0A%20%20%20%20POP_PALETTE%20%3D%20%7B%22Harri%22%3A%22%23E63946%22%2C%22Naemi%22%3A%22%23457B9D%22%2C%22Najdi%22%3A%22%232A9D8F%22%7D%0A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(13%2C5))%0A%20%20%20%20_fig.suptitle(%22PCA%20%E2%80%94%20Sheep%20Breed%20Differentiation%20(r%C2%B2%20%3D%200.3%20LD-pruned)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D13%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20for%20_ax%2C%20(_px%2C_py)%20in%20zip(_axes%2C%20%5B(1%2C2)%2C(1%2C3)%5D)%3A%0A%20%20%20%20%20%20%20%20for%20_pop%2C%20_grp%20in%20_ev.groupby(%22POP%22)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.scatter(_grp%5Bf%22PC%7B_px%7D%22%5D%2C%20_grp%5Bf%22PC%7B_py%7D%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%3DPOP_PALETTE%5B_pop%5D%2C%20label%3D_pop%2C%20s%3D65%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.85%2C%20edgecolors%3D%22white%22%2C%20linewidth%3D0.4)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(f%22PC%7B_px%7D%20(%7Bpct_var%5B_px-1%5D%3A.1f%7D%20%25)%22%2C%20fontsize%3D10)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(f%22PC%7B_py%7D%20(%7Bpct_var%5B_py-1%5D%3A.1f%7D%20%25)%22%2C%20fontsize%3D10)%0A%20%20%20%20%20%20%20%20_ax.set_title(f%22PC%7B_px%7D%20vs%20PC%7B_py%7D%22)%0A%20%20%20%20%20%20%20%20_ax.grid(alpha%3D0.25%2C%20linestyle%3D%22--%22)%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20_handles%20%3D%20%5Bmpatches.Patch(color%3Dc%2Clabel%3Dp)%20for%20p%2Cc%20in%20POP_PALETTE.items()%5D%0A%20%20%20%20_fig.legend(handles%3D_handles%2C%20loc%3D%22lower%20center%22%2C%20ncol%3D3%2C%20fontsize%3D10%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bbox_to_anchor%3D(0.5%2C-0.04)%2C%20frameon%3DFalse)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20pop_palette%20%3D%20POP_PALETTE%0A%20%20%20%20mo.vstack(%5Bmo.md(%22%23%23%23%20PC1%20vs%20PC2%20%C2%B7%20PC1%20vs%20PC3%22)%2C%20mo.as_html(_fig)%5D)%0A%20%20%20%20return%20pct_var%2C%20pop_palette%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20np%2C%20pct_var%2C%20pd%2C%20plt)%3A%0A%20%20%20%20_pcs%20%20%3D%20np.arange(1%2C%20len(pct_var)%2B1)%0A%20%20%20%20_fig2%2C%20_ax2%20%3D%20plt.subplots(figsize%3D(7%2C4))%0A%20%20%20%20_ax2.bar(_pcs%2C%20pct_var%2C%20color%3D%22%23457B9D%22%2C%20edgecolor%3D%22white%22%2C%20width%3D0.7)%0A%20%20%20%20_ax2.set_xlabel(%22Principal%20Component%22)%3B%20_ax2.set_ylabel(%22%25%20Variance%20Explained%22)%0A%20%20%20%20_ax2.set_title(%22Scree%20Plot%22%2C%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax2.set_xticks(_pcs)%0A%20%20%20%20_ax2.grid(axis%3D%22y%22%2C%20alpha%3D0.35%2C%20linestyle%3D%22--%22)%0A%20%20%20%20_ax2.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20_scree%20%3D%20pd.DataFrame(%7B%22PC%22%3A_pcs%2C%22%25%20Variance%22%3Apct_var.round(2)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Cumulative%20%25%22%3Apct_var.cumsum().round(2)%7D)%0A%20%20%20%20mo.vstack(%5Bmo.md(%22%23%23%23%20Scree%20plot%22)%2C%20mo.as_html(_fig2)%2C%20mo.ui.table(_scree)%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%209%20%C2%B7%20Population%20Structure%20%E2%80%94%20NMF%20K%20%3D%201%E2%80%936%0A%0A%20%20%20%20The%20reference%20study%20used%20**sNMF**%20(LEA%20R%20package)%20testing%20K%20%3D%201%E2%80%9310%20and%0A%20%20%20%20selected%20the%20optimal%20K%20via%20minimum%20cross-entropy%20across%20replicate%20runs.%0A%0A%20%20%20%20Here%20we%20use%20**NMF**%20(scikit-learn)%20on%20the%20LD-pruned%20(r%C2%B2%20%3D%200.3)%20genotype%0A%20%20%20%20dosage%20matrix.%20%20We%20test%20K%20%3D%201%E2%80%936%20and%20use%20NMF%20**reconstruction%20error**%0A%20%20%20%20(analogous%20to%20cross-entropy)%20to%20identify%20the%20optimal%20K%20at%20the%20%22elbow%22.%0A%0A%20%20%20%20**Steps%3A**%0A%20%20%20%201.%20%60plink%20--recodeA%60%20exports%20LD-pruned%20genotypes%20as%200%2F1%2F2%20dosages%0A%20%20%20%202.%20Missing%20values%20imputed%20with%20SNP-mean%0A%20%20%20%203.%20NMF%20fit%20for%20K%20%3D%201%E2%80%936%3B%20reconstruction%20error%20plotted%0A%20%20%20%204.%20Admixture%20bar%20charts%20for%20all%20K%20stacked%20vertically%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20mo%2C%20prune_prefix%2C%20subprocess)%3A%0A%20%20%20%20raw_file%20%3D%20OUT%20%2F%20%22sheep_geno_r03.raw%22%0A%20%20%20%20_cmd%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_pruned_r03%20%5C%5C%5Cn%20%20%20%20--recodeA%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20sheep_geno_r03%22)%0A%20%20%20%20if%20not%20raw_file.exists()%3A%0A%20%20%20%20%20%20%20%20_r%20%3D%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cprune_prefix%2C%22--recodeA%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2Cstr(OUT%2F%22sheep_geno_r03%22)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20%20%20%20%20_st%20%3D%20%22%E2%9C%93%22%20if%20_r.returncode%3D%3D0%20else%20f%22%E2%9C%97%20%7B_r.stderr%5B%3A150%5D%7D%22%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20_st%20%3D%20%22skipped%20(file%20exists)%22%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Command%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22Status%3A%20%7B_st%7D%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(raw_file%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20np%2C%20pd%2C%20plt%2C%20pop_df%2C%20raw_file)%3A%0A%20%20%20%20from%20sklearn.decomposition%20import%20NMF%20as%20_NMF%0A%20%20%20%20from%20sklearn.preprocessing%20import%20normalize%20as%20_norm%0A%0A%20%20%20%20_raw%20%3D%20pd.read_csv(raw_file%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20_ids%20%3D%20_raw%5B%22IID%22%5D.values%0A%20%20%20%20_G%20%20%20%3D%20_raw.iloc%5B%3A%2C6%3A%5D.values.astype(float)%0A%20%20%20%20_cm%20%20%3D%20np.nanmean(_G%2C%20axis%3D0)%0A%20%20%20%20_nm%20%20%3D%20np.isnan(_G)%0A%20%20%20%20_G%5B_nm%5D%20%3D%20np.take(_cm%2C%20np.where(_nm)%5B1%5D)%0A%0A%20%20%20%20_K_range%20%20%20%3D%20list(range(1%2C%207))%0A%20%20%20%20_errors%20%20%20%20%3D%20%5B%5D%0A%20%20%20%20_Q_all%20%20%20%20%20%3D%20%7B%7D%0A%20%20%20%20_pop_order%20%3D%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%0A%20%20%20%20_palettes%20%20%3D%20%7B%0A%20%20%20%20%20%20%20%201%3A%5B%22%23E63946%22%5D%2C%0A%20%20%20%20%20%20%20%202%3A%5B%22%23E63946%22%2C%22%23457B9D%22%5D%2C%0A%20%20%20%20%20%20%20%203%3A%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%5D%2C%0A%20%20%20%20%20%20%20%204%3A%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%5D%2C%0A%20%20%20%20%20%20%20%205%3A%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%2C%22%238338EC%22%5D%2C%0A%20%20%20%20%20%20%20%206%3A%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%2C%22%238338EC%22%2C%22%23FB8500%22%5D%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20for%20_K%20in%20_K_range%3A%0A%20%20%20%20%20%20%20%20_m%20%3D%20_NMF(n_components%3D_K%2C%20init%3D%22nndsvda%22%2C%20max_iter%3D600%2C%20random_state%3D42)%0A%20%20%20%20%20%20%20%20_W%20%3D%20_m.fit_transform(_G)%0A%20%20%20%20%20%20%20%20_errors.append(_m.reconstruction_err_)%0A%20%20%20%20%20%20%20%20_rs%20%3D%20_W.sum(axis%3D1%2C%20keepdims%3DTrue)%3B%20_rs%5B_rs%3D%3D0%5D%3D1%0A%20%20%20%20%20%20%20%20_Q%20%20%3D%20_W%20%2F%20_rs%0A%20%20%20%20%20%20%20%20_qdf%20%3D%20pd.DataFrame(_Q%2C%20columns%3D%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D)%0A%20%20%20%20%20%20%20%20_qdf%5B%22IID%22%5D%20%3D%20_ids%0A%20%20%20%20%20%20%20%20_qdf%20%3D%20_qdf.merge(pop_df%5B%5B%22IID%22%2C%22POP%22%5D%5D%2C%20on%3D%22IID%22%2C%20how%3D%22left%22)%0A%20%20%20%20%20%20%20%20_qdf%5B%22dom%22%5D%20%3D%20_qdf%5B%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D%5D.idxmax(axis%3D1)%0A%20%20%20%20%20%20%20%20_Q_all%5B_K%5D%20%3D%20pd.concat(%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B_qdf%5B_qdf%5B%22POP%22%5D%3D%3Dp%5D.sort_values(%22dom%22)%20for%20p%20in%20_pop_order%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20ignore_index%3DTrue)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20cross-entropy%20(reconstruction%20error)%20plot%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig_ce%2C%20_ax_ce%20%3D%20plt.subplots(figsize%3D(6%2C3))%0A%20%20%20%20_ax_ce.plot(_K_range%2C%20_errors%2C%20%22o-%22%2C%20color%3D%22%23E63946%22%2C%20linewidth%3D2%2C%20markersize%3D7)%0A%20%20%20%20_ax_ce.set_xlabel(%22K%20(number%20of%20ancestral%20components)%22)%0A%20%20%20%20_ax_ce.set_ylabel(%22Reconstruction%20error%22)%0A%20%20%20%20_ax_ce.set_title(%22NMF%20Reconstruction%20Error%20vs%20K%22%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax_ce.set_xticks(_K_range)%0A%20%20%20%20_ax_ce.grid(alpha%3D0.3%2C%20linestyle%3D%22--%22)%0A%20%20%20%20_ax_ce.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20admixture%20bar%20charts%20K%3D2%E2%80%936%20(stacked%20vertically)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_show_K%20%3D%20%5B2%2C3%2C4%2C5%2C6%5D%0A%20%20%20%20_fig_s%2C%20_axes_s%20%3D%20plt.subplots(len(_show_K)%2C%201%2C%20figsize%3D(14%2C%203*len(_show_K)))%0A%20%20%20%20_fig_s.suptitle(%22Population%20Structure%20%E2%80%94%20NMF%20Admixture%20(K%20%3D%202%E2%80%936)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D13%2C%20fontweight%3D%22bold%22%2C%20y%3D1.01)%0A%0A%20%20%20%20for%20_ax%2C%20_K%20in%20zip(_axes_s%2C%20_show_K)%3A%0A%20%20%20%20%20%20%20%20_sorted%20%3D%20_Q_all%5B_K%5D%0A%20%20%20%20%20%20%20%20_btm%20%3D%20np.zeros(len(_sorted))%0A%20%20%20%20%20%20%20%20for%20_ki%2C%20_col%20in%20enumerate(%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.bar(range(len(_sorted))%2C%20_sorted%5B_col%5D.values%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bottom%3D_btm%2C%20color%3D_palettes%5B_K%5D%5B_ki%5D%2C%20width%3D1.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Anc.%20%7B_ki%2B1%7D%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20_btm%20%2B%3D%20_sorted%5B_col%5D.values%0A%20%20%20%20%20%20%20%20_sz%20%3D%20%5Blen(_sorted%5B_sorted%5B%22POP%22%5D%3D%3Dp%5D)%20for%20p%20in%20_pop_order%5D%0A%20%20%20%20%20%20%20%20_bd%20%3D%20np.cumsum(%5B0%5D%2B_sz)%0A%20%20%20%20%20%20%20%20for%20_b%20in%20_bd%5B1%3A-1%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.axvline(_b-0.5%2C%20color%3D%22black%22%2C%20linewidth%3D1.5)%0A%20%20%20%20%20%20%20%20_md%20%3D%20(_bd%5B%3A-1%5D%2B_bd%5B1%3A%5D)%2F2%0A%20%20%20%20%20%20%20%20for%20_pos%2C_pop%20in%20zip(_md%2C_pop_order)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.text(_pos%2C1.04%2C_pop%2Cha%3D%22center%22%2Cva%3D%22bottom%22%2Cfontsize%3D10%2Cfontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlim(-0.5%2Clen(_sorted)-0.5)%3B%20_ax.set_ylim(0%2C1.18)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(%22Ancestry%5Cnproportion%22%2C%20fontsize%3D9)%3B%20_ax.set_xticks(%5B%5D)%0A%20%20%20%20%20%20%20%20_ax.set_title(f%22K%20%3D%20%7B_K%7D%22%2C%20fontsize%3D11%2C%20fontweight%3D%22bold%22%2C%20loc%3D%22left%22)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22upper%20right%22%2Cfontsize%3D8%2Cframeon%3DFalse%2Cncol%3D_K%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bbox_to_anchor%3D(1%2C1.2))%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%22right%22%2C%22bottom%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20_err_df%20%3D%20pd.DataFrame(%7B%22K%22%3A_K_range%2C%20%22Reconstruction%20Error%22%3A%5Bround(e%2C1)%20for%20e%20in%20_errors%5D%7D)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Cross-entropy%20proxy%20(reconstruction%20error)%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22The%20optimal%20K%20is%20at%20the%20**elbow**%20%E2%80%94%20where%20error%20stops%20decreasing%20steeply.%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_ce)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_err_df)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Admixture%20bar%20charts%20(K%20%3D%202%E2%80%936)%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22Each%20bar%20%3D%20one%20individual.%20%20Breeds%20are%20separated%20by%20vertical%20lines.%20%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Individuals%20sorted%20by%20dominant%20component%20within%20each%20breed.%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_s)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%20%C2%A79b%20%C2%B7%20Population%20Structure%20%E2%80%94%20ADMIXTURE%20(Alternative)%0A%0A%20%20%20%20K%20%3D%201%E2%80%936%20was%20tested%20with%20**3%20independent%20replicates**%20(seeds%2042%2C%20123%2C%20456).%0A%20%20%20%20The%20replicate%20with%20the%20lowest%20CV%20error%20is%20shown%20for%20each%20K.%0A%0A%20%20%20%20%3E%20**Result%3A%20K%20%3D%203%20is%20optimal**%20%E2%80%94%20consistent%20with%20the%20three%20known%20breeds%0A%20%20%20%20%3E%20(Harri%2C%20Naemi%2C%20Najdi).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20mo%2C%20np%2C%20pd%2C%20plt%2C%20pop_df)%3A%0A%20%20%20%20import%20re%20as%20_re2%0A%0A%20%20%20%20_cmd_admix%20%3D%20(%0A%20%20%20%20%20%20%20%20%22%23%20Run%20for%20K%3D1%E2%80%936%20with%203%20replicates%20(seeds%2042%2C%20123%2C%20456)%5Cn%22%0A%20%20%20%20%20%20%20%20%22for%20K%20in%201%202%203%204%205%206%3B%20do%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20for%20SEED%20in%2042%20123%20456%3B%20do%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20%20%20%20%20admixture%20--cv%3D5%20--seed%3D%24SEED%20-j4%20sheep_pruned_r03.bed%20%24K%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20%20%20%20%20%20%20%20%20%3E%20admixture_K%24%7BK%7D_S%24%7BSEED%7D.log%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20%20%20%20%20mv%20sheep_pruned_r03.%24%7BK%7D.Q%20admixture_K%24%7BK%7D_S%24%7BSEED%7D.Q%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20%20%20%20%20mv%20sheep_pruned_r03.%24%7BK%7D.P%20admixture_K%24%7BK%7D_S%24%7BSEED%7D.P%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20done%5Cn%22%0A%20%20%20%20%20%20%20%20%22done%22%0A%20%20%20%20)%0A%0A%20%20%20%20_seeds_a%20%20%20%3D%20%5B42%2C%20123%2C%20456%5D%0A%20%20%20%20_K_range_a%20%3D%20list(range(1%2C%207))%0A%20%20%20%20_pal_a%20%3D%20%7B%0A%20%20%20%20%20%20%20%201%3A%20%5B%22%23E63946%22%5D%2C%0A%20%20%20%20%20%20%20%202%3A%20%5B%22%23E63946%22%2C%22%23457B9D%22%5D%2C%0A%20%20%20%20%20%20%20%203%3A%20%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%5D%2C%0A%20%20%20%20%20%20%20%204%3A%20%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%5D%2C%0A%20%20%20%20%20%20%20%205%3A%20%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%2C%22%238338EC%22%5D%2C%0A%20%20%20%20%20%20%20%206%3A%20%5B%22%23E63946%22%2C%22%23457B9D%22%2C%22%232A9D8F%22%2C%22%23F4A261%22%2C%22%238338EC%22%2C%22%23FB8500%22%5D%2C%0A%20%20%20%20%7D%0A%20%20%20%20_pop_order_a%20%3D%20%5B%22Harri%22%2C%20%22Naemi%22%2C%20%22Najdi%22%5D%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20read%20CV%20errors%20from%20logs%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_cv_data%20%3D%20%7BK%3A%20%5B%5D%20for%20K%20in%20_K_range_a%7D%0A%20%20%20%20_best_a%20%20%3D%20%7B%7D%0A%20%20%20%20for%20_K%20in%20_K_range_a%3A%0A%20%20%20%20%20%20%20%20for%20_s%20in%20_seeds_a%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_log%20%3D%20(OUT%20%2F%20f%22admixture_K%7B_K%7D_S%7B_s%7D.log%22).read_text()%0A%20%20%20%20%20%20%20%20%20%20%20%20_m%20%20%20%3D%20_re2.search(r%22CV%20error%20%5C(K%3D%5Cd%2B%5C)%3A%5Cs%2B(%5B%5Cd.%5D%2B)%22%2C%20_log)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_m%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_cv_data%5B_K%5D.append((_s%2C%20float(_m.group(1))))%0A%20%20%20%20%20%20%20%20_best_a%5B_K%5D%20%3D%20min(_cv_data%5B_K%5D%2C%20key%3Dlambda%20x%3A%20x%5B1%5D)%0A%0A%20%20%20%20_cv_means_a%20%3D%20%5Bnp.mean(%5Bv%20for%20_%2C%20v%20in%20_cv_data%5BK%5D%5D)%20for%20K%20in%20_K_range_a%5D%0A%20%20%20%20_cv_sds_a%20%20%20%3D%20%5Bnp.std(%20%5Bv%20for%20_%2C%20v%20in%20_cv_data%5BK%5D%5D)%20for%20K%20in%20_K_range_a%5D%0A%20%20%20%20_best_K_a%20%20%20%3D%20_K_range_a%5Bint(np.argmin(_cv_means_a))%5D%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20CV%20error%20plot%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig_cv%2C%20_ax_cv%20%3D%20plt.subplots(figsize%3D(6%2C%203))%0A%20%20%20%20_ax_cv.errorbar(_K_range_a%2C%20_cv_means_a%2C%20yerr%3D_cv_sds_a%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fmt%3D%22o-%22%2C%20color%3D%22%23457B9D%22%2C%20linewidth%3D2%2C%20markersize%3D7%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capsize%3D4%2C%20label%3D%22mean%20%C2%B1%20SD%20(3%20replicates)%22)%0A%20%20%20%20_ax_cv.scatter(%5B_best_K_a%5D%2C%20%5B_cv_means_a%5B_best_K_a%20-%201%5D%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20color%3D%22%23E63946%22%2C%20s%3D120%2C%20zorder%3D5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Optimal%20K%20%3D%20%7B_best_K_a%7D%22)%0A%20%20%20%20_ax_cv.set_xlabel(%22K%20(number%20of%20ancestral%20components)%22)%0A%20%20%20%20_ax_cv.set_ylabel(%22Cross-validation%20error%22)%0A%20%20%20%20_ax_cv.set_title(%22ADMIXTURE%20Cross-Validation%20Error%22%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax_cv.set_xticks(_K_range_a)%0A%20%20%20%20_ax_cv.legend(fontsize%3D9%2C%20frameon%3DFalse)%0A%20%20%20%20_ax_cv.grid(alpha%3D0.3%2C%20linestyle%3D%22--%22)%0A%20%20%20%20_ax_cv.spines%5B%5B%22top%22%2C%20%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20read%20FAM%20for%20sample%20order%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fam%20%3D%20pd.read_csv(OUT%20%2F%20%22sheep_pruned_r03.fam%22%2C%20sep%3Dr%22%5Cs%2B%22%2C%20header%3DNone%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20usecols%3D%5B1%5D%2C%20names%3D%5B%22IID%22%5D)%0A%20%20%20%20_fam%20%3D%20_fam.merge(pop_df%5B%5B%22IID%22%2C%20%22POP%22%5D%5D%2C%20on%3D%22IID%22%2C%20how%3D%22left%22)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20admixture%20bar%20charts%20K%3D2%E2%80%936%20(best%20replicate%20per%20K)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_show_Ka%20%3D%20%5B2%2C%203%2C%204%2C%205%2C%206%5D%0A%20%20%20%20_fig_a%2C%20_axes_a%20%3D%20plt.subplots(len(_show_Ka)%2C%201%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20figsize%3D(14%2C%203%20*%20len(_show_Ka)))%0A%20%20%20%20_fig_a.suptitle(%22Population%20Structure%20%E2%80%94%20ADMIXTURE%20(K%20%3D%202%E2%80%936%2C%20best%20replicate)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D13%2C%20fontweight%3D%22bold%22%2C%20y%3D1.01)%0A%0A%20%20%20%20for%20_ax%2C%20_K%20in%20zip(_axes_a%2C%20_show_Ka)%3A%0A%20%20%20%20%20%20%20%20_bs%20%3D%20_best_a%5B_K%5D%5B0%5D%0A%20%20%20%20%20%20%20%20_Q%20%20%3D%20np.loadtxt(OUT%20%2F%20f%22admixture_K%7B_K%7D_S%7B_bs%7D.Q%22)%0A%20%20%20%20%20%20%20%20_qdf%20%3D%20pd.DataFrame(_Q%2C%20columns%3D%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D)%0A%20%20%20%20%20%20%20%20_qdf%5B%22IID%22%5D%20%3D%20_fam%5B%22IID%22%5D.values%0A%20%20%20%20%20%20%20%20_qdf%5B%22POP%22%5D%20%3D%20_fam%5B%22POP%22%5D.values%0A%20%20%20%20%20%20%20%20_qdf%5B%22dom%22%5D%20%3D%20_qdf%5B%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D%5D.idxmax(axis%3D1)%0A%20%20%20%20%20%20%20%20_sorted%20%3D%20pd.concat(%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B_qdf%5B_qdf%5B%22POP%22%5D%20%3D%3D%20p%5D.sort_values(%22dom%22)%20for%20p%20in%20_pop_order_a%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20ignore_index%3DTrue)%0A%20%20%20%20%20%20%20%20_btm%20%3D%20np.zeros(len(_sorted))%0A%20%20%20%20%20%20%20%20for%20_ki%2C%20_col%20in%20enumerate(%5Bf%22K%7Bk%2B1%7D%22%20for%20k%20in%20range(_K)%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.bar(range(len(_sorted))%2C%20_sorted%5B_col%5D.values%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bottom%3D_btm%2C%20color%3D_pal_a%5B_K%5D%5B_ki%5D%2C%20width%3D1.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Anc.%20%7B_ki%2B1%7D%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20_btm%20%2B%3D%20_sorted%5B_col%5D.values%0A%20%20%20%20%20%20%20%20_sz%20%3D%20%5Blen(_sorted%5B_sorted%5B%22POP%22%5D%20%3D%3D%20p%5D)%20for%20p%20in%20_pop_order_a%5D%0A%20%20%20%20%20%20%20%20_bd%20%3D%20np.cumsum(%5B0%5D%20%2B%20_sz)%0A%20%20%20%20%20%20%20%20for%20_b%20in%20_bd%5B1%3A-1%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.axvline(_b%20-%200.5%2C%20color%3D%22black%22%2C%20linewidth%3D1.5)%0A%20%20%20%20%20%20%20%20_md%20%3D%20(_bd%5B%3A-1%5D%20%2B%20_bd%5B1%3A%5D)%20%2F%202%0A%20%20%20%20%20%20%20%20for%20_pos%2C%20_pop%20in%20zip(_md%2C%20_pop_order_a)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.text(_pos%2C%201.04%2C%20_pop%2C%20ha%3D%22center%22%2C%20va%3D%22bottom%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D10%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlim(-0.5%2C%20len(_sorted)%20-%200.5)%0A%20%20%20%20%20%20%20%20_ax.set_ylim(0%2C%201.18)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(%22Ancestry%5Cnproportion%22%2C%20fontsize%3D9)%0A%20%20%20%20%20%20%20%20_ax.set_xticks(%5B%5D)%0A%20%20%20%20%20%20%20%20_ax.set_title(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22K%20%3D%20%7B_K%7D%20%20(seed%20%7B_bs%7D%2C%20CV%20%3D%20%7B_best_a%5B_K%5D%5B1%5D%3A.5f%7D)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D11%2C%20fontweight%3D%22bold%22%2C%20loc%3D%22left%22)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22upper%20right%22%2C%20fontsize%3D8%2C%20frameon%3DFalse%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ncol%3D_K%2C%20bbox_to_anchor%3D(1%2C%201.2))%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%20%22right%22%2C%20%22bottom%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20CV%20summary%20table%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_cv_df%20%3D%20pd.DataFrame(%7B%0A%20%20%20%20%20%20%20%20%22K%22%3A%20%20%20%20%20%20%20%20%20%20_K_range_a%2C%0A%20%20%20%20%20%20%20%20%22CV%20mean%22%3A%20%20%20%20%5Bround(v%2C%205)%20for%20v%20in%20_cv_means_a%5D%2C%0A%20%20%20%20%20%20%20%20%22CV%20SD%22%3A%20%20%20%20%20%20%5Bround(v%2C%205)%20for%20v%20in%20_cv_sds_a%5D%2C%0A%20%20%20%20%20%20%20%20%22Best%20seed%22%3A%20%20%5B_best_a%5BK%5D%5B0%5D%20for%20K%20in%20_K_range_a%5D%2C%0A%20%20%20%20%20%20%20%20%22Best%20CV%22%3A%20%20%20%20%5Bround(_best_a%5BK%5D%5B1%5D%2C%205)%20for%20K%20in%20_K_range_a%5D%2C%0A%20%20%20%20%7D)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Command%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_admix%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22**Optimal%20K%20%3D%20%7B_best_K_a%7D**%20%E2%80%94%20lowest%20mean%20CV%20error%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22(%7B_cv_means_a%5B_best_K_a-1%5D%3A.5f%7D).%20%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Consistent%20with%20the%203%20known%20breeds%20(Harri%20%C2%B7%20Naemi%20%C2%B7%20Najdi).%22%0A%20%20%20%20%20%20%20%20)%2C%20kind%3D%22success%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Cross-validation%20error%20(mean%20%C2%B1%20SD%2C%203%20replicates)%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_cv)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_cv_df)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Admixture%20bar%20charts%20(K%20%3D%202%E2%80%936%2C%20best%20replicate%20per%20K)%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22Each%20bar%20%3D%20one%20individual.%20%20Breeds%20separated%20by%20vertical%20lines.%20%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Individuals%20sorted%20by%20dominant%20ancestry%20component%20within%20each%20breed.%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_a)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%2010%20%C2%B7%20Genetic%20Diversity%20%E2%80%94%20Ho%2C%20He%2C%20F_IS%20and%20Ne%20Trajectory%0A%0A%20%20%20%20%7C%20Metric%20%7C%20Source%20%7C%20Formula%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20**Ho**%20%7C%20%60plink%20--het%60%20%7C%20(N_NM%20%E2%88%92%20O_HOM)%20%2F%20N_NM%20%7C%0A%20%20%20%20%7C%20**He**%20%7C%20%60plink%20--freq%20--keep%60%20%7C%20mean(2pq)%20across%20SNPs%20%7C%0A%20%20%20%20%7C%20**F_IS**%20%7C%20derived%20%7C%20(He%20%E2%88%92%20Ho)%20%2F%20He%20%7C%0A%20%20%20%20%7C%20**Ne(t)**%20%7C%20%60plink%20--r2%20--ld-window-kb%201000%60%20%7C%20Sved%20(1971)%3A%20Ne%20%3D%20(1%2Fr%C2%B2%20%E2%88%92%201)%20%2F%204c%20%7C%0A%0A%20%20%20%20The%20**Ne%20trajectory**%20shows%20historical%20effective%20population%20size%20at%0A%20%20%20%20multiple%20generations%20in%20the%20past%2C%20estimated%20from%20pairwise%20LD%20at%20different%0A%20%20%20%20chromosomal%20distances%20using%20a%201%20Mb%20window%20(matching%20the%20reference%20study).%0A%0A%20%20%20%20Distance%20class%20%E2%86%92%20generation%3A%20t%20%E2%89%88%201%20%2F%20(2c)%20where%20c%20%3D%20distance%20in%20Morgans.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20OUT%2C%0A%20%20%20%20PLINK%2C%0A%20%20%20%20het_file%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20pd%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20pop_df%2C%0A%20%20%20%20pop_palette%2C%0A%20%20%20%20qc_prefix%2C%0A%20%20%20%20sns%2C%0A%20%20%20%20subprocess%2C%0A)%3A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20commands%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_cmd_het2%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--het%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20sheep_het%22)%0A%20%20%20%20_cmd_frq2%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--keep%20keep_%3CBreed%3E.txt%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--freq%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20%3CBreed%3E_freq%22)%0A%20%20%20%20_cmd_ld2%20%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--keep%20keep_%3CBreed%3E.txt%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--r2%20%5C%5C%5Cn%20%20%20%20--ld-window%20999999%20%5C%5C%5Cn%20%20%20%20--ld-window-kb%201000%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--ld-window-r2%200%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%20%20%20%20--out%20%3CBreed%3E_ld1mb%22)%0A%0A%20%20%20%20for%20_b%20in%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%3A%0A%20%20%20%20%20%20%20%20_frq_p%20%3D%20OUT%20%2F%20f%22%7B_b%7D_freq.frq%22%0A%20%20%20%20%20%20%20%20if%20not%20_frq_p.exists()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%22--keep%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20str(OUT%2Ff%22keep_%7B_b%7D.txt%22)%2C%22--freq%22%2C%22--chr-set%22%2C%2226%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--out%22%2Cstr(OUT%2Ff%22%7B_b%7D_freq%22)%5D%2Ccapture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20%20%20%20%20_ld_p%20%3D%20OUT%20%2F%20f%22%7B_b%7D_ld1mb.ld%22%0A%20%20%20%20%20%20%20%20if%20not%20_ld_p.exists()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%22--keep%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20str(OUT%2Ff%22keep_%7B_b%7D.txt%22)%2C%22--r2%22%2C%22--ld-window%22%2C%22999999%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--ld-window-kb%22%2C%221000%22%2C%22--ld-window-r2%22%2C%220%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2Cstr(OUT%2Ff%22%7B_b%7D_ld1mb%22)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Ho%20from%20--het%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_het%20%3D%20pd.read_csv(het_file%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20_het.columns%20%3D%20%5Bc.strip()%20for%20c%20in%20_het.columns%5D%0A%20%20%20%20_het%20%3D%20_het.merge(pop_df%5B%5B%22IID%22%2C%22POP%22%5D%5D%2C%20on%3D%22IID%22%2C%20how%3D%22left%22)%0A%20%20%20%20_het%5B%22Ho%22%5D%20%3D%20(_het%5B%22N(NM)%22%5D%20-%20_het%5B%22O(HOM)%22%5D)%20%2F%20_het%5B%22N(NM)%22%5D%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20He%2C%20F_IS%2C%20single-point%20Ne%20(50%20kb)%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_div_rows%20%3D%20%5B%5D%0A%20%20%20%20for%20_b%20in%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%3A%0A%20%20%20%20%20%20%20%20_ho%20%20%3D%20_het%5B_het%5B%22POP%22%5D%3D%3D_b%5D%5B%22Ho%22%5D.mean()%0A%20%20%20%20%20%20%20%20_frq%20%3D%20pd.read_csv(OUT%2Ff%22%7B_b%7D_freq.frq%22%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20%20%20%20%20_he%20%20%3D%20(2*_frq%5B%22MAF%22%5D*(1-_frq%5B%22MAF%22%5D)).mean()%0A%20%20%20%20%20%20%20%20_fis%20%3D%20(_he%20-%20_ho)%20%2F%20_he%20if%20_he%20%3E%200%20else%20float(%22nan%22)%0A%20%20%20%20%20%20%20%20_div_rows.append(%7B%22Breed%22%3A_b%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Ho%22%3Around(_ho%2C4)%2C%20%22He%22%3Around(_he%2C4)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22F_IS%22%3Around(_fis%2C4)%7D)%0A%20%20%20%20_div_df%20%3D%20pd.DataFrame(_div_rows)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20violin%20plot%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig_h%2C%20(_a1%2C_a2)%20%3D%20plt.subplots(1%2C2%2Cfigsize%3D(13%2C5))%0A%20%20%20%20_fig_h.suptitle(%22Observed%20Heterozygosity%20(Ho)%22%2C%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ord%20%3D%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%0A%20%20%20%20sns.violinplot(data%3D_het%2Cx%3D%22POP%22%2Cy%3D%22Ho%22%2Cpalette%3Dpop_palette%2Cax%3D_a1%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20inner%3D%22box%22%2Corder%3D_ord)%0A%20%20%20%20_a1.set_xlabel(%22Breed%22)%3B%20_a1.set_ylabel(%22Ho%22)%0A%20%20%20%20_a1.grid(axis%3D%22y%22%2Calpha%3D0.3%2Clinestyle%3D%22--%22)%3B%20_a1.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20sns.stripplot(data%3D_het%2Cx%3D%22POP%22%2Cy%3D%22Ho%22%2Cpalette%3Dpop_palette%2Cjitter%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.55%2Cax%3D_a2%2Csize%3D6%2Corder%3D_ord)%0A%20%20%20%20_a2.set_xlabel(%22Breed%22)%3B%20_a2.set_ylabel(%22Ho%22)%0A%20%20%20%20_a2.grid(axis%3D%22y%22%2Calpha%3D0.3%2Clinestyle%3D%22--%22)%3B%20_a2.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Ne%20trajectory%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20%23%20distance%20bins%20(bp)%20and%20midpoints%20in%20Morgans%20(1%20cM%20per%20Mb)%0A%20%20%20%20_dist_bins%20%3D%20%5B%0A%20%20%20%20%20%20%20%20(%20%20%20%20%200%2C%20%2025_000)%2C%20%20(%2025_000%2C%20%2050_000)%2C%20(%2050_000%2C%20100_000)%2C%0A%20%20%20%20%20%20%20%20(100_000%2C%20200_000)%2C%20%20(200_000%2C%20400_000)%2C%20(400_000%2C%20700_000)%2C%0A%20%20%20%20%20%20%20%20(700_000%2C1_000_000)%2C%0A%20%20%20%20%5D%0A%20%20%20%20_fig_ne%2C%20_ax_ne%20%3D%20plt.subplots(figsize%3D(9%2C5))%0A%20%20%20%20_ne_all%20%3D%20%7B%7D%0A%20%20%20%20for%20_b%20in%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%3A%0A%20%20%20%20%20%20%20%20_ld%20%3D%20pd.read_csv(OUT%2Ff%22%7B_b%7D_ld1mb.ld%22%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20%20%20%20%20_ld%5B%22dist%22%5D%20%3D%20(_ld%5B%22BP_B%22%5D-_ld%5B%22BP_A%22%5D).abs()%0A%20%20%20%20%20%20%20%20_ne_pts%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20_lo%2C_hi%20in%20_dist_bins%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_sub%20%3D%20_ld%5B(_ld%5B%22dist%22%5D%3E%3D_lo)%20%26%20(_ld%5B%22dist%22%5D%3C_hi)%5D%5B%22R2%22%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20len(_sub)%20%3C%2010%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20%20%20%20%20_r2m%20%3D%20_sub.mean()%0A%20%20%20%20%20%20%20%20%20%20%20%20_c%20%20%20%3D%20(_lo%2B_hi)%2F(2*1e8)%20%20%20%20%20%20%20%20%20%20%23%20midpoint%20in%20Morgans%0A%20%20%20%20%20%20%20%20%20%20%20%20_t%20%20%20%3D%20int(round(1%2F(2*_c)))%20%20%20%20%20%20%20%20%23%20generations%20ago%0A%20%20%20%20%20%20%20%20%20%20%20%20_ne%20%20%3D%20(1%2F_r2m%20-%201)%2F(4*_c)%20if%20_r2m%3E0%20else%20float(%22nan%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ne_pts.append((_t%2C%20_ne))%0A%20%20%20%20%20%20%20%20_ne_pts.sort(key%3Dlambda%20x%3Ax%5B0%5D)%0A%20%20%20%20%20%20%20%20_ts%20%20%3D%20%5Bp%5B0%5D%20for%20p%20in%20_ne_pts%5D%0A%20%20%20%20%20%20%20%20_nes%20%3D%20%5Bp%5B1%5D%20for%20p%20in%20_ne_pts%5D%0A%20%20%20%20%20%20%20%20_ne_all%5B_b%5D%20%3D%20(_ts%2C%20_nes)%0A%20%20%20%20%20%20%20%20_ax_ne.plot(_ts%2C%20_nes%2C%20%22o-%22%2C%20color%3Dpop_palette%5B_b%5D%2C%20label%3D_b%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20linewidth%3D2%2C%20markersize%3D6)%0A%0A%20%20%20%20_ax_ne.set_xlabel(%22Generations%20ago%20(t%20%3D%201%2F2c)%22)%0A%20%20%20%20_ax_ne.set_ylabel(%22Effective%20population%20size%20(Ne)%22)%0A%20%20%20%20_ax_ne.set_title(%22Historical%20Ne%20Trajectory%20(LD-based%2C%20Sved%201971)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax_ne.legend(fontsize%3D10%2C%20frameon%3DFalse)%0A%20%20%20%20_ax_ne.grid(alpha%3D0.3%2C%20linestyle%3D%22--%22)%0A%20%20%20%20_ax_ne.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_het2%7D%5Cn%5Cn%7B_cmd_frq2%7D%5Cn%5Cn%7B_cmd_ld2%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_h)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%23%20Genetic%20diversity%20summary%0A%20%20%20%20%20%20%20%20Ho%2C%20observed%20heterozygosity%3B%20He%2C%20expected%20heterozygosity%3B%0A%20%20%20%20%20%20%20%20F_IS%20%3D%20(He%20%E2%88%92%20Ho)%20%2F%20He%20(Wright's%20inbreeding%20coefficient%20within%20populations).%0A%20%20%20%20%20%20%20%20%22%22%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_div_df)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Historical%20Ne%20trajectory%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22Ne%20estimated%20at%20multiple%20time%20points%20using%20Sved's%20(1971)%20formula%20from%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22mean%20pairwise%20r%C2%B2%20in%20each%20distance%20class%20(1%20Mb%20window).%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_ne)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%2011%20%C2%B7%20Pairwise%20F_ST%20%2B%20Manhattan%20Plot%0A%0A%20%20%20%20**Genome-wide%20F_ST**%20(Weir%20%26%20Cockerham%201984)%20is%20estimated%20pairwise%20for%0A%20%20%20%20each%20of%20the%20three%20breed%20pairs%20using%20%60plink%20--fst%20--within%60.%0A%0A%20%20%20%20The%20**Manhattan%20plot**%20of%20per-SNP%20F_ST%20highlights%20genomic%20regions%20with%0A%20%20%20%20unusually%20high%20differentiation%20%E2%80%94%20potential%20signatures%20of%20selection%20or%0A%20%20%20%20local%20adaptation.%20%20SNPs%20with%20F_ST%20%3E%20mean%20%2B%203%20SD%20are%20highlighted.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20Path%2C%20mo%2C%20pd%2C%20plt%2C%20pop_df%2C%20qc_prefix%2C%20re%2C%20sns%2C%20subprocess)%3A%0A%20%20%20%20import%20pandas%20as%20_pd3%0A%0A%20%20%20%20_clust%20%3D%20OUT%2F%22sheep_clusters.txt%22%0A%20%20%20%20pop_df%5B%5B%22FID%22%2C%22IID%22%2C%22POP%22%5D%5D.to_csv(_clust%2Csep%3D%22%5Ct%22%2Cindex%3DFalse%2Cheader%3DFalse)%0A%0A%20%20%20%20_cmd_fst%20%3D%20(%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--keep%20keep_%3CA%3E_%3CB%3E.txt%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--fst%20%5C%5C%5Cn%20%20%20%20--within%20sheep_clusters.txt%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%20%20%20%20--out%20fst_%3CA%3E_%3CB%3E%22)%0A%0A%20%20%20%20_pairs%20%3D%20%5B(%22Harri%22%2C%22Naemi%22)%2C(%22Harri%22%2C%22Najdi%22)%2C(%22Naemi%22%2C%22Najdi%22)%5D%0A%20%20%20%20_rows%20%20%3D%20%5B%5D%0A%20%20%20%20for%20_a%2C_b%20in%20_pairs%3A%0A%20%20%20%20%20%20%20%20_kf%20%20%3D%20OUT%2Ff%22keep_%7B_a%7D_%7B_b%7D.txt%22%0A%20%20%20%20%20%20%20%20pop_df%5Bpop_df%5B%22POP%22%5D.isin(%5B_a%2C_b%5D)%5D%5B%5B%22FID%22%2C%22IID%22%5D%5D.to_csv(%0A%20%20%20%20%20%20%20%20%20%20%20%20_kf%2Csep%3D%22%5Ct%22%2Cindex%3DFalse%2Cheader%3DFalse)%0A%20%20%20%20%20%20%20%20_op%20%20%3D%20str(OUT%2Ff%22fst_%7B_a%7D_%7B_b%7D%22)%0A%20%20%20%20%20%20%20%20_lp%20%20%3D%20Path(_op%2B%22.log%22)%0A%20%20%20%20%20%20%20%20if%20not%20_lp.exists()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%22--keep%22%2Cstr(_kf)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--fst%22%2C%22--within%22%2Cstr(_clust)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--chr-set%22%2C%2226%22%2C%22--out%22%2C_op%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%20%20%20%20%20%20%20%20_lt%20%20%3D%20_lp.read_text()%20if%20_lp.exists()%20else%20%22%22%0A%20%20%20%20%20%20%20%20_m2%20%20%3D%20re.search(r%22Mean%20Fst%20estimate%3A%5Cs%2B(%5B%5Cd.%5D%2B)%22%2C_lt)%0A%20%20%20%20%20%20%20%20_w2%20%20%3D%20re.search(r%22Weighted%20Fst%20estimate%3A%5Cs%2B(%5B%5Cd.%5D%2B)%22%2C_lt)%0A%20%20%20%20%20%20%20%20_rows.append(%7B%22Pair%22%3Af%22%7B_a%7D%20vs%20%7B_b%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Mean%20F_ST%22%3Afloat(_m2.group(1))%20if%20_m2%20else%20float(%22nan%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Weighted%20F_ST%22%3Afloat(_w2.group(1))%20if%20_w2%20else%20float(%22nan%22)%7D)%0A%0A%20%20%20%20fst_df%20%3D%20_pd3.DataFrame(_rows)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20FST%20heatmap%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_pops%20%3D%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%0A%20%20%20%20_mat%20%20%3D%20%7Bp%3A%7Bq%3Afloat(%22nan%22)%20for%20q%20in%20_pops%7D%20for%20p%20in%20_pops%7D%0A%20%20%20%20for%20_%2C_row%20in%20fst_df.iterrows()%3A%0A%20%20%20%20%20%20%20%20_a2%2C_b2%20%3D%20%5Bx.strip()%20for%20x%20in%20_row%5B%22Pair%22%5D.split(%22%20vs%20%22)%5D%0A%20%20%20%20%20%20%20%20_mat%5B_a2%5D%5B_b2%5D%20%3D%20_row%5B%22Weighted%20F_ST%22%5D%0A%20%20%20%20%20%20%20%20_mat%5B_b2%5D%5B_a2%5D%20%3D%20_row%5B%22Weighted%20F_ST%22%5D%0A%20%20%20%20for%20_p%20in%20_pops%3A%0A%20%20%20%20%20%20%20%20_mat%5B_p%5D%5B_p%5D%20%3D%200.0%0A%20%20%20%20_hm%20%3D%20_pd3.DataFrame(_mat%2Cindex%3D_pops%2Ccolumns%3D_pops).astype(float)%0A%0A%20%20%20%20_fig_fh%2C%20_ax_fh%20%3D%20plt.subplots(figsize%3D(5.5%2C4.5))%0A%20%20%20%20sns.heatmap(_hm%2C%20annot%3DTrue%2C%20fmt%3D%22.4f%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cmap%3D%22RdBu%22%2C%20vmin%3D0%2C%20vmax%3Dfloat(fst_df%5B%22Weighted%20F_ST%22%5D.max())%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20linewidths%3D0.8%2C%20annot_kws%3D%7B%22fontsize%22%3A12%7D%2C%20ax%3D_ax_fh)%0A%20%20%20%20_ax_fh.set_title(%22Pairwise%20F_ST%20%E2%80%94%20Identity-by-State%20Relatedness%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D11%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20helper%3A%20build%20Manhattan%20PNG%20(bytes)%20for%20one%20FST%20file%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20import%20io%20as%20_io%0A%0A%20%20%20%20def%20_manhattan(fst_file%2C%20title)%3A%0A%20%20%20%20%20%20%20%20_df%20%3D%20pd.read_csv(fst_file%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20%20%20%20%20_df.columns%20%3D%20%5Bc.strip()%20for%20c%20in%20_df.columns%5D%0A%20%20%20%20%20%20%20%20_df%20%3D%20_df%5B_df%5B%22FST%22%5D.notna()%20%26%20(_df%5B%22FST%22%5D%20%3E%3D%200)%5D.copy()%0A%20%20%20%20%20%20%20%20_df%5B%22FST%22%5D%20%3D%20_df%5B%22FST%22%5D.clip(0%2C%201)%0A%20%20%20%20%20%20%20%20_df%20%3D%20_df.sort_values(%5B%22CHR%22%2C%20%22POS%22%5D)%0A%20%20%20%20%20%20%20%20_chrs%20%3D%20sorted(_df%5B%22CHR%22%5D.unique())%0A%20%20%20%20%20%20%20%20_pal%20%20%3D%20sns.color_palette(%22tab20%22%2C%20len(_chrs))%0A%20%20%20%20%20%20%20%20_offset%20%3D%200%3B%20_cum%20%3D%20%5B%5D%3B%20_tick_pos%20%3D%20%5B%5D%3B%20_tick_lab%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20_c%2C%20_col%20in%20zip(_chrs%2C%20_pal)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_s%20%3D%20_df%5B_df%5B%22CHR%22%5D%20%3D%3D%20_c%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_cum.extend(_s%5B%22POS%22%5D.values%20%2B%20_offset)%0A%20%20%20%20%20%20%20%20%20%20%20%20_tick_pos.append(_offset%20%2B%20(_s%5B%22POS%22%5D.max()%20-%20_s%5B%22POS%22%5D.min())%20%2F%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20_tick_lab.append(str(_c))%0A%20%20%20%20%20%20%20%20%20%20%20%20_offset%20%2B%3D%20_s%5B%22POS%22%5D.max()%20%2B%201_000_000%0A%20%20%20%20%20%20%20%20_df%5B%22cum_pos%22%5D%20%3D%20_cum%0A%20%20%20%20%20%20%20%20_thr%20%3D%20_df%5B%22FST%22%5D.mean()%20%2B%203%20*%20_df%5B%22FST%22%5D.std()%0A%20%20%20%20%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(14%2C%204)%2C%20dpi%3D80)%0A%20%20%20%20%20%20%20%20for%20_c%2C%20_col%20in%20zip(_chrs%2C%20_pal)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_s%20%3D%20_df%5B_df%5B%22CHR%22%5D%20%3D%3D%20_c%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.scatter(_s.loc%5B_s%5B%22FST%22%5D%20%3C%20%20_thr%2C%20%22cum_pos%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_s.loc%5B_s%5B%22FST%22%5D%20%3C%20%20_thr%2C%20%22FST%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%3D%5B_col%5D%20*%20(_s%5B%22FST%22%5D%20%3C%20_thr).sum()%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20s%3D3%2C%20alpha%3D0.5%2C%20linewidths%3D0%2C%20rasterized%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.scatter(_s.loc%5B_s%5B%22FST%22%5D%20%3E%3D%20_thr%2C%20%22cum_pos%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_s.loc%5B_s%5B%22FST%22%5D%20%3E%3D%20_thr%2C%20%22FST%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%3D%22red%22%2C%20s%3D12%2C%20alpha%3D0.9%2C%20linewidths%3D0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20zorder%3D5%2C%20rasterized%3DTrue)%0A%20%20%20%20%20%20%20%20_ax.axhline(_thr%2C%20color%3D%22red%22%2C%20linestyle%3D%22--%22%2C%20linewidth%3D1%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22mean%20%2B%203%20SD%20(%7B_thr%3A.3f%7D)%22)%0A%20%20%20%20%20%20%20%20_ax.set_xticks(_tick_pos)%3B%20_ax.set_xticklabels(_tick_lab%2C%20fontsize%3D7)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22Chromosome%22)%3B%20_ax.set_ylabel(%22F_ST%22)%0A%20%20%20%20%20%20%20%20_ax.set_title(title%2C%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax.legend(fontsize%3D9%2C%20frameon%3DFalse)%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%20%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20%20%20%20%20plt.tight_layout()%0A%20%20%20%20%20%20%20%20_buf%20%3D%20_io.BytesIO()%0A%20%20%20%20%20%20%20%20_fig.savefig(_buf%2C%20format%3D%22png%22%2C%20dpi%3D80%2C%20bbox_inches%3D%22tight%22)%0A%20%20%20%20%20%20%20%20plt.close(_fig)%0A%20%20%20%20%20%20%20%20_buf.seek(0)%0A%20%20%20%20%20%20%20%20return%20mo.image(_buf.read())%0A%0A%20%20%20%20_man_items%20%3D%20%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Per-SNP%20F_ST%20Manhattan%20plots%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22Red%20dots%20exceed%20mean%20%2B%203%20SD%20%E2%80%94%20candidate%20regions%20of%20strong%20differentiation.%22)%2C%0A%20%20%20%20%5D%0A%20%20%20%20for%20_a%2C%20_b%20in%20_pairs%3A%0A%20%20%20%20%20%20%20%20_man_items.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20_manhattan(OUT%20%2F%20f%22fst_%7B_a%7D_%7B_b%7D.fst%22%2C%20f%22Per-SNP%20F_ST%20%E2%80%94%20%7B_a%7D%20vs%20%7B_b%7D%22)%0A%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Command%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_fst%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20F_ST%20heatmap%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22F_ST%20values%20with%20**higher%20genetic%20similarity**%20in%20**red**%2C%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22**lower%20similarity**%20in%20**blue**.%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_fh)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(fst_df)%2C%0A%20%20%20%20%20%20%20%20*_man_items%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%2012%20%C2%B7%20Runs%20of%20Homozygosity%20(ROH)%20and%20F_ROH%0A%0A%20%20%20%20Parameters%20are%20aligned%20with%20the%20reference%20study%3A%0A%0A%20%20%20%20%7C%20Flag%20%7C%20Reference%20%7C%20Our%20analysis%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20%60--homozyg-snp%60%20%7C%2050%20%7C%2050%20%7C%0A%20%20%20%20%7C%20%60--homozyg-kb%60%20%7C%201000%20(1%20Mb)%20%7C%20**1000**%20(corrected%20from%201500)%20%7C%0A%20%20%20%20%7C%20%60--homozyg-gap%60%20%7C%20100%20kb%20%7C%20**100**%20(corrected%20from%201000)%20%7C%0A%20%20%20%20%7C%20%60--homozyg-window-missing%60%20%7C%202%20%7C%20**2**%20(added)%20%7C%0A%0A%20%20%20%20Robustness%20analyses%20repeat%20the%20ROH%20detection%20at%20**300%20kb**%20and%20**500%20kb**%0A%20%20%20%20minimum%20lengths%2C%20following%20Ceballos%20et%20al.%0A%0A%20%20%20%20%60%60%60%0A%20%20%20%20F_ROH%20%3D%20%CE%A3%20L_ROH%20%2F%20L_autosome%20%20%20%20%20(sheep%20autosomal%20genome%20%E2%89%88%202%2C600%20Mb)%0A%20%20%20%20%60%60%60%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20PLINK%2C%20mo%2C%20qc_prefix%2C%20subprocess)%3A%0A%20%20%20%20_cmd_roh%20%3D%20(%0A%20%20%20%20%20%20%20%20%22%23%20Primary%20analysis%20(%E2%89%A5%201%20Mb)%5Cn%22%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--homozyg%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-snp%2050%20%5C%5C%5Cn%20%20%20%20--homozyg-kb%201000%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-density%2050%20%5C%5C%5Cn%20%20%20%20--homozyg-gap%20100%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-window-missing%202%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--out%20sheep_roh2%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%22%23%20Robustness%3A%20300%20kb%5Cn%22%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--homozyg%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-snp%2030%20%5C%5C%5Cn%20%20%20%20--homozyg-kb%20300%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-density%2050%20%5C%5C%5Cn%20%20%20%20--homozyg-gap%20100%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-window-missing%202%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--out%20sheep_roh_300kb%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%22%23%20Robustness%3A%20500%20kb%5Cn%22%0A%20%20%20%20%20%20%20%20%22plink%20%5C%5C%5Cn%20%20%20%20--bfile%20sheep_qc%20%5C%5C%5Cn%20%20%20%20--homozyg%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-snp%2040%20%5C%5C%5Cn%20%20%20%20--homozyg-kb%20500%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-density%2050%20%5C%5C%5Cn%20%20%20%20--homozyg-gap%20100%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--homozyg-window-missing%202%20%5C%5C%5Cn%20%20%20%20--chr-set%2026%20%5C%5C%5Cn%22%0A%20%20%20%20%20%20%20%20%22%20%20%20%20--out%20sheep_roh_500kb%22%0A%20%20%20%20)%0A%0A%20%20%20%20for%20_tag%2C%20_snp%2C%20_kb%20in%20%5B(%22roh2%22%2C50%2C1000)%2C(%22roh_300kb%22%2C30%2C300)%2C(%22roh_500kb%22%2C40%2C500)%5D%3A%0A%20%20%20%20%20%20%20%20_f%20%3D%20OUT%2Ff%22sheep_%7B_tag%7D.hom.indiv%22%0A%20%20%20%20%20%20%20%20if%20not%20_f.exists()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20subprocess.run(%5BPLINK%2C%22--bfile%22%2Cqc_prefix%2C%22--homozyg%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--homozyg-snp%22%2Cstr(_snp)%2C%22--homozyg-kb%22%2Cstr(_kb)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--homozyg-density%22%2C%2250%22%2C%22--homozyg-gap%22%2C%22100%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--homozyg-window-missing%22%2C%222%22%2C%22--chr-set%22%2C%2226%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22--out%22%2Cstr(OUT%2Ff%22sheep_%7B_tag%7D%22)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20capture_output%3DTrue%2Ctext%3DTrue)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Commands%22)%2C%0A%20%20%20%20%20%20%20%20mo.md(f%22%60%60%60bash%5Cn%7B_cmd_roh%7D%5Cn%60%60%60%22)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20mo%2C%20pd%2C%20plt%2C%20pop_df%2C%20pop_palette%2C%20sns)%3A%0A%20%20%20%20_AUTOSOME_KB%20%3D%202_600_000%0A%0A%20%20%20%20def%20_load_roh(tag)%3A%0A%20%20%20%20%20%20%20%20_f%20%3D%20OUT%20%2F%20f%22sheep_%7Btag%7D.hom.indiv%22%0A%20%20%20%20%20%20%20%20_r%20%3D%20pd.read_csv(_f%2C%20sep%3Dr%22%5Cs%2B%22)%0A%20%20%20%20%20%20%20%20_r.columns%20%3D%20%5Bc.strip()%20for%20c%20in%20_r.columns%5D%0A%20%20%20%20%20%20%20%20_r%20%3D%20_r.merge(pop_df%5B%5B%22IID%22%2C%22POP%22%5D%5D%2C%20on%3D%22IID%22%2C%20how%3D%22left%22)%0A%20%20%20%20%20%20%20%20_r%5B%22F_ROH%22%5D%20%20%20%20%20%20%20%20%3D%20_r%5B%22KB%22%5D%20%20%20%20%2F%20_AUTOSOME_KB%0A%20%20%20%20%20%20%20%20_r%5B%22Total_ROH_Mb%22%5D%20%3D%20_r%5B%22KB%22%5D%20%20%20%20%2F%201_000%0A%20%20%20%20%20%20%20%20_r%5B%22Mean_ROH_Mb%22%5D%20%20%3D%20_r%5B%22KBAVG%22%5D%20%2F%201_000%0A%20%20%20%20%20%20%20%20return%20_r%0A%0A%20%20%20%20_roh_main%20%3D%20_load_roh(%22roh2%22)%0A%20%20%20%20_roh_300%20%20%3D%20_load_roh(%22roh_300kb%22)%0A%20%20%20%20_roh_500%20%20%3D%20_load_roh(%22roh_500kb%22)%0A%20%20%20%20_ord%20%20%20%20%20%20%3D%20%5B%22Harri%22%2C%22Naemi%22%2C%22Najdi%22%5D%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20main%20ROH%20plots%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig_r%2C%20_axs_r%20%3D%20plt.subplots(1%2C3%2C%20figsize%3D(15%2C5))%0A%20%20%20%20_fig_r.suptitle(%22ROH%20(%E2%89%A5%201%20Mb)%20per%20Individual%22%2C%20fontsize%3D12%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20for%20_ax%2C_col%2C_title%2C_yl%20in%20%5B%0A%20%20%20%20%20%20%20%20(_axs_r%5B0%5D%2C%22F_ROH%22%2C%20%20%20%20%20%20%22Inbreeding%20(F_ROH)%22%2C%20%20%20%20%20%20%20%22F_ROH%22)%2C%0A%20%20%20%20%20%20%20%20(_axs_r%5B1%5D%2C%22Total_ROH_Mb%22%2C%22Total%20ROH%20length%20(Mb)%22%2C%20%20%20%20%22Mb%22)%2C%0A%20%20%20%20%20%20%20%20(_axs_r%5B2%5D%2C%22NSEG%22%2C%20%20%20%20%20%20%20%20%22N%20ROH%20segments%22%2C%20%20%20%20%20%20%20%20%20%20%20%22N%22)%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20sns.boxplot(data%3D_roh_main%2Cx%3D%22POP%22%2Cy%3D_col%2Cpalette%3Dpop_palette%2Cax%3D_ax%2Corder%3D_ord)%0A%20%20%20%20%20%20%20%20sns.stripplot(data%3D_roh_main%2Cx%3D%22POP%22%2Cy%3D_col%2Cpalette%3Dpop_palette%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20jitter%3DTrue%2Calpha%3D0.5%2Cax%3D_ax%2Csize%3D5%2Corder%3D_ord)%0A%20%20%20%20%20%20%20%20_ax.set_title(_title)%3B%20_ax.set_xlabel(%22Breed%22)%3B%20_ax.set_ylabel(_yl)%0A%20%20%20%20%20%20%20%20_ax.grid(axis%3D%22y%22%2Calpha%3D0.3%2Clinestyle%3D%22--%22)%0A%20%20%20%20%20%20%20%20_ax.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20robustness%3A%20F_ROH%20at%20300%20%2F%20500%20%2F%201000%20kb%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20_fig_rb%2C%20_ax_rb%20%3D%20plt.subplots(figsize%3D(9%2C4))%0A%20%20%20%20_colors_rb%20%3D%20%7B%22Harri%22%3A%22%23E63946%22%2C%22Naemi%22%3A%22%23457B9D%22%2C%22Najdi%22%3A%22%232A9D8F%22%7D%0A%20%20%20%20_markers%20%20%20%3D%20%7B%22300%20kb%22%3A%22%5E%22%2C%22500%20kb%22%3A%22s%22%2C%221000%20kb%22%3A%22o%22%7D%0A%20%20%20%20for%20_breed%20in%20_ord%3A%0A%20%20%20%20%20%20%20%20for%20_roh_df%2C%20_label%20in%20%5B(_roh_300%2C%22300%20kb%22)%2C(_roh_500%2C%22500%20kb%22)%2C(_roh_main%2C%221000%20kb%22)%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_val%20%3D%20_roh_df%5B_roh_df%5B%22POP%22%5D%3D%3D_breed%5D%5B%22F_ROH%22%5D.mean()%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax_rb.scatter(%5B_label%5D%2C%5B_val%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20color%3D_colors_rb%5B_breed%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20marker%3D_markers%5B_label%5D%2C%20s%3D80%2C%20zorder%3D5)%0A%20%20%20%20import%20matplotlib.lines%20as%20_mlines%0A%20%20%20%20_hb%20%3D%20%5B_mlines.Line2D(%5B%5D%2C%5B%5D%2C%20marker%3D%22o%22%2C%20color%3D%22w%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20markerfacecolor%3Dc%2C%20markersize%3D9%2C%20label%3Db)%0A%20%20%20%20%20%20%20%20%20%20%20for%20b%2Cc%20in%20_colors_rb.items()%5D%0A%20%20%20%20_ax_rb.legend(handles%3D_hb%2C%20title%3D%22Breed%22%2C%20fontsize%3D9)%0A%20%20%20%20_ax_rb.set_xlabel(%22Minimum%20ROH%20length%22)%3B%20_ax_rb.set_ylabel(%22Mean%20F_ROH%22)%0A%20%20%20%20_ax_rb.set_title(%22Robustness%20of%20F_ROH%20Across%20Minimum%20ROH%20Length%20Thresholds%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D11%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20_ax_rb.grid(alpha%3D0.3%2Clinestyle%3D%22--%22)%3B%20_ax_rb.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20summary%20tables%20%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%E2%94%80%0A%20%20%20%20def%20_roh_tbl(df%2C%20label)%3A%0A%20%20%20%20%20%20%20%20_t%20%3D%20(df.groupby(%22POP%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20.agg(N%3D(%22IID%22%2C%22count%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Total_ROH_Mb%3D(%22Total_ROH_Mb%22%2C%22mean%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Mean_ROH_Mb%3D(%22Mean_ROH_Mb%22%2C%22mean%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20F_ROH%3D(%22F_ROH%22%2C%22mean%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20N_ROH%3D(%22NSEG%22%2C%22mean%22))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20.round(4).reset_index()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20.rename(columns%3D%7B%22POP%22%3A%22Breed%22%2C%22Total_ROH_Mb%22%3A%22Mean%20Total%20ROH%20(Mb)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Mean_ROH_Mb%22%3A%22Mean%20ROH%20length%20(Mb)%22%2C%22F_ROH%22%3A%22Mean%20F_ROH%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22N_ROH%22%3A%22Mean%20N_ROH%22%7D))%0A%20%20%20%20%20%20%20%20_t.insert(0%2C%22Min%20ROH%22%2C%20label)%0A%20%20%20%20%20%20%20%20return%20_t%0A%0A%20%20%20%20import%20pandas%20as%20_pd4%0A%20%20%20%20_tbl_all%20%3D%20_pd4.concat(%5B_roh_tbl(_roh_300%2C%22300%20kb%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_roh_tbl(_roh_500%2C%22500%20kb%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_roh_tbl(_roh_main%2C%221000%20kb%22)%5D%2C%20ignore_index%3DTrue)%0A%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_r)%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20ROH%20summary%20%E2%80%94%20main%20analysis%20(%E2%89%A5%201%20Mb)%22)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_roh_tbl(_roh_main%2C%221000%20kb%22))%2C%0A%20%20%20%20%20%20%20%20mo.md(%22%23%23%23%20Robustness%20%E2%80%94%20F_ROH%20across%20minimum%20ROH%20length%20thresholds%22)%2C%0A%20%20%20%20%20%20%20%20mo.as_html(_fig_rb)%2C%0A%20%20%20%20%20%20%20%20mo.ui.table(_tbl_all)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%0A%20%20%20%20%23%23%2013%20%C2%B7%20Per-Chromosome%20SNP%20Density%20(post-QC)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OUT%2C%20mo%2C%20pd%2C%20plt%2C%20sns)%3A%0A%20%20%20%20_bim%20%3D%20pd.read_csv(OUT%2F%22sheep_qc.bim%22%2C%20sep%3D%22%5Ct%22%2C%20header%3DNone%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20names%3D%5B%22CHR%22%2C%22SNP%22%2C%22CM%22%2C%22BP%22%2C%22A1%22%2C%22A2%22%5D)%0A%20%20%20%20_cnt%20%3D%20_bim.groupby(%22CHR%22).size().reset_index(name%3D%22N_SNPs%22).sort_values(%22CHR%22)%0A%20%20%20%20_fig7%2C_ax7%20%3D%20plt.subplots(figsize%3D(12%2C4))%0A%20%20%20%20_ax7.bar(_cnt%5B%22CHR%22%5D.astype(str)%2C_cnt%5B%22N_SNPs%22%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20color%3Dsns.color_palette(%22viridis%22%2Clen(_cnt)))%0A%20%20%20%20_ax7.set_xlabel(%22Chromosome%22)%3B%20_ax7.set_ylabel(%22SNPs%20after%20QC%22)%0A%20%20%20%20_ax7.set_title(%22Post-QC%20SNP%20Density%20per%20Chromosome%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D12%2Cfontweight%3D%22bold%22)%0A%20%20%20%20_ax7.grid(axis%3D%22y%22%2Calpha%3D0.3%2Clinestyle%3D%22--%22)%0A%20%20%20%20_ax7.spines%5B%5B%22top%22%2C%22right%22%5D%5D.set_visible(False)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20mo.vstack(%5Bmo.as_html(_fig7)%2C%20mo.ui.table(_cnt)%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Summary%20%E2%80%94%20Comparison%20with%20Reference%20Study%0A%0A%20%20%20%20%7C%20Step%20%7C%20Reference%20(OvineSNP50)%20%7C%20This%20analysis%20(OvineSNPIII)%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20Raw%20SNPs%20%7C%2064%2C734%20%7C%2053%2C370%20%7C%0A%20%20%20%20%7C%20After%20missingness%20QC%20%7C%2062%2C702%20%7C%20%E2%80%94%20(VCF%20pre-processed)%20%7C%0A%20%20%20%20%7C%20After%20MAF%20%E2%89%A5%200.01%20%7C%2058%2C116%20%7C%2051%2C457%20%7C%0A%20%20%20%20%7C%20Heterozygosity%20outlier%20%7C%20H16%20removed%20%E2%86%92%2095%20samples%20%7C%20No%20outliers%2C%2096%20samples%20%7C%0A%20%20%20%20%7C%20After%20LD%20pruning%20(r%C2%B2%20%3D%200.3)%20%7C%2034%2C079%20%7C%20**40%2C930**%20%7C%0A%20%20%20%20%7C%20MAF%20distribution%20%7C%20%E2%9C%93%20(4%20bins%20per%20breed)%20%7C%20%E2%9C%93%207%20%7C%0A%20%20%20%20%7C%20PCA%20%7C%20%E2%9C%93%20GRM-based%20%7C%20%E2%9C%93%208%20%7C%0A%20%20%20%20%7C%20Structure%20%7C%20sNMF%20K%3D1%E2%80%9310%2C%20cross-entropy%20%7C%20NMF%20K%3D1%E2%80%936%2C%20reconstruction%20error%209%20%7C%0A%20%20%20%20%7C%20Ho%20%C2%B7%20He%20%C2%B7%20F_IS%20%7C%20%E2%9C%93%20%7C%20%E2%9C%93%2010%20%7C%0A%20%20%20%20%7C%20Ne%20trajectory%20%7C%20%E2%9C%93%20LD%201%20Mb%20window%20%7C%20%E2%9C%93%2010%20%7C%0A%20%20%20%20%7C%20F_ST%20%2B%20Manhattan%20%7C%20%E2%9C%93%20heatmap%20%2B%20plot%20%7C%20%E2%9C%93%2011%20%7C%0A%20%20%20%20%7C%20ROH%20(%E2%89%A5%201%20Mb%2C%20gap%20100%20kb)%20%7C%20%E2%9C%93%20%2B%20robustness%20300%2F500%20kb%20%7C%20%E2%9C%93%2012%20%7C%0A%20%20%20%20%7C%20IBD%20(Beagle%20%2B%20Refined%20IBD)%20%7C%20%E2%9C%93%20%7C%20%E2%9C%97%20requires%20phased%20data%20%7C%0A%20%20%20%20%7C%20Phylogenetic%20network%20%7C%20%E2%9C%93%20SplitsTree%20%7C%20%E2%9C%97%20requires%20SplitsTree%20%7C%0A%0A%20%20%20%20%3E%20**Breeds%3A**%20Harri%20(H%2C%20n%20%3D%2039)%20%C2%B7%20Naemi%20(M%2C%20n%20%3D%2024)%20%C2%B7%20Najdi%20(N%2C%20n%20%3D%2033)%0A%20%20%20%20%3E%20All%20files%20in%20%60plink_analysis%2F%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
d50d2249c3366fd4efd95dfcbca6278e