<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <h2><span style="color
Step1:  Development of the Processor class to calculate all the stats
Step2:  Prototyping the dcons function to split alleles per base
Step3:  Loading pickled results
Step4:  Prototype Dxy
Step5:  Prototype Fst
Step6:  Prototyping collating stats across runs
<ASSISTANT_TASK:>
Python Code:
!hostname
%load_ext autoreload
%autoreload 2
%matplotlib inline
import ipyrad
import ipyrad.analysis as ipa
import ipyparallel as ipp
from ipyrad.analysis.popgen import Popgen
from ipyrad import Assembly
from ipyrad.analysis.locus_extracter import LocusExtracter
ipyclient = ipp.Client(cluster_id="popgen")
print(len(ipyclient))
# popgen tools can accept either an ipyrad assembly
data = ipyrad.load_json("/tmp/ipyrad-test/watdo.json")
# or alternatively the path to your VCF or HDF5 formatted snps file
#data = "/tmp/ipyrad-test/watdo_outfiles/watdo.snps.hdf5"
imap = {
    "pop1" : ["1A_0", "1B_0", "1C_0", "1D_0"],
    "pop2" : ["2E_0", "2F_0", "2G_0", "2H_0"],
    "pop3" : ["3I_0", "3J_0", "3K_0", "3L_0"],
}
popgen = Popgen(data=data, imap=imap)
popgen.samples
popgen.params
from IPython.display import display
popgen.run(ipyclient=ipyclient)
popgen.results
display(popgen.results)
from ipyrad.analysis.locus_extracter import LocusExtracter
import ipyparallel as ipp
ipyclient = ipp.Client(cluster_id="popgen")
print(len(ipyclient))
lex = LocusExtracter(
            data=data.seqs_database,
            imap=imap,
            mincov=len(imap),  # ENFORCE at least 1 per spp.
        )
lex.run(ipyclient=ipyclient)
print(len(popgen.lex.loci))
popgen.lex.get_locus(1, as_df=True)
import pandas as pd
wat = pd.DataFrame()
with h5py.File(data.snps_database, 'r') as io5:
    diffs = io5["snps"][0] != io5["snps"][1]
    for idx, name in enumerate(io5["snps"].attrs["names"]):
        wat[name.decode("utf-8")] = io5["snps"][idx]
wat["1A_0"]
import h5py
with h5py.File(data.seqs_database, 'r') as io5:
    print(io5.keys())
    print(io5["phymap"].attrs.keys())
    print(io5["phymap"].attrs["phynames"])
    print(io5["phy"][0])
Explanation: <h2><span style="color:gray">ipyrad-analysis toolkit:</span> Popgen summary statistics</h2>
Calculate summary statistics such as pi, Tajima's D, Fst
End of explanation
from collections import Counter
from itertools import combinations
import numpy as np
# Make a processor and give it some data
loci = [lex.get_locus(x, as_df=True) for x in range(2)]
proc = ipa.popgen.Processor(popgen.params, 0, loci)
proc.run()
#proc._pi()
locus = loci[0]
#locus[10][:5] = 82
#display(locus)
#%timeit proc.pi(locus)
print(proc.pi(locus))
print(proc.Watterson(locus))
print(proc.TajimasD(locus))
print(proc.results)
Explanation: Development of the Processor class to calculate all the stats
End of explanation
from ipyrad.assemble.utils import DCONS
import pandas as pd
import itertools
p1 = popgen.imap["pop1"]
#locus.loc[p1, :].apply(lambda x: [DCONS[y] for y in x])
cts = np.array(locus.apply(lambda bases:\
                Counter(x for x in bases if x not in [45, 78])))
snps = np.array([len(x) for x in cts]) > 1
cts = cts[snps]
def dcons(counter):
    new = list(itertools.chain(*[DCONS[x]*ct for x, ct in counter.items()]))
    return Counter(new)
print(cts)
%timeit list(map(dcons, cts))
Explanation: Prototyping the dcons function to split alleles per base
End of explanation
ipyrad.analysis.popgen._calc_sumstats(popgen, 10, loci)
import pickle
!ls analysis-popgen/
with open("analysis-popgen/0.p", 'rb') as inp:
    dat = pickle.load(inp)
dat
Explanation: Loading pickled results
End of explanation
proc._process_locus_pops(locus, ["pop1", "pop3"])
pop_cts, sidxs = proc._process_locus_pops(locus, ["pop1", "pop2"])
# Between population summary statistics
def _dxy(cts_a, cts_b):
    Dxy = 0
    ncomps = 0
    for cta, ctb in zip(cts_a, cts_b):
        ncomps += sum(list(cta.values())) *\
                    sum(list(ctb.values()))
        for ka, va in cta.items():
            for kb, vb in ctb.items():
                if ka == kb: continue
                Dxy += va*vb
            print(Dxy, ncomps)
    return Dxy/ncomps
