20161002

.envrc

@@ -1,6 +1,6 @@
 set +u
 export PROJECTHOME="/home/ywtang/project/biyelunwen"
-export SCRIPTS="$PROJECTHOME/99.scripts"
+export SCRIPTS="$PROJECTHOME/99.scripts_v2"
 export PUEUE_CONFIG_PATH="$PROJECTHOME/.pueue.yml"
 watch_file pixi.lock
 eval "$(pixi shell-hook)"

.gitignore

@@ -1,6 +1,7 @@
 # pixi environments
 .pixi/*
 !.pixi/config.toml
+**/__pycache__
 .pueue/*
 00.rawdata/*
 01.reference/*
@@ -18,9 +19,12 @@
 11.reference_v2/*
 12.orthofinder/*
 98.results/*
-99.scripts/bucky/
-99.scripts/phyparts/
-99.scripts/bpp/
+99.scripts*/bucky/
+99.scripts*/phyparts/
+99.scripts*/bpp/
 run.status.json
 !*/_description.md
 tmp/
+13.gene_phylogeny/*
+14.taxa_phylogeny/*
+15.dsuite/*

.vscode/settings.json

@@ -1,5 +1,6 @@
 {
     "julia.executablePath": "/home/ywtang/project/biyelunwen/.pixi/envs/default/bin/julia",
     "julia.enableCrashReporter": false,
-    "julia.enableTelemetry": false
+    "julia.enableTelemetry": false,
+    "python-envs.defaultEnvManager": "ms-python.python:system"
 }

99.scripts/workflow/phylogeny_reconstruction/06.snaq.jl

@@ -1,5 +1,5 @@
 #! /usr/bin/env julia
-## Installing Dependencies
+# Installing Dependencies
 ## using Pkg
 ## Pkg.add("Distributed")
 ## Pkg.add("DataFrames")

99.scripts_v2/01.prepare_sequences.sh

@@ -11,21 +11,22 @@ seqkit replace -k HG.map.txt -p "^(\S+).*$" -r "{kv}" GCA_030763125.1_ASM3076312
 # MG & ZG use hapA genome
 grep -Fwf MG.target.txt MG.gff > MG.target.gff
 grep -Fwf ZG.target.txt ZG.gff > ZG.target.gff
+cd ..
 
 # Extract CDS sequences from GFF3 files using gffread
-pueue add -- gffread -g raw/HN.fa -x HN.cds.fa -y HN.pep.fa raw/HN.gff3
-pueue add -- gffread -g raw/HG.fa -x HG.cds.fa -y HG.pep.fa raw/Hippophae_gyantsensis.gff
-pueue add -- gffread -g raw/HS.fa -x HS.cds.fa -y HS.pep.fa raw/HS.gff3
-pueue add -- gffread -g raw/YN.genome.fa -x YN.cds.fa -y YN.pep.fa raw/YN.gff3
-pueue add -- gffread -g raw/ZG.genome.fa -x ZG.cds.fa -y ZG.pep.fa raw/ZG.target.gff
-pueue add -- gffread -g raw/ZJ.genome.fa -x ZJ.cds.fa -y ZJ.pep.fa raw/ZJ.gff
-pueue add -- gffread -g raw/ZY.fa -x ZY.cds.fa -y ZY.pep.fa raw/ZY.gff3
-pueue add -- gffread -g raw/Hippophae_tibetana.genomic.fasta -x HT.cds.fa -y HT.pep.fa raw/Hippophae_tibetana.mRNA.gff
-pueue add -- gffread -g raw/MG.genome.fasta -x MG.cds.fa -y MG.pep.fa raw/MG.target.gff
+gffread -V -g raw/HN.fa -x HN.cds.fa -y HN.pep.fa raw/HN.gff3
+gffread -V -g raw/HG.fa -x HG.cds.fa -y HG.pep.fa raw/Hippophae_gyantsensis.gff
+gffread -V -g raw/HS.fa -x HS.cds.fa -y HS.pep.fa raw/HS.gff3
+gffread -V -g raw/YN.genome.fa -x YN.cds.fa -y YN.pep.fa raw/YN.gff3
+gffread -V -g raw/ZG.genome.fa -x ZG.cds.fa -y ZG.pep.fa raw/ZG.target.gff
+gffread -V -g raw/ZJ.genome.fa -x ZJ.cds.fa -y ZJ.pep.fa raw/ZJ.gff
+gffread -V -g raw/ZY.fa -x ZY.cds.fa -y ZY.pep.fa raw/ZY.gff3
+gffread -V -g raw/Hippophae_tibetana.genomic.fasta -x HT.cds.fa -y HT.pep.fa raw/Hippophae_tibetana.mRNA.gff
+gffread -V -g raw/MG.genome.fasta -x MG.cds.fa -y MG.pep.fa raw/MG.target.gff
 cp raw/EA.gff3.cds EA.cds.fa
 cp raw/EM.genome.cds.fasta EM.cds.fa
-cp raw/EA.gff3.pep EA.pep.fa
-cp raw/EM.pep.fasta EM.pep.fa
+seqkit grep -s -v -r -p '[.*][A-Z]' raw/EA.gff3.pep > EA.pep.fa
+seqkit grep -s -v -r -p '[.*][A-Z]' raw/EM.pep.fasta > EM.pep.fa
 
 mkdir pep
 mkdir cds
@@ -33,8 +34,15 @@ mv ./*.cds.fa cds
 mv ./*.pep.fa pep
 
 cd pep
+# 保留主要转录本
 primary_transcript ./ last_dot
-
 # EM基因命名特殊性，导致直接使用primary_transcript会出现问题，直接复制即可
 cp EM.pep.fa primary_transcripts/EM.pep.fa
-
+# 移除终止密码子和序列描述
+cd primary_transcripts
+for i in *.pep.fa; do
+    base=$(basename "$i" .pep.fa)
+    echo "Processing $base"
+    sed -i "s/\.$//g" "$i"
+    sed -i "/^>/ s/[[:space:]].*//" "$i"
+done

