20261002

99.scripts_v2/utils/asm_qc.sh

@@ -0,0 +1,22 @@
+#!/usr/bin/env bash
+# 轻量组装 QC：contig 数 / 总长度 / N50 / N 占比 / GC%
+# 用法: bash asm_qc.sh <fasta文件或目录...> > qc.tsv
+set -euo pipefail
+
+files=()
+for p in "$@"; do
+    if [ -d "$p" ]; then
+        while IFS= read -r f; do files+=("$f"); done \
+            < <(find "$p" -type f \( -name '*.fa' -o -name '*.fasta' -o -name '*.fna' -o -name '*.fa.gz' -o -name '*.fasta.gz' -o -name '*.fna.gz' \))
+    else
+        files+=("$p")
+    fi
+done
+
+seqkit stats -a -T -j 8 "${files[@]}" | awk 'BEGIN{OFS="\t"}
+NR==1 {print $0, "N%"; next}
+{
+    # $5=sum_len, $19=sum_n（N 碱基数）
+    n = ($5 > 0) ? sprintf("%.4f", $19 / $5 * 100) : "NA"
+    print $0, n
+}'

99.scripts_v2/utils/classify_frameshifts.py

@@ -0,0 +1,269 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+from collections import Counter
+from pathlib import Path
+
+
+LOCUS_HEADER = [
+    "OG",
+    "sequence_count",
+    "independent_event_count",
+    "max_single_sequence_unique_events",
+    "carrier_pattern",
+    "severe",
+]
+
+
+def read_fasta(fasta_path: Path) -> dict[str, str]:
+    sequences: dict[str, list[str]] = {}
+    sequence_id: str | None = None
+
+    with fasta_path.open() as fasta:
+        for line_number, raw_line in enumerate(fasta, start=1):
+            line = raw_line.strip()
+            if not line:
+                continue
+            if line.startswith(">"):
+                sequence_id = line[1:].split()[0]
+                if not sequence_id:
+                    raise ValueError(
+                        f"{fasta_path}:{line_number}: empty FASTA identifier"
+                    )
+                if sequence_id in sequences:
+                    raise ValueError(
+                        f"{fasta_path}:{line_number}: duplicate FASTA identifier "
+                        f"{sequence_id!r}"
+                    )
+                sequences[sequence_id] = []
+            elif sequence_id is None:
+                raise ValueError(
+                    f"{fasta_path}:{line_number}: sequence before first FASTA header"
+                )
+            else:
+                sequences[sequence_id].append(line)
+
+    if not sequences:
+        raise ValueError(f"{fasta_path}: no FASTA sequences found")
+
+    joined = {name: "".join(parts) for name, parts in sequences.items()}
+    lengths = {len(sequence) for sequence in joined.values()}
+    if len(lengths) != 1:
+        raise ValueError(f"{fasta_path}: aligned sequences have different lengths")
+    return joined
+
+
+def event_mode(carrier_count: int, sequence_count: int) -> str:
+    if carrier_count == 1:
+        return "single_species"
+    if carrier_count == sequence_count - 1:
+        return "all_but_one"
+    if carrier_count == sequence_count:
+        return "all_taxa"
+    return "branch_shared"
+
+
+def locus_pattern(
+    event_modes: list[str],
+    event_count: int,
+    max_unique_events: int,
+) -> str:
+    modes = set(event_modes)
+    if max_unique_events >= 2:
+        return "single_species_scattered_multiple_events"
+    if event_count == 1 and modes == {"single_species"}:
+        return "single_species_single_event"
+    if modes == {"all_but_one"}:
+        return "all_but_one"
+    if modes == {"branch_shared"}:
+        return "branch_shared"
+    if modes == {"all_taxa"}:
+        return "all_taxa"
+    if not modes:
+        return "no_frameshift_detected"
+    return "mixed"
+
+
+def merge_overlapping_triplets(
+    triplets: list[tuple[tuple[int, int, int], str]],
+) -> list[tuple[tuple[int, ...], list[str]]]:
+    components: list[tuple[set[int], set[str]]] = []
+
+    for columns, carrier in triplets:
+        column_set = set(columns)
+        matching = [
+            index
+            for index, (component_columns, _) in enumerate(components)
+            if component_columns & column_set
+        ]
+
+        if not matching:
+            components.append((column_set, {carrier}))
+            continue
+
+        merged_columns = set(column_set)
+        merged_carriers = {carrier}
+        for index in reversed(matching):
+            component_columns, component_carriers = components.pop(index)
+            merged_columns.update(component_columns)
+            merged_carriers.update(component_carriers)
+        components.append((merged_columns, merged_carriers))
+
+    return [
+        (tuple(sorted(columns)), sorted(carriers))
+        for columns, carriers in components
+    ]
+
+
+def classify_fasta(fasta_path: Path) -> tuple[list[str], list[list[str]]]:
+    sequences = read_fasta(fasta_path)
+    sequence_ids = list(sequences)
+    sequence_count = len(sequence_ids)
+    unique_events_by_sequence: Counter[str] = Counter()
+    event_rows: list[list[str]] = []
+    modes: list[str] = []
+    complete_triplets: list[tuple[tuple[int, int, int], str]] = []
+    events: list[tuple[tuple[int, ...], list[str]]] = []
+
+    for sequence_id, sequence in sequences.items():
+        marker_columns = [
+            column_index + 1
+            for column_index, character in enumerate(sequence)
+            if character == "!"
+        ]
+        for index in range(0, len(marker_columns), 3):
+            event_columns = tuple(marker_columns[index : index + 3])
+            if len(event_columns) == 3:
+                complete_triplets.append((event_columns, sequence_id))
+            else:
+                # A trailing 1-2 markers still form one unique event, but cannot
+                # be classified as shared because they are not a complete triplet.
+                events.append((event_columns, [sequence_id]))
+
+    events.extend(merge_overlapping_triplets(complete_triplets))
+    events.sort(key=lambda event: (event[0], event[1]))
+
+    for event_columns, carriers in events:
+        carriers = [
+            sequence_id for sequence_id in sequence_ids if sequence_id in carriers
+        ]
+        mode = event_mode(len(carriers), sequence_count)
+        modes.append(mode)
+        if len(carriers) == 1:
+            unique_events_by_sequence[carriers[0]] += 1
+        event_rows.append(
+            [
+                fasta_path.stem,
+                ",".join(str(column) for column in event_columns),
+                str(len(carriers)),
+                ",".join(carriers),
+                mode,
+            ]
+        )
+
+    event_count = len(event_rows)
+    max_unique_events = max(unique_events_by_sequence.values(), default=0)
+    severe = event_count >= 3 or max_unique_events >= 2
+    locus_row = [
+        fasta_path.stem,
+        str(sequence_count),
+        str(event_count),
+        str(max_unique_events),
+        locus_pattern(modes, event_count, max_unique_events),
+        "yes" if severe else "no",
+    ]
+    return locus_row, event_rows
+
+
+def fasta_paths_from_list(list_path: Path) -> list[Path]:
+    paths: list[Path] = []
+    seen: set[Path] = set()
+    with list_path.open() as fasta_list:
+        for line_number, raw_line in enumerate(fasta_list, start=1):
+            value = raw_line.strip()
+            if not value:
+                continue
+            path = Path(value)
+            if not path.is_absolute():
+                path = list_path.parent / path
+            path = path.resolve()
+            if not path.is_file():
+                raise FileNotFoundError(
+                    f"{list_path}:{line_number}: FASTA file not found: {path}"
+                )
+            if path not in seen:
+                paths.append(path)
+                seen.add(path)
+    return paths
+
+
+def write_results(
+    output_dir: Path,
+    locus_rows: list[list[str]],
+    event_rows: list[list[str]],
+) -> None:
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    with (output_dir / "frameshift_events.tsv").open("w", newline="") as output:
+        writer = csv.writer(output, delimiter="\t", lineterminator="\n")
+        writer.writerow(["OG", "alignment_column", "carrier_count", "carriers", "mode"])
+        writer.writerows(event_rows)
+
+    with (output_dir / "frameshift_locus_classification.tsv").open(
+        "w", newline=""
+    ) as output:
+        writer = csv.writer(output, delimiter="\t", lineterminator="\n")
+        writer.writerow(LOCUS_HEADER)
+        writer.writerows(locus_rows)
+
+    severe_rows = [row for row in locus_rows if row[-1] == "yes"]
+    with (output_dir / "frameshift_severe_subset.tsv").open(
+        "w", newline=""
+    ) as output:
+        writer = csv.writer(output, delimiter="\t", lineterminator="\n")
+        writer.writerow(LOCUS_HEADER)
+        writer.writerows(severe_rows)
+
+    pattern_counts = Counter(row[4] for row in locus_rows)
+    with (output_dir / "frameshift_classification_summary.tsv").open(
+        "w", newline=""
+    ) as output:
+        writer = csv.writer(output, delimiter="\t", lineterminator="\n")
+        writer.writerow(["metric", "count"])
+        writer.writerow(["total_FS_loci", len(locus_rows)])
+        writer.writerow(["severe_loci", len(severe_rows)])
+        for pattern, count in sorted(pattern_counts.items()):
+            writer.writerow([pattern, count])
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description="Classify frameshift loci by event count and carrier pattern."
+    )
+    parser.add_argument(
+        "--fasta-list",
+        required=True,
+        type=Path,
+        help="File containing one aligned FASTA path per line.",
+    )
+    parser.add_argument(
+        "--output-dir",
+        required=True,
+        type=Path,
+        help="Directory for classification tables.",
+    )
+    args = parser.parse_args()
+
+    fasta_paths = fasta_paths_from_list(args.fasta_list.resolve())
+    locus_rows: list[list[str]] = []
+    event_rows: list[list[str]] = []
+    for fasta_path in fasta_paths:
+        locus_row, rows = classify_fasta(fasta_path)
+        locus_rows.append(locus_row)
+        event_rows.extend(rows)
+    write_results(args.output_dir, locus_rows, event_rows)
+
+
+if __name__ == "__main__":
+    main()

