a peak around 10³ in phlash results

[YakunWang]

Hi
I'm trying phlash on 5–6 different plant species (both trees and herbs), and I noticed something interesting. In every dataset, there is a very obvious peak around (10^3) in the result (see the attached figure). Is this expected, or could it be caused by some common factor in the data or analysis? I'm wondering what might explain this pattern.
I also have a question about VCF filtering. For phlash, is it enough to apply the basic quality filters and keep only biallelic SNPs? Do you recommend filtering by MAF ？
Here's the command I used:

import os
import logging
import pickle
from collections import defaultdict
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import phlash
os.environ["CUDA_VISIBLE_DEVICES"] = "0"
logging.basicConfig(
level=logging.INFO,
format="%(asctime)s - %(levelname)s - %(message)s"
)

onekg_base = "/mnt/raid1/wyk/ulmus/phlash/"
template = "filtered_391samples.snp.vcf.gz"
pop_file = "pop_info.txt"
fai_file = "../snp/uelo.chr.fa.fai"
pop_dict = defaultdict(list)
with open(pop_file) as f:
for line in f:
if line.strip():
sample, pop = line.strip().split()
pop_dict[pop].append(sample)

logging.info(f"读取到 {len(pop_dict)} 个群体")

fai = pd.read_csv(
fai_file,
sep="\t",
header=None,
names=["chrom", "length", "offset", "linebases", "linewidth"],
)
chr_len_dict = dict(zip(fai["chrom"], fai["length"]))

chrom_range = range(1, 15)
chrom_prefix = "Chr" #

def run_population(pop, samples):
logging.info(f"Processing population: {pop} ({len(samples)} samples)")

chroms_1kg = []
for chrom in chrom_range:
    chrom_str = f"{chrom_prefix}{chrom:02d}"
    length = chr_len_dict.get(chrom_str)
    if length is None:
        logging.warning(f"{chrom_str} not found in fai file, skipping")
        continue

    region_str = f"{chrom_str}:1-{length}"
    path = os.path.join(onekg_base, template)

    logging.info(f"Loading {region_str} from {path}")
    contig = phlash.contig(path, samples=samples, region=region_str)
    chroms_1kg.append(contig)
if not chroms_1kg:
    logging.error(f"No valid chromosomes for {pop}, skipping.")
    return
with open(f"{pop}_chroms_1kg.pkl", "wb") as f:
    pickle.dump(chroms_1kg, f)

logging.info(f"Fitting model for {pop}")
results = phlash.fit(
    chroms_1kg,
    mutation_rate=3.75e-8,
    niter=1000
)
with open(f"{pop}_phlash_results.pkl", "wb") as f:
    pickle.dump(results, f)

times = np.array([dm.eta.t[1:] for dm in results])
T = np.geomspace(times.min(), times.max(), 1000)
Nes = np.array([dm.eta(T, Ne=True) for dm in results])

generation_times = [6, 10, 1]
for gen in generation_times:
    time_in_years = T * gen
    plt.figure(figsize=(6, 4))
    plt.plot(time_in_years, np.median(Nes, axis=0), color="blue")
    plt.xscale("log")
    plt.yscale("log")
    plt.xlabel("Time (years)")
    plt.ylabel("Effective population size (Ne)")
    plt.title(f"Ne history ({gen} years) - {pop}")
    plt.tight_layout()
    plt.savefig(f"Ne_history_{gen}_years_{pop}.pdf")
    plt.close()

df_dict = {"Time_generations": T, "Ne_median": np.median(Nes, axis=0)}
for gen in generation_times:
    df_dict[f"Time_years_{gen}"] = T * gen
df = pd.DataFrame(df_dict)
df.to_csv(f"Ne_history_data_{pop}.csv", index=False)

logging.info(f"✅ {pop} completed.")

if name == "main":
for pop, samples in pop_dict.items():

    if len(samples) == 0:
        logging.warning(f"Skipping population {pop}, no samples found")
        continue

    result_file = f"{pop}_phlash_results.pkl"

    if os.path.exists(result_file):
        logging.info(f"✓ {pop} already finished, skipping.")
        continue

    try:
        run_population(pop, samples)
    except Exception as e:
        logging.exception(f"❌ {pop} failed: {e}")
        continue

logging.info("All populations processed.")