Dxy = _dxy(pop_cts["pop1"], pop_cts["pop2"])
Dxy/len(locus)
%timeit proc.Dxy(locus, ["pop1", "pop2"])
Explanation: Prototype Dxy
End of explanation
proc._fst_full(locus)
print(np.zeros(len(proc.data.imap), len(proc.data.imap)))
Dxy_arr = pd.DataFrame(
    data=np.zeros(len(proc.data.imap), len(proc.data.imap)),
    index=proc.data.imap.keys(),
    columns=proc.data.imap.keys(),
)
loci = [lex.get_locus(x, as_df=True) for x in range(100)]
proc = ipa.popgen.Processor(popgen.params, 0, loci)
proc.run()
Explanation: Prototype Fst
End of explanation
import glob
pickles = glob.glob(os.path.join(popgen.workdir, "*.p"))
sorted(pickles, key=lambda x: int(x.rsplit("/", 1)[-1][:-2]))
#pickles[0].rsplit("/", 1)[-1][:-2]
pdicts = {}
for pkl in pickles:
    with open(pkl, 'rb') as inp:
        pdicts[int(pkl.rsplit("/", 1)[-1][:-2])] = pickle.load(inp)
pdicts[0]["pi"]
#print(pdicts[0]["pi"])
pdicts[0]["Fst"].keys()
full_res = {}
for d in [pdicts]: full_res.update(d)
full_res.keys()
pidx = sorted(full_res.keys())
pi_dict = {}
w_theta_dict = {}
tajd_dict = {}
for idx in pidx:
    pi_dict.update(full_res[idx]["pi"])
    w_theta_dict.update(full_res[idx]["Watterson"])
    tajd_dict.update(full_res[idx]["TajimasD"])
popstats = {}
for pop in proc.imap:
    popstats[pop] = pd.DataFrame([], columns=["pi",
                                                "raw_pi",
                                                "Watterson",
                                                "raw_Watterson",
                                                "TajimasD"], index=range(len(popgen.lex.loci)))
    for lidx in range(len(popgen.lex.loci)):
        popstats[pop]["pi"].loc[lidx] = pi_dict[lidx][pop]["pi_per_base"]
        popstats[pop]["raw_pi"].loc[lidx] = pi_dict[lidx][pop]["pi"]
        popstats[pop]["Watterson"].loc[lidx] = w_theta_dict[lidx][pop]["w_theta_per_base"]
        popstats[pop]["raw_Watterson"].loc[lidx] = w_theta_dict[lidx][pop]["w_theta"]
        popstats[pop]["TajimasD"].loc[lidx] = tajd_dict[lidx][pop]
                
lidx = sorted(full_res.keys())
for idx in lidx[:1]:
    for pop in proc.imap:
        for bidx in full_res[idx]["pi"]:
            print(full_res[idx]["pi"][bidx][pop]["pi_per_base"])
#        pi_per_base = np.mean(full_res[idx]["pi"][idx][pop]["pi_per_base"])
#        print(pop, pi_per_base)
pi_dict[0]
#popstats["pop1"].mean()
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
df = pd.read_csv("/tmp/gorg-tropics_sags_tableS2 - SAGs.csv")
print(df.columns)
set(df["Lineage"])
#df[["Genome completeness (%)", "Lineage"]]
df[["Raw read count", "Lineage"]]
plt.hist(df[df["Lineage"] == "AEGEAN-169"]["Genome completeness (%)"])
import locale
locale.setlocale(locale.LC_ALL, 'en_US.UTF8')
fig, ax = plt.subplots(figsize=(10, 10))
for l in set(df["Lineage"]):
    #print(l, np.mean(df[df["Lineage"] == l]["Genome completeness (%)"]))
    #print(l, np.mean([locale.atoi(x) for x in df[df["Lineage"] == l]["Raw read count"]]))
    #print(l, np.std([locale.atoi(x) for x in df[df["Lineage"] == l]["Assembly size (bp)"]]))
    lmask = df[np.array(df["Lineage"] == l) + np.array(df["Genome completeness (%)"]>80)]
#        cmask = df[df["Genome completeness (%)" > 80]]
    try:
        alpha=0.05
        if l == "AEGEAN-169": alpha=1
        plt.hist(lmask["Genome completeness (%)"], alpha=alpha, label=l, bins=40)    
    except: pass
plt.xlim(80, 100)
plt.legend()
    
nsamps=10
nspecies=3
dfs = []
for idx in range(nspecies):
    df = pd.DataFrame([
                [idx] * nsamps,
                range(nsamps),
                np.random.normal(0, 10, nsamps),
                np.random.normal(0, 1, nsamps),
                np.random.randint(0, 100, nsamps),
                np.random.choice(["small", "medium", "large"], nsamps),
                ],
                index=["Species_ID", "Sample_ID", "Trait1", "Trait2", "Trait3", "Trait4"]).T
    dfs.append(df)
df = pd.concat(dfs)
df.to_csv("/tmp/watdo.csv", index=False)
!cat /tmp/watdo.csv
pd.set_option('display.max_rows', 999)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
inf = "/home/isaac/Continuosity/NEON/NEON_seq-metabarcode-zooplankton/NEON.D03.BARC.DP1.20221.001.2019-07.expanded.20210123T023002Z.RELEASE-2021/NEON.D03.BARC.DP1.20221.001.zoo_metabarcodeTaxonomy.2019-07.expanded.20201218T153238Z.csv"
df = pd.read_csv(inf)
df
Explanation: Prototyping collating stats across runs
End of explanation
<END_TASK>