99.scripts_v2/02.orthofinder.sh

@@ -1,3 +1,7 @@
+#!/bin/bash
+
+set -euo pipefail
+
 # prepare input files for orthofinder
 mkdir input_pep
 cp ../11.reference_v2/pep/primary_transcripts/*.pep.fa input_pep
@@ -17,5 +21,18 @@ done
 pueue add -- orthofinder -t 12 -a 4 -f input_pep/core -o orthofinder_out
 
 # assign outgroups
-pueue add -- orthofinder -t 12 -a 4 --assign input_pep/outgroup --core orthofinder_out/Results_Hippophae_only
+pueue add -- orthofinder -t 12 -a 4 --assign input_pep/outgroup --core orthofinder_out/Results_Aug03
 
+ln -s orthofinder_out/Results_Aug03 Hippophae_only
+ln -s orthofinder_out/Results_Aug03_1/ Hippophae_outgroup
+
+# get all cds into a file
+echo -ne >all_cds.fa
+for i in ../11.reference_v2/cds/*.cds.fa; do
+    base=$(basename "$i" .cds.fa)
+    echo "Processing $base"
+    sed "s/>/>${base}@/g" "$i" >>all_cds.fa
+done
+
+# get og list and extract cds sequences for each og
+pueue add -- python3 "$PROJECTHOME"/99.scripts_v2/utils/get_og_seqs.py -d Hippophae_outgroup/Single_Copy_Orthologue_Sequences -a all_cds.fa -o ogs_cds

99.scripts_v2/03.gene_phylogeny.sh

@@ -0,0 +1,290 @@
+#!/bin/bash
+
+set -euo pipefail
+
+ln -s ../12.orthofinder/ogs_cds/ .
+
+# parallel MACSE alignment for each orthologous group
+mkdir alignments
+echo -n >macse.cmds
+for og_fasta in ogs_cds/*.fa; do
+    og_name=$(basename "$og_fasta" .fa)
+    out_stem="alignments/$og_name"
+    echo "macse -prog alignSequences -seq $og_fasta -out_AA ${out_stem}.pal -out_NT ${out_stem}.nal > ${out_stem}.log 2>&1" >>macse.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a macse.cmds bash -c "cmd"
+
+# "!" --> "-"
+rg -l "!" alignments/*.nal >fs.list
+mkdir alignments_clean
+for og_fasta in ogs_cds/*.fa; do
+    og_name=$(basename "$og_fasta" .fa)
+    out_stem="alignments/$og_name"
+    echo "Processing $og_name"
+    sed 's/!/-/g' "${out_stem}.nal" >"alignments_clean/${og_name}.nal"
+done
+
+# classify frameshift loci by independent event count and carrier pattern
+python3 "$PROJECTHOME"/99.scripts_v2/utils/classify_frameshifts.py \
+    --fasta-list fs.list \
+    --output-dir frameshift_summary
+
+# trim alignments with trimAl
+mkdir -p alignments_trimmed
+echo -n >trimal.cmds
+for i in alignments_clean/*.nal; do
+    j=$(basename "$i" .nal)
+    out_dir="alignments_trimmed"
+    echo "trimal -in $i -out $out_dir/${j}.trimed.fa -automated1" >>trimal.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a trimal.cmds bash -c "cmd"
+
+# retain MSAs at least 300 bp long
+min_alignment_length=300
+mkdir -p alignments_filtered
+printf "gene\talignment_length_bp\n" >alignments_shorter_than_300bp.tsv
+printf "gene\tremoved_sequence\n" >trimal_removed_sequences.tsv
+for i in alignments_trimmed/*.fa; do
+    [[ -e "$i" ]] || continue
+    j=$(basename "$i" .trimed.fa)
+    untrimmed_alignment="alignments_clean/${j}.nal"
+    filtered_alignment="alignments_filtered/${j}.fa"
+    if [[ ! -s "$untrimmed_alignment" ]]; then
+        printf "Missing non-empty pre-trim alignment: %s\n" "$untrimmed_alignment" >&2
+        exit 1
+    fi
+
+    removed_sequences=()
+    mapfile -t removed_sequences < <(
+        comm -23 \
+            <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' \
+                "$untrimmed_alignment" | sort -u) \
+            <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' \
+                "$i" | sort -u)
+    )
+    if ((${#removed_sequences[@]} > 0)); then
+        rm -f "$filtered_alignment"
+        for sequence in "${removed_sequences[@]}"; do
+            printf "%s\t%s\n" "$j" "$sequence" >>trimal_removed_sequences.tsv
+        done
+        continue
+    fi
+
+    alignment_length=$(awk '
+        /^>/ {
+            if (found_sequence)
+                exit
+            found_sequence = 1
+            next
+        }
+        found_sequence {
+            gsub(/[[:space:]]/, "")
+            sequence_length += length($0)
+        }
+        END { print sequence_length + 0 }
+    ' "$i")
+
+    if ((alignment_length >= min_alignment_length)); then
+        if [[ ! "$i" -ef "$filtered_alignment" ]]; then
+            ln -f "$i" "$filtered_alignment"
+        fi
+    else
+        rm -f "$filtered_alignment"
+        printf "%s\t%d\n" "$j" "$alignment_length" >>alignments_shorter_than_300bp.tsv
+    fi
+done
+
+# raxml-ng for each filtered orthologous group
+mkdir -p raxml
+echo -n >raxml.cmds
+for i in alignments_filtered/*.fa; do
+    [[ -e "$i" ]] || continue
+    j=$(basename "$i" .fa)
+    out_prefix="raxml/${j}/${j}"
+    mkdir -p "raxml/${j}"
+    echo "$SCRIPTS/utils/raxml.sh $i $out_prefix 1 ZJ" >>raxml.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a raxml.cmds bash -c "cmd"
+
+# mrbayes for each filtered orthologous group
+mkdir -p mrbayes
+echo -n >mrbayes.cmds
+for i in alignments_filtered/*.fa; do
+    [[ -e "$i" ]] || continue
+    j=$(basename "$i" .fa)
+    out_prefix="mrbayes/${j}/${j}"
+    mkdir -p "mrbayes/${j}"
+    echo "bash $SCRIPTS/utils/mrbayes.sh $i ${out_prefix}.nex" >>mrbayes.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a mrbayes.cmds bash -c "cmd"
+
+# mbsum
+mkdir -p mbsum
+echo -n "" >mbsum.skip.log
+echo -n "" >mbsum.cmds
+for i in mrbayes/*/mrbayes.log; do
+    og_dir=$(dirname "$i")
+    base=$(basename "$og_dir")
+    ## skip Average standard deviation of split frequencies > 0.01
+    dsf=$(awk '/Average standard deviation of split frequencies:/ {out=$7} END{print out+0}' "$i" 2>/dev/null)
+    if awk -v d="$dsf" -v min_dsf=0.01 'BEGIN{if (d >= min_dsf) exit 0; exit 1}'; then
+        echo "Skipping ${base} due to high DSF: ${dsf}" >>mbsum.skip.log
+        continue
+    else
+        echo "mbsum ${og_dir}/*.t -n 1000 -o mbsum/${base}.in >> $og_dir/mbsum.log 2>&1" >>mbsum.cmds
+    fi
+done
+export PATH=$PATH:$SCRIPTS/bucky/bin
+pueue add -- xargs -t -P 12 -I cmd -a mbsum.cmds bash -c "cmd"
+
+#########################################################################
+#                                                                       #
+#  DEPRECATED SECTION: TreeShrink and FastTree-based filtering          #
+#                                                                       #
+#########################################################################
+
+# remove long branches with TreeShrink and FastTree
+mkdir -p treeshrink
+echo -n >fasttree.cmds
+for i in alignments_trimmed/*.trimed.fa; do
+    j=$(basename "$i" .trimed.fa)
+    mkdir -p treeshrink/"${j}"
+    absolute_path=$(readlink -f "$i")
+    ln -s "$absolute_path" "treeshrink"/"${j}"/input.fasta
+    echo "FastTree -nt -gtr -quiet treeshrink/${j}/input.fasta > treeshrink/${j}/input.tree" >>fasttree.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a fasttree.cmds bash -c "cmd"
+
+pueue add -- run_treeshrink.py -f -i treeshrink/ -t input.tree -a input.fasta
+
+# filter orthologous groups after TreeShrink based on minimum taxon number and minimum sequence length
+min_taxon=11
+min_seq_length=300
+mkdir -p final_alignments
+echo -n "" >treeshrink_filtered.list
+for i in treeshrink/*/output.fasta; do
+    j=$(dirname "$i")
+    j=$(basename "$j")
+    seqlen=$(seqkit fx2tab -C ATCG "$i" | awk '{print $3}' | sort -n | head -n 1)
+    seqnum=$(grep -c ">" "$i")
+    if [[ $seqnum -eq $min_taxon && $seqlen -ge $min_seq_length ]]; then
+        #cp -l "$i" final_alignments/"${j}.fa"
+        echo -n ""
+    else
+        echo "$i" >>treeshrink_filtered.list
+    fi
+done
+
+# summarize sequences removed by TreeShrink
+mkdir -p treeshrink_summary
+removed_sequences="treeshrink_summary/treeshrink_removed_sequences.tsv"
+removed_by_species="treeshrink_summary/treeshrink_removed_by_species.tsv"
+removed_group_summary="treeshrink_summary/treeshrink_removed_group_summary.tsv"
+printf "OG\tspecies\n" >"$removed_sequences"
+for output_fasta in treeshrink/*/output.fasta; do
+    og_name=$(basename "$(dirname "$output_fasta")")
+    input_fasta="$(dirname "$output_fasta")/input.fasta"
+    comm -23 \
+        <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' "$input_fasta" | sort -u) \
+        <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' "$output_fasta" | sort -u) |
+        while IFS= read -r species; do
+            printf "%s\t%s\n" "$og_name" "$species"
+        done >>"$removed_sequences"
+done
+printf "species\tfiltered_OG_count\n" >"$removed_by_species"
+awk -F '\t' 'NR > 1 {count[$2]++} END {for (species in count) print species, count[species]}' \
+    OFS='\t' "$removed_sequences" | sort -k2,2nr -k1,1 >>"$removed_by_species"
+printf "group\tfiltered_OG_count\n" >"$removed_group_summary"
+awk -F '\t' '
+    NR > 1 {
+        if ($2 == "EA" || $2 == "EM" || $2 == "ZJ")
+            outgroup[$1] = 1
+        else
+            ingroup[$1] = 1
+    }
+    END {
+        for (og in outgroup) outgroup_count++
+        for (og in ingroup) ingroup_count++
+        print "outgroup_EA_EM_ZJ", outgroup_count + 0
+        print "ingroup", ingroup_count + 0
+    }
+' OFS='\t' "$removed_sequences" >>"$removed_group_summary"
+
+# compare OGs containing frameshifts with OGs filtered by TreeShrink
+fs_overlap_detail="treeshrink_summary/treeshrink_fs_overlap.tsv"
+fs_overlap_summary="treeshrink_summary/treeshrink_fs_overlap_summary.tsv"
+printf "OG\tcategory\n" >"$fs_overlap_detail"
+awk -F '\t' '
+    NR == FNR {
+        og = $0
+        sub(/^.*\//, "", og)
+        sub(/\.nal$/, "", og)
+        frameshift[og] = 1
+        next
+    }
+    FNR > 1 {
+        treeshrink[$1] = 1
+    }
+    END {
+        for (og in frameshift) {
+            if (og in treeshrink)
+                print og, "overlap"
+            else
+                print og, "fs.list_only"
+        }
+        for (og in treeshrink) {
+            if (!(og in frameshift))
+                print og, "TreeShrink_only"
+        }
+    }
+' OFS='\t' fs.list "$removed_sequences" | sort -k2,2 -k1,1 >>"$fs_overlap_detail"
+
+printf "category\tOG_count\n" >"$fs_overlap_summary"
+awk -F '\t' '
+    NR > 1 {
+        count[$2]++
+    }
+    END {
+        overlap = count["overlap"] + 0
+        fs_only = count["fs.list_only"] + 0
+        treeshrink_only = count["TreeShrink_only"] + 0
+        print "fs.list_total", overlap + fs_only
+        print "TreeShrink_total", overlap + treeshrink_only
+        print "overlap", overlap
+        print "fs.list_only", fs_only
+        print "TreeShrink_only", treeshrink_only
+    }
+' OFS='\t' "$fs_overlap_detail" >>"$fs_overlap_summary"
+
+# raxml-ng for each filtered orthologous group
+mkdir -p raxml
+echo -n >raxml.cmds
+for i in final_alignments/*.fa; do
+    j=$(basename "$i" .fa)
+    out_prefix="raxml/${j}/${j}"
+    mkdir -p "raxml/${j}"
+    echo "$SCRIPTS/utils/raxml.sh $i $out_prefix 1 ZJ" >>raxml.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a raxml.cmds bash -c "cmd"
+
+# mrbayes for each filtered orthologous group
+mkdir -p mrbayes
+echo -n >mrbayes.cmds
+for i in final_alignments/*.fa; do
+    j=$(basename "$i" .fa)
+    out_prefix="mrbayes/${j}/${j}"
+    mkdir -p "mrbayes/${j}"
+    echo "bash $SCRIPTS/utils/mrbayes.sh $i ${out_prefix}.nex" >>mrbayes.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a mrbayes.cmds bash -c "cmd"
+
+## treeshrink with per-gene
+mkdir -p treeshrink_per_gene
+echo -n >treeshrink_per_gene.cmds
+for i in final_alignments/*.fa; do
+    j=$(basename "$i" .fa)
+    out_prefix="treeshrink_per_gene/${j}/${j}"
+    mkdir -p "treeshrink_per_gene/${j}"
+    echo "bash $SCRIPTS/utils/treeshrink.sh $i $out_prefix" >>treeshrink_per_gene.cmds
+done
+pueue add -- xargs -t -P 12 -I cmd -a treeshrink_per_gene.cmds bash -c "cmd"

99.scripts_v2/04.taxa_phylogeny.sh

@@ -0,0 +1,174 @@
+#!/bin/bash
+
+set -euo pipefail
+export PATH="$PATH:${SCRIPTS}/bucky/bin"
+
+# treeshrink
+mkdir -p treeshrink
+for i in ../13.gene_phylogeny/alignments_filtered/*.fa; do
+    j=$(basename "$i" .fa)
+    mkdir -p treeshrink/"${j}"
+    absolute_path=$(readlink -f "$i")
+    ln -s "$absolute_path" "treeshrink"/"${j}"/input.fasta
+    ln -s "$(readlink -f ../13.gene_phylogeny/raxml)"/"${j}/${j}".raxml.bestTree "treeshrink"/"${j}"/input.tree
+done
+pueue add -- run_treeshrink.py -f -i treeshrink/ -t input.tree -a input.fasta -m per-gene
+
+# treeshrink summary
+mkdir -p treeshrink_summary
+removed_sequences="treeshrink_summary/treeshrink_removed_sequences.tsv"
+removed_by_species="treeshrink_summary/treeshrink_removed_by_species.tsv"
+removed_group_summary="treeshrink_summary/treeshrink_removed_group_summary.tsv"
+printf "OG\tspecies\n" >"$removed_sequences"
+for output_fasta in treeshrink/*/output.fasta; do
+    og_name=$(basename "$(dirname "$output_fasta")")
+    input_fasta="$(dirname "$output_fasta")/input.fasta"
+    comm -23 \
+        <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' "$input_fasta" | sort -u) \
+        <(awk '/^>/{sub(/^>/, ""); sub(/[[:space:]].*$/, ""); print}' "$output_fasta" | sort -u) |
+        while IFS= read -r species; do
+            printf "%s\t%s\n" "$og_name" "$species"
+        done >>"$removed_sequences"
+done
+printf "species\tfiltered_OG_count\n" >"$removed_by_species"
+awk -F '\t' 'NR > 1 {count[$2]++} END {for (species in count) print species, count[species]}' \
+    OFS='\t' "$removed_sequences" | sort -k2,2nr -k1,1 >>"$removed_by_species"
+printf "group\tfiltered_OG_count\n" >"$removed_group_summary"
+awk -F '\t' '
+    NR > 1 {
+        if ($2 == "EA" || $2 == "EM" || $2 == "ZJ")
+            outgroup[$1] = 1
+        else
+            ingroup[$1] = 1
+    }
+    END {
+        for (og in outgroup) outgroup_count++
+        for (og in ingroup) ingroup_count++
+        print "outgroup_EA_EM_ZJ", outgroup_count + 0
+        print "ingroup", ingroup_count + 0
+    }
+' OFS='\t' "$removed_sequences" >>"$removed_group_summary"
+
+# LBS filter
+min_taxon=11
+echo -n "" >treeshrink_filtered.list
+for i in treeshrink/*/output.fasta; do
+    j=$(dirname "$i")
+    j=$(basename "$j")
+    seqnum=$(grep -c ">" "$i")
+    if [[ $seqnum -ge $min_taxon ]]; then
+        echo -n ""
+    else
+        echo "$j" >>treeshrink_filtered.list
+    fi
+done
+
+# data preparation
+ln -s ../13.gene_phylogeny/raxml .
+ln -s ../13.gene_phylogeny/mbsum .
+
+##
+## 第一套，所有基因保留，仅移除treeshrink过滤掉的基因
+##
+
+mkdir -p all_genes
+cd all_genes
+
+# aster coalescence tree
+for file in ../raxml/*/*.support; do
+    gene=$(basename "$(dirname "$file")")
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list; then
+        cat "$file"
+    fi
+done >all.ml.tre
+mkdir -p species_tree
+pueue add -- astral4 -t 6 -i all.ml.tre -o species_tree/coal.tre --root ZJ
+
+# bucky
+mkdir bucky
+for mbfile in ../mbsum/*.in; do
+    gene=$(basename "$mbfile" .in)
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list; then
+        cp "$mbfile" bucky/"$gene".in
+    fi
+done
+cd bucky
+tar -cf input.mb.tar -- ./*.in
+rm -- ./*.in
+pueue add -- "$SCRIPTS"/utils/bucky.pl -a 1 -o bucky_1 -s input.mb.tar
+
+# SNaQ
+mkdir snaq
+cd snaq
+## run snaq.jl with interactive mode
+# julia snaq.jl
+cd ../..
+
+##
+## 第二套，移除treeshrink过滤掉的基因和移码突变基因
+mkdir -p exclude_all_fs
+cd exclude_all_fs
+cut -f 1 ../../13.gene_phylogeny/frameshift_summary/frameshift_locus_classification.tsv | sed 1d > fs_exclude.list
+
+# aster coalescence tree
+for file in ../raxml/*/*.support; do
+    gene=$(basename "$(dirname "$file")")
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list fs_exclude.list; then
+        cat "$file"
+    fi
+done >all.ml.tre
+mkdir -p species_tree
+pueue add -- astral4 -t 6 -i all.ml.tre -o species_tree/coal.tre --root ZJ
+
+# bucky
+mkdir bucky
+for mbfile in ../mbsum/*.in; do
+    gene=$(basename "$mbfile" .in)
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list fs_exclude.list; then
+        cp "$mbfile" bucky/"$gene".in
+    fi
+done
+cd bucky
+tar -cf input.mb.tar -- ./*.in
+rm -- ./*.in
+pueue add -- "$SCRIPTS"/utils/bucky.pl -a 1 -o bucky_1 -s input.mb.tar
+
+mkdir snaq
+cd snaq
+## run snaq.jl with interactive mode
+# julia snaq.jl
+cd ../..
+
+## 第三套，移除严重移码突变基因和treeshrink过滤掉的基因
+mkdir -p exclude_severe_fs
+cd exclude_severe_fs
+grep -wv single_species_single_event ../../13.gene_phylogeny/frameshift_summary/frameshift_locus_classification.tsv | cut -f 1 | sed 1d > fs_exclude.list
+
+# aster coalescence tree
+for file in ../raxml/*/*.support; do
+    gene=$(basename "$(dirname "$file")")
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list fs_exclude.list; then
+        cat "$file"
+    fi
+done >all.ml.tre
+mkdir -p species_tree
+pueue add -- astral4 -t 6 -i all.ml.tre -o species_tree/coal.tre --root ZJ
+
+# bucky
+mkdir bucky
+for mbfile in ../mbsum/*.in; do
+    gene=$(basename "$mbfile" .in)
+    if ! grep -Fqx -- "$gene" ../treeshrink_filtered.list fs_exclude.list; then
+        cp "$mbfile" bucky/"$gene".in
+    fi
+done
+cd bucky
+tar -cf input.mb.tar -- ./*.in
+rm -- ./*.in
+pueue add -- "$SCRIPTS"/utils/bucky.pl -a 1 -o bucky_1 -s input.mb.tar
+
+mkdir snaq
+cd snaq
+## run snaq.jl with interactive mode
+# julia snaq.jl
+cd ../..

99.scripts_v2/05.dsuite.sh

@@ -0,0 +1,143 @@
+#!/bin/bash
+
+set -euo pipefail
+
+ln -sf ../11.reference_v2/raw .
+cd raw
+ln -sf Hippophae_tibetana.genomic.fasta HT.fa
+ln -sf YN.genome.fa YN.fa
+seqkit grep -f ZG.target.txt ZG.genome.fa >  ZG.fa
+ln -sf ZJ.genome.fa ZJ.fa
+seqkit grep -f MG.target.txt MG.genome.fasta >  MG.fa
+ln -sf EM.genome.fasta EM.fa
+cd ..
+
+bash "$SCRIPTS"/utils/asm_qc.sh raw/ZJ.fa \
+  raw/EM.fa \
+  raw/EA.fa \
+  raw/HT.fa \
+  raw/ZY.fa \
+  raw/ZG.fa \
+  raw/MG.fa \
+  raw/YN.fa \
+  raw/HG.fa \
+  raw/HS.fa \
+  raw/HN.fa > qc.tsvWDw
+
+##
+## minimap2 比对：YN 为参考（ZG 含单倍型出局），asm5 档位（属内分歧 0.001-0.007，EA/EM 约 0.03 也覆盖）
+## ZJ 已剔除（跨科，分歧 0.19）；YN 不做自比对
+##
+REF=raw/YN.fa
+mkdir -p paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HG.fa '>' paf/HG.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HN.fa '>' paf/HN.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HS.fa '>' paf/HS.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/HT.fa '>' paf/HT.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/MG.fa '>' paf/MG.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/ZG.fa '>' paf/ZG.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/ZY.fa '>' paf/ZY.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/EA.fa '>' paf/EA.paf
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/EM.fa '>' paf/EM.paf
+
+## 覆盖速览（pilot 判档用；asm5 若 EA/EM 覆盖塌，对这两个换 asm10 重跑）
+for f in paf/*.paf; do
+    awk -v OFS='\t' '{alen += $11; qlen = $2} END {printf "%s\tqlen=%d\taligned=%d\tcov=%.2f%%\n", FILENAME, qlen, alen, alen / qlen * 100}' "$f"
+done
+
+## YN 自比对：Dsuite 需要 YN 作样本（自比对产生全纯合参考基因型 VCF）
+pueue add -- minimap2 -cx asm5 --cs=long -t 8 "$REF" raw/YN.fa '>' paf/YN.paf
+
+##
+## PAF 过滤：只留 primary 链（tp:A:P）并按参考坐标排序（paftools call 要求排序输入）
+## primary 链是一对一共线块的近似；如需严格双向最佳命中，在此之上再加过滤
+##
+mkdir -p paf_primary
+grep 'tp:A:P' paf/HG.paf | sort -k6,6 -k8,8n > paf_primary/HG.paf
+grep 'tp:A:P' paf/HN.paf | sort -k6,6 -k8,8n > paf_primary/HN.paf
+grep 'tp:A:P' paf/HS.paf | sort -k6,6 -k8,8n > paf_primary/HS.paf
+grep 'tp:A:P' paf/HT.paf | sort -k6,6 -k8,8n > paf_primary/HT.paf
+grep 'tp:A:P' paf/MG.paf | sort -k6,6 -k8,8n > paf_primary/MG.paf
+grep 'tp:A:P' paf/ZG.paf | sort -k6,6 -k8,8n > paf_primary/ZG.paf
+grep 'tp:A:P' paf/ZY.paf | sort -k6,6 -k8,8n > paf_primary/ZY.paf
+grep 'tp:A:P' paf/YN.paf | sort -k6,6 -k8,8n > paf_primary/YN.paf
+grep 'tp:A:P' paf/EA.paf | sort -k6,6 -k8,8n > paf_primary/EA.paf
+grep 'tp:A:P' paf/EM.paf | sort -k6,6 -k8,8n > paf_primary/EM.paf
+
+##
+## paftools call：每样本一个 VCF（-s 样本名；-l/-L 最小块长阈值，可调）
+##
+mkdir -p vcf
+paftools.js call -s HG -l 5000 -L 10000 paf_primary/HG.paf > vcf/HG.vcf
+paftools.js call -s HN -l 5000 -L 10000 paf_primary/HN.paf > vcf/HN.vcf
+paftools.js call -s HS -l 5000 -L 10000 paf_primary/HS.paf > vcf/HS.vcf
+paftools.js call -s HT -l 5000 -L 10000 paf_primary/HT.paf > vcf/HT.vcf
+paftools.js call -s MG -l 5000 -L 10000 paf_primary/MG.paf > vcf/MG.vcf
+paftools.js call -s ZG -l 5000 -L 10000 paf_primary/ZG.paf > vcf/ZG.vcf
+paftools.js call -s ZY -l 5000 -L 10000 paf_primary/ZY.paf > vcf/ZY.vcf
+paftools.js call -s YN -l 5000 -L 10000 paf_primary/YN.paf > vcf/YN.vcf
+paftools.js call -s EA -l 5000 -L 10000 paf_primary/EA.paf > vcf/EA.vcf
+paftools.js call -s EM -l 5000 -L 10000 paf_primary/EM.paf > vcf/EM.vcf
+
+bgzip -f vcf/HG.vcf && tabix -f -p vcf vcf/HG.vcf.gz
+bgzip -f vcf/HN.vcf && tabix -f -p vcf vcf/HN.vcf.gz
+bgzip -f vcf/HS.vcf && tabix -f -p vcf vcf/HS.vcf.gz
+bgzip -f vcf/HT.vcf && tabix -f -p vcf vcf/HT.vcf.gz
+bgzip -f vcf/MG.vcf && tabix -f -p vcf vcf/MG.vcf.gz
+bgzip -f vcf/ZG.vcf && tabix -f -p vcf vcf/ZG.vcf.gz
+bgzip -f vcf/ZY.vcf && tabix -f -p vcf vcf/ZY.vcf.gz
+bgzip -f vcf/YN.vcf && tabix -f -p vcf vcf/YN.vcf.gz
+bgzip -f vcf/EA.vcf && tabix -f -p vcf vcf/EA.vcf.gz
+bgzip -f vcf/EM.vcf && tabix -f -p vcf vcf/EM.vcf.gz
+
+##
+## 合并为 10 类群多样本 VCF
+##
+pueue add -- bcftools merge -O z -o vcf/dsuite.10taxa.vcf.gz vcf/HG.vcf.gz vcf/HN.vcf.gz vcf/HS.vcf.gz vcf/HT.vcf.gz vcf/MG.vcf.gz vcf/ZG.vcf.gz vcf/ZY.vcf.gz vcf/YN.vcf.gz vcf/EA.vcf.gz vcf/EM.vcf.gz
+tabix -f -p vcf vcf/dsuite.10taxa.vcf.gz
+
+##
+## 物种树：主数据集第二套（exclude_all_fs）ASTRAL 树剔除 ZJ
+## TODO: 剪树工具任选（nw_prune / gotree / R ape），产物放 tree/dsuite.nwk
+##
+mkdir -p tree
+# nw_prune ../14.taxa_phylogeny/exclude_all_fs/species_tree/coal.tre ZJ > tree/dsuite.nwk
+
+##
+## SETS 文件：双外类群交叉验证，EA / EM 各标一次 Outgroup
+##
+mkdir -p sets dsuite
+cat > sets/EA_out.txt <<'EOF'
+HG	HG
+HN	HN
+HS	HS
+HT	HT
+MG	MG
+ZG	ZG
+ZY	ZY
+YN	YN
+EA	Outgroup
+EM	EM
+EOF
+cat > sets/EM_out.txt <<'EOF'
+HG	HG
+HN	HN
+HS	HS
+HT	HT
+MG	MG
+ZG	ZG
+ZY	ZY
+YN	YN
+EA	EA
+EM	Outgroup
+EOF
+
+##
+## Dtrios + Fbranch（EA / EM 各跑一遍）
+##
+pueue add -- Dsuite Dtrios -t tree/dsuite.nwk -o dsuite/EA_out vcf/dsuite.10taxa.vcf.gz sets/EA_out.txt
+pueue add -- Dsuite Dtrios -t tree/dsuite.nwk -o dsuite/EM_out vcf/dsuite.10taxa.vcf.gz sets/EM_out.txt
+Dsuite Fbranch tree/dsuite.nwk dsuite/EA_out_tree.txt > dsuite/EA_out_fbranch.txt
+Dsuite Fbranch tree/dsuite.nwk dsuite/EM_out_tree.txt > dsuite/EM_out_fbranch.txt
+# dtools.py dsuite/EA_out_fbranch.txt tree/dsuite.nwk
+# dtools.py dsuite/EM_out_fbranch.txt tree/dsuite.nwk

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);

pixi.toml

@@ -1,8 +1,7 @@
 [workspace]
 authors = ["IvisTang <me@ivistang.com>"]
 channels = [
-    "https://mirrors4.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/",
-    "https://mirrors4.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/",
+    "https://prefix.dev/julia-forge",
     "conda-forge",
     "bioconda"
 ]
@@ -57,19 +56,29 @@ r-phytools = ">=2.5_2,<3"
 r-argparser = ">=0.7.2,<0.8"
 r-processx = ">=3.8.6,<4"
 
-[feature.beast.dependencies]
-beagle-lib = ">=4.0.1,<5"
-beast = ">=10.5.0,<11"
-
 [feature.genome.dependencies]
 gffread = ">=0.12.9,<0.13"
 orthofinder = ">=3.1.5,<4"
 seqkit = ">=2.13.0,<3"
 tmux = ">=3.7b_,<4"
-numpy = "<2"
+numpy = "<2.5"
+python = "3.12.*"
+biopython = ">=1.87,<2"
+macse = ">=2.7,<3"
+trimal = ">=1.5.1,<2"
+treeshrink = ">=1.4.0,<2"
+raxml-ng = ">=2.0.2,<3"
+modeltest-ng = ">=0.1.7,<0.2"
+mrbayes = ">=3.2.7,<4"
+seqmagick = ">=0.8.6,<0.9"
+aster = ">=1.25,<2"
+julia = ">=1.12.6,<1.13"
+r-ape = ">=5.8_1,<6"
+r-languageserver = ">=0.3.18,<0.4"
+conda = ">=26.7.0,<27"
+minimap2 = ">=2.31,<3"
 
 
 [environments]
-base = ["base"]
-beast = ["beast"]
+# base = ["base"]
 default = ["genome"]