使用iqtree基于SNP构建种群/物种树

1. 格式转换
   pop.list的第一列是个体名、第二列是种群名，tab分隔符.

python3 vcf.to.pomo.py  in.vcf  pop.list  pomo.cf

cat vcf.to.pomo.py

#!/usr/bin/env python3

import argparse
import gzip
import os
import shutil
import sys
from collections import OrderedDict

BASES = "ACGT"
BASE_INDEX = {b: i for i, b in enumerate(BASES)}


def open_text(path):
    if path.endswith(".gz"):
        return gzip.open(path, "rt")
    return open(path, "r")


def read_popmap(path):
    sample_to_pop = OrderedDict()
    pop_order = OrderedDict()

    with open(path) as f:
        for line in f:
            if not line.strip() or line.startswith("#"):
                continue
            parts = line.split()
            if len(parts) < 2:
                continue
            sample, pop = parts[0], parts[1]
            sample_to_pop[sample] = pop
            pop_order.setdefault(pop, 0)
            pop_order[pop] += 1

    if not sample_to_pop:
        sys.exit("ERROR: pop list is empty or incorrectly formatted.")

    return sample_to_pop, list(pop_order.keys())


def parse_gt(sample_field):
    gt = sample_field.split(":", 1)[0]
    if gt in {".", "./.", ".|."}:
        return []
    gt = gt.replace("|", "/")
    alleles = []
    for x in gt.split("/"):
        if x == ".":
            continue
        try:
            alleles.append(int(x))
        except ValueError:
            continue
    return alleles


def main():
    parser = argparse.ArgumentParser(
        description="Convert a biallelic SNP VCF and a sample to population map into a PoMo counts file."
    )
    parser.add_argument("vcf", help="Input VCF, plain text or gzipped")
    parser.add_argument("popmap", help="Two column file: sample population")
    parser.add_argument("output", help="Output PoMo counts file")
    parser.add_argument(
        "--allow-extra-vcf-samples",
        action="store_true",
        help="Ignore VCF samples that are not present in the popmap"
    )
    args = parser.parse_args()

    sample_to_pop, pop_order = read_popmap(args.popmap)
    tmp = args.output + ".tmp"

    n_sites = 0
    header_seen = False

    with open_text(args.vcf) as fin, open(tmp, "w") as fout:
        for line in fin:
            if line.startswith("##"):
                continue

            if line.startswith("#CHROM"):
                header_seen = True
                fields = line.rstrip("\n").split("\t")
                vcf_samples = fields[9:]

                vcf_sample_set = set(vcf_samples)
                popmap_sample_set = set(sample_to_pop)

                missing_in_vcf = sorted(popmap_sample_set - vcf_sample_set)
                if missing_in_vcf:
                    sys.exit(
                        "ERROR: These popmap samples are absent from the VCF:\n"
                        + "\n".join(missing_in_vcf[:20])
                    )

                extra_in_vcf = sorted(vcf_sample_set - popmap_sample_set)
                if extra_in_vcf and not args.allow_extra_vcf_samples:
                    sys.exit(
                        "ERROR: These VCF samples are absent from the popmap:\n"
                        + "\n".join(extra_in_vcf[:20])
                        + "\nUse --allow-extra-vcf-samples to ignore them."
                    )

                sample_pop_pairs = []
                for i, sample in enumerate(vcf_samples, start=9):
                    if sample in sample_to_pop:
                        sample_pop_pairs.append((i, sample_to_pop[sample]))
                continue

            if not header_seen:
                continue

            fields = line.rstrip("\n").split("\t")
            if len(fields) < 10:
                continue

            chrom = fields[0]
            pos = fields[1]
            ref = fields[3].upper()
            alt = fields[4].upper()

            if "," in alt:
                continue
            if ref not in BASE_INDEX or alt not in BASE_INDEX:
                continue
            if len(ref) != 1 or len(alt) != 1:
                continue

            allele_to_base = [ref, alt]
            counts = {pop: [0, 0, 0, 0] for pop in pop_order}

            for col_i, pop in sample_pop_pairs:
                for allele_i in parse_gt(fields[col_i]):
                    if allele_i >= len(allele_to_base):
                        continue
                    base = allele_to_base[allele_i]
                    counts[pop][BASE_INDEX[base]] += 1

            row_counts = []
            bad_site = False
            for pop in pop_order:
                c = counts[pop]
                if sum(c) == 0:
                    bad_site = True
                    break
                row_counts.append(",".join(str(x) for x in c))

            if bad_site:
                continue

            fout.write(chrom + "\t" + pos + "\t" + "\t".join(row_counts) + "\n")
            n_sites += 1

    if n_sites == 0:
        if os.path.exists(tmp):
            os.remove(tmp)
        sys.exit("ERROR: No usable biallelic A/C/G/T SNP sites were written.")

    with open(args.output, "w") as out, open(tmp, "r") as inp:
        out.write(f"COUNTSFILE\tNPOP\t{len(pop_order)}\tNSITES\t{n_sites}\n")
        out.write("CHROM\tPOS\t" + "\t".join(pop_order) + "\n")
        shutil.copyfileobj(inp, out)

    os.remove(tmp)

    sys.stderr.write(f"Done. Populations: {len(pop_order)}. Sites: {n_sites}.\n")


if __name__ == "__main__":
    main()

2. 运行iqtree

iqtree3  -s pomo.cf  -m GTR+P+N19+WB  -B 1000  -T AUTO  --prefix output