99.scripts_v2/utils/get_og_seqs.py

@@ -0,0 +1,148 @@
+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Extract Single Copy OG sequences from a comprehensive FASTA file based on OG definitions from results of orthofinder.
+"""
+
+import sys
+from Bio import SeqIO
+from Bio.SeqRecord import SeqRecord
+from Bio.Seq import Seq
+from pathlib import Path
+import argparse
+
+
+def get_og_fastas(og_dir: Path, ext: str) -> dict[str, list[str]]:
+    """
+    Get a dictionary of OGs and their corresponding sequence IDs from FASTA files in a directory.
+
+    Args:
+        og_dir (Path): Path to the directory containing OG FASTA files.
+        ext (str): Extension of the FASTA files. Default is ".fa".
+
+    Returns:
+        dict: Dictionary with OG names as keys and lists of sequence IDs as values.
+    """
+    og_dict: dict[str, list[str]] = {}
+    try:
+        for fasta_file in og_dir.glob(f"*{ext}"):
+            og_name = fasta_file.stem
+            seq_ids: list[str] = []
+            for record in SeqIO.parse(fasta_file, "fasta"):
+                seq_ids.append(record.id)
+            og_dict[og_name] = seq_ids
+    except Exception as e:
+        print(f"Error processing OG FASTA files in {og_dir}: {e}")
+        sys.exit(1)
+    return og_dict
+
+
+def parse_all_fasta(fasta_file: Path) -> dict[str, Seq]:
+    """
+    Parse a FASTA file and return a dictionary of sequence IDs and their corresponding Seq objects.
+
+    Args:
+        fasta_file (Path): Path to the FASTA file.
+    Returns:
+        dict: Dictionary with sequence IDs as keys and Seq objects as values.
+    """
+    seq_dict: dict[str, Seq] = {}
+    try:
+        for record in SeqIO.parse(fasta_file, "fasta"):
+            seq_dict[record.id] = record.seq
+    except Exception as e:
+        print(f"Error parsing FASTA file {fasta_file}: {e}")
+        sys.exit(1)
+    return seq_dict
+
+
+def output_og_seqs(
+    all_seq_dict: dict[str, Seq], og_dict: dict[str, list[str]], output_dir: Path
+):
+    """
+    Output sequences for each OG into separate FASTA files, remove gene id and only keep taxon name.
+
+    Args:
+        all_seq_dict (dict): Dictionary of all sequences with sequence IDs as keys.
+        og_dict (dict): Dictionary of OGs with OG names as keys and lists of sequence IDs as values.
+        output_dir (Path): Path to the output directory.
+    """
+    for og_name, seq_ids in og_dict.items():
+        og_seqs: list[SeqRecord] = []
+        for seq_id in seq_ids:
+            if seq_id in all_seq_dict:
+                seq_record = SeqRecord(
+                    all_seq_dict[seq_id], id=seq_id.split("@")[0], description=""
+                )
+                og_seqs.append(seq_record)
+            else:
+                print(f"Warning: Sequence ID {seq_id} not found in all sequences.")
+        output_fasta = output_dir / f"{og_name}.fa"
+        try:
+            SeqIO.write(og_seqs, output_fasta, "fasta")
+        except Exception as e:
+            print(f"Error writing to FASTA file {output_fasta}: {e}")
+            sys.exit(1)
+
+
+def write_og_list(og_dict: dict[str, list[str]]):
+    """
+    Write the OG list to a text file.
+
+    Args:
+        og_dict (dict): Dictionary of OGs with OG names as keys and lists of sequence IDs as values.
+    """
+    try:
+        with open("og_list.tsv", "w") as f:
+            for og_name, seq_ids in og_dict.items():
+                line = f"{og_name}\t" + "\t".join(seq_ids) + "\n"
+                f.write(line)
+        print("OG list written to og_list.tsv")
+    except Exception as e:
+        print(f"Error writing OG list to file: {e}")
+        sys.exit(1)
+
+
+def main(og_dir_path: str, all_fasta_path: str, output_dir_path: str, ext: str = ".fa"):
+    og_dir = Path(og_dir_path)
+    all_fasta = Path(all_fasta_path)
+    output_dir = Path(output_dir_path)
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    og_dict = get_og_fastas(og_dir, ext)
+    all_seqs = parse_all_fasta(all_fasta)
+    output_og_seqs(all_seqs, og_dict, output_dir)
+    write_og_list(og_dict)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="Extract OG sequences from a comprehensive FASTA file based on OG definitions."
+    )
+    parser.add_argument(
+        "-d",
+        "--og_dir",
+        required=True,
+        help="Directory containing OG FASTA files.",
+    )
+    parser.add_argument(
+        "-a",
+        "--all_fasta",
+        required=True,
+        help="FASTA file containing all sequences.",
+    )
+    parser.add_argument(
+        "-o",
+        "--output_dir",
+        required=True,
+        help="Output directory for OG FASTA files.",
+    )
+    parser.add_argument(
+        "-e",
+        "--ext",
+        default=".fa",
+        help="Extension of OG FASTA files (default: .fa).",
+    )
+    args = parser.parse_args()
+    main(args.og_dir, args.all_fasta, args.output_dir, args.ext)

99.scripts_v2/utils/mbblock.txt

@@ -0,0 +1,12 @@
+begin mrbayes;
+set nowarnings=yes;
+set usebeagle=yes;
+set autoclose=yes;
+set seed=12345;
+set swapseed=12345;
+lset nst=6 rates=gamma;
+mcmcp ngen=2000000 burninfrac=.25 samplefreq=1000 printfreq=10000 checkpoint=no
+diagnfreq=10000 nruns=3 nchains=3 temp=0.40 swapfreq=10 stoprule=no;
+mcmc;
+sumt;
+end;

99.scripts_v2/utils/mrbayes.sh

@@ -0,0 +1,20 @@
+#! /bin/bash
+set -e
+
+dirname=$(dirname "$0")
+if [ "$#" -ne 2 ]; then
+    echo "Usage: $0 <input_fasta> <output_nexus>"
+    exit 1
+fi
+
+in_fasta=$(readlink -f "$1")
+out_nex=$2
+
+# Convert fasta to nexus
+out_dir=$(dirname "$out_nex")
+mkdir -p "$out_dir"
+seqmagick convert --output-format nexus --alphabet dna --input-format fasta "$in_fasta" "$out_nex"
+# append MrBayes block to nexus file
+cat "$dirname/mbblock.txt" >>"$out_nex"
+# Run MrBayes
+mb "$out_nex" >"${out_dir}"/mrbayes.log 2>&1

99.scripts_v2/utils/raxml.sh

@@ -0,0 +1,20 @@
+#! /bin/bash
+set -e
+if [ "$#" -ne 4 ]; then
+    echo "Usage: $0 <alignment.fa> <output_prefix> <threads> <outgroup>"
+    echo "Run modeltest-ng and raxml-ng on a given alignment"
+    exit 1
+fi
+
+aln=$(readlink -f "$1")
+out=$2
+threads=$3
+outgroup=$4
+outdir=$(dirname "$out")
+mkdir -p "$outdir"
+# run modeltest-ng
+modeltest-ng -p "$threads" -r 12345 --force -i "$aln" -d nt -t ml -o "$out".modeltest
+# run raxml-ng
+cmd=$(grep "raxml-ng" "$out".modeltest.out | tail -n 1 | sed 's/> //')
+params="--blopt nr_safe --all --bs-trees 1000 --outgroup $outgroup --redo --threads $threads --seed 12345 --prefix $out"
+bash -c "$cmd $params" >/dev/null

99.scripts_v2/utils/snaq.jl

@@ -0,0 +1,195 @@
+#! /usr/bin/env julia
+# Installing Dependencies
+## using Pkg
+## Pkg.add("Distributed")
+## Pkg.add("DataFrames")
+## Pkg.add("CSV")
+## Pkg.add("SNaQ")
+## Pkg.add("PhyloNetworks")
+## Pkg.add("RCall")
+## Pkg.add("PhyloPlots")
+## Pkg.add("QuartetNetworkGoodnessFit")
+
+# Running SNaQ Analysis
+using PhyloNetworks, SNaQ;
+using Distributed;
+addprocs(9);
+@everywhere using PhyloNetworks, SNaQ;
+nruns = 100; # number of runs for each hmax
+astralfile = joinpath("..", "species_tree", "coal.tre");
+astraltree = readnewick(astralfile);
+
+### Reading RAxML gene trees and ASTRAL species tree
+### running in raxml_snaq/ folder
+# raxmltrees = joinpath("..", "..", "species_tree", "all.trees");
+# inputCF = readtrees2CF(raxmltrees);
+# net0 = snaq!(astraltree, inputCF, hmax=0, filename="net0", seed=123, outgroup="ZJ", runs=nruns);
+# net1 = snaq!(net0, inputCF, hmax=1, filename="net1", seed=123, outgroup="ZJ", runs=nruns);
+# net2 = snaq!(net1, inputCF, hmax=2, filename="net2", seed=123, outgroup="ZJ", runs=nruns);
+# net3 = snaq!(net2, inputCF, hmax=3, filename="net3", seed=123, outgroup="ZJ", runs=nruns);
+# net4 = snaq!(net3, inputCF, hmax=4, filename="net4", seed=123, outgroup="ZJ", runs=nruns);
+
+### Alternatively, reading in the input files from Bucky
+### running in input_snaq/ folder
+inputCFfile = joinpath("..", "bucky", "bucky_1", "input.mb.CFs.csv");
+inputCF = readtableCF(inputCFfile);
+net0 = snaq!(astraltree, inputCF, hmax=0, filename="net0", seed=123, outgroup="ZJ", runs=nruns);
+net1 = snaq!(net0, inputCF, hmax=1, filename="net1", seed=123, outgroup="ZJ", runs=nruns);
+net2 = snaq!(net1, inputCF, hmax=2, filename="net2", seed=123, outgroup="ZJ", runs=nruns);
+net3 = snaq!(net2, inputCF, hmax=3, filename="net3", seed=123, outgroup="ZJ", runs=nruns);
+net4 = snaq!(net3, inputCF, hmax=4, filename="net4", seed=12345, outgroup="ZJ", runs=nruns);
+
+
+# Plotting the SNaQ results
+using PhyloPlots, RCall;
+using PhyloNetworks, SNaQ;
+## Network scores vs. hmax
+net0 = readsnaqnetwork("net0.out");
+net1 = readsnaqnetwork("net1.out");
+net2 = readsnaqnetwork("net2.out");
+net3 = readsnaqnetwork("net3.out");
+net4 = readsnaqnetwork("net4.out");
+scores = [loglik(net0), loglik(net1), loglik(net2), loglik(net3), loglik(net4)];
+hmax = collect(0:4);
+R"pdf"("snaq_network_scores.pdf", width=12, height=8);
+R"plot"(hmax, scores, type="b", ylab="network score", xlab="hmax", col="blue");
+R"dev.off"();
+
+## Rerooting and rotating the networks for better visualization
+rootatnode!(net1, "ZJ");
+rootatnode!(net2, "ZJ");
+rootatnode!(net3, "ZJ");
+rootatnode!(net4, "ZJ");
+### rotate!(net1, -2);
+### rotate!(net2, -2);
+### rotate!(net3, -2);
+### rotate!(net4, -2);
+
+## Plotting the networks
+R"pdf"("snaq_networks_fulltree.pdf", width=14, height=10);
+R"layout(matrix(1:4, 2, 2, byrow=TRUE))"; # to get 4 plots into a single figure: 2 row, 2 columns
+R"par"(mar=[0, 0, 1.5, 0]);        # for smaller margins
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net1, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net1, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=1, loglik=-", round(loglik(net1), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net2, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net2, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=2, loglik=-", round(loglik(net2), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net3, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net3, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=3, loglik=-", round(loglik(net3), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net4, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net4, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=4, loglik=-", round(loglik(net4), digits=2)), font=2);
+R"dev.off"();
+
+R"pdf"("snaq_networks_majortree.pdf", width=14, height=10);
+R"layout(matrix(1:4, 2, 2, byrow=TRUE))"; # to get 4 plots into a single figure: 2 row, 2 columns
+R"par"(mar=[0, 0, 1.5, 0]);        # for smaller margins
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net1, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net1, style=:majortree, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=1, loglik=-", round(loglik(net1), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net2, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net2, style=:majortree, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=2, loglik=-", round(loglik(net2), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net3, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net3, style=:majortree, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=3, loglik=-", round(loglik(net3), digits=2)), font=2);
+xmin, xmax = PhyloPlots.PhyloPlots.edgenode_coordinates(net4, false, false)[13:14];
+xmax += (xmax - xmin) * 0.3;
+plot(net4, style=:majortree, showgamma=true, tipoffset=0.1, xlim=[xmin, xmax]);
+R"mtext"(string("hmax=4, loglik=-", round(loglik(net4), digits=2)), font=2);
+R"dev.off"();
+
+## expected vs. observed quartet concordance factors
+using CSV, DataFrames, Distributions, Random, RCall;
+using PhyloNetworks, SNaQ;
+using QuartetNetworkGoodnessFit;
+using Statistics, Printf;
+
+inputCFfile = joinpath("..", "bucky", "bucky_1", "input.mb.CFs.csv");
+inputCF = readtableCF(inputCFfile);
+Random.seed!(123);
+jitter = 0.005;
+hmax_list = 0:3;
+
+R"pdf"("expCFs_obsvsfitted.pdf", width=6, height=6);
+R"par(mfrow=c(2,2), mar=c(4,4,3,1))"
+gof_results = DataFrame(h=Int[], pvalue_normal=Float64[], pvalue_empirical=Float64[], z_uncorrected=Float64[], sigma=Float64[], nsim=Int[])
+
+for h in hmax_list
+    net = readsnaqnetwork("net$(h).out")
+    topologymaxQpseudolik!(net, inputCF)
+
+    df_long = fittedquartetCF(inputCF, :long)
+    df_jittered = deepcopy(df_long)
+    df_jittered.obsCF .+= rand(Uniform(-jitter/2, jitter/2), nrow(df_long))
+    df_jittered.expCF .+= rand(Uniform(-jitter/2, jitter/2), nrow(df_long))
+
+    res = quarnetGoFtest!(net, inputCF, false; seed=123, nsim=1000)
+    z_obs = res[2]
+    z_sim = sort!(res[6])
+    p_emp = mean(z_sim .>= z_obs)      # 经验 p 值：模拟 z 超过观测 z 的比例
+
+    # 存表时两种 p 值都留，加上模拟分布的特征备查
+    push!(gof_results, (h=h,
+        pvalue_normal=res[1],
+        pvalue_empirical=p_emp,
+        z_uncorrected=z_obs,
+        sigma=res[3],
+        nsim=1000))
+
+
+    ptxt = p_emp == 0 ? "p < 0.001 (emp.)" : @sprintf("p = %.3f (emp.)", p_emp)
+    ttl = "h = $h, $ptxt"
+    R"plot"(df_jittered.obsCF, df_jittered.expCF,
+        xlim=[0, 1.05], ylim=[0, 1.05],
+        xlab="quartet CF observed in gene trees",
+        ylab="quartet CF expected from network",
+        main=ttl,
+        pch=16, col="#008080", cex=0.5)
+    R"abline"(0, 1, col="#008080", lwd=0.5)
+end
+
+R"dev.off"();
+
+CSV.write("gof_results.csv", gof_results)
+println(gof_results)
+
+
+net_final = readsnaqnetwork("net0.out");
+topologymaxQpseudolik!(net_final, inputCF);
+df_long = fittedquartetCF(inputCF, :long);
+
+### Adding jitter to the points for better visualization
+df_jittered = deepcopy(df_long);
+df_jittered.obsCF .+= rand(Uniform(-jitter / 2, jitter / 2), nrow(df_long));
+df_jittered.expCF .+= rand(Uniform(-jitter / 2, jitter / 2), nrow(df_long));
+df_wide = fittedquartetCF(inputCF, :wide)
+println(names(df_wide))
+df_wide.outlier_p = res1[4]
+outliers = sort(df_wide[df_wide.outlier_p .< 0.01, :], :outlier_p)
+CSV.write("outlier_quartets_h1.csv", outliers)
+
+
+R"plot"(df_jittered.obsCF, df_jittered.expCF, xlab="quartet CF observed in gene trees",
+    ylab="quartet CF expected from network", pch=16,
+    col="#008080",
+    cex=0.7);
+R"abline"(0, 1, col="#008080", lwd=0.5);
+R"dev.off"();
+
+# Goodness of fit of the SNaQ networks
+using QuartetNetworkGoodnessFit;
+res1 = quarnetGoFtest!(net_final, inputCF, true; seed=123, nsim=1000);
+res1[[1, 2, 3]] # p-value, uncorrected z, σ
+z_obs = res1[2]
+z_sim = sort!(res1[6])
+p_emp = mean(z_sim .>= z_obs)

99.scripts_v2/utils/snaq_raw.jl

@@ -0,0 +1,19 @@
+#! /usr/bin/env julia
+
+
+# Running SNaQ Analysis
+using PhyloNetworks, SNaQ;
+using Distributed;
+addprocs(9);
+@everywhere using PhyloNetworks, SNaQ;
+nruns = 100; # number of runs for each hmax
+astralfile = joinpath("..", "species_tree", "coal.tre");
+astraltree = readnewick(astralfile);
+
+inputCFfile = joinpath("..", "bucky", "bucky_1", "input.mb.CFs.csv");
+inputCF = readtableCF(inputCFfile);
+net0 = snaq!(astraltree, inputCF, hmax=0, filename="net0", seed=123, outgroup="ZJ", runs=nruns);
+net1 = snaq!(net0, inputCF, hmax=1, filename="net1", seed=123, outgroup="ZJ", runs=nruns);
+net2 = snaq!(net1, inputCF, hmax=2, filename="net2", seed=123, outgroup="ZJ", runs=nruns);
+net3 = snaq!(net2, inputCF, hmax=3, filename="net3", seed=123, outgroup="ZJ", runs=nruns);
+net4 = snaq!(net3, inputCF, hmax=4, filename="net4", seed=12345, outgroup="ZJ", runs=nruns);