PHYLOGENOMICS OF Brachyotum

From assembly of raw data to phylogenetic analyses
For nuclear (target capture) and plastid (genome skimming) data sets

Pipeline for phylogenomic analyses for the paper Paredes-Burneo et al. (2026). Evidence of introgression amid phylogenetic conflict in Brachyotum, a plant radiation from the Tropical Andes. American Journal of Botany. From data assembly to multiple phylogenetic analyses.

Raw data and their resulting processed data sets stored on Zenodo. You may use code in this GitHub repository to replicate the analyses.

Software dependencies:

ASTRAL (ASTER: ASTRAL IV and weighted ASTRAL)
BEAST 2.7.6 (and the accompanying programs, such as BEAUti)
CAPTUS 1.0.1
Geneious 2024.0
HybPhaser 2.0
IQ-TREE 2.2.0.3
julia
phyx
python
R
SNaQ (PhyloNetworks)
Trimmomatic 0.39

ASSEMBLY (AND PREVIOUS CLEANING)

Assembly of nuclear loci (sequenced with target capture)

Previous cleaning of barcodes: trimming_Brach.sh
If there were unpaired files after running trimming_Brach.sh, then proceed with:

for infile in ./*R1un.trim.fq
do
base=$(basename ${infile} R1un.trim.fq)
cat ${infile} ${base}R2un.trim.fq > ${base}UN.fq
done

Performing assembly

while read name;
do hybpiper assemble -t_dna seqMELAS.fasta -r $name*.trim.fq --prefix $name --bwa --unpaired $name*UN.fq --cpu 10 --run_intronerate;
done < namelist.txt

Retrieving sequences

hybpiper retrieve_sequences -t_dna seqMELAS.fasta dna --sample_names namelist.txt

Retrieving long putative paralogs (including non-chimeric)

hybpiper paralog_retriever namelist.txt -t_dna seqMELAS.fasta

Assembly stats

hybpiper stats -t_dna seqMELAS.fasta gene namelist.txt

Recovery visualization heatmap_paralogs.R

Assembly of plastid loci (sequenced with genome skimming)

ALIGNMENT, FILTERING AND HETEROZYGOSITY ASSESSMENT

Alignment of nuclear loci

As a parallel task:

for filename in $(ls *.fasta)
do
echo macse -prog alignSequences -seq $filename -out_NT $(cut -d'.' -f1 <<<"$filename").NT.aln -out_AA $(cut -d'.' -f1 <<<"$filename").AA.aln > $(cut -d'.' -f1 <<<"$filename").sh
done


parallel -j 10 bash ::: *.sh

In case one locus does not align:

mafft --auto NM_179236.3.FNA > NM_179236.3.FNA.aln

Data filtering of the nuclear data set

For inference of orthologous sequences: custom script based on Morales-Briones' script

To rename putative copies within alignment files for the outgroup:

for filename in ${filename#/}
do
awk '/^>/ {printf("\n%s@",$0);next; } { printf("%s",$0);} END {printf("\n");}' $filename |
grep -v "Tlepidota.[0-9]\+" |
sed 's/@/\n/' > @$filename;
done

To remove sequences for one contaminated sample across loci

for file in *.fna;do
pxrms -s "$file" -r *Brachyotum_russatum* -o "$file".pxrms;
done

To remove missing data per site from alignments

pxclsq -s *pxrms -o *pxclsq -p 0.6

Heterozygosity (and allele divergence) assessment of the nuclear data set

Following the HybPhaser pipeline.

Alignment and data filtering of plastid loci

Following the CAPTUS pipeline.

PHYLOGENETIC INFERENCE

Under a concatenation approach (nuclear)

Alignments concatenated with Geneious.

iqtree -s nuclear_concat_brachyotum.fa -m GTR+G -bb 1000 -alrt 1000 -bnni -nt AUTO

Under a summary-method approach (nuclear)

First, Gene trees with IQ-TREE

for filename in $(ls *.fa-cln);
do echo iqtree -s $filename -m GTR+G -bb 1000 -alrt 1000 -bnni -nt AUTO > $filename.sh;
done
parallel -j 15 bash ::: *.sh

Then, gene tree reconciliation with ASTRAL

~/ASTER-Linux/bin/wastral -t 15 -R -i INPUT -o OUTPUT 2> LOG
~/ASTER-Linux/bin/astral4 -u 2 -R -i INPUT -o OUTPUT 2> LOG

Under a co-estimation approach (nuclear)

To start an independent run

./beast/bin/beast brachyotum.xml -working -seed 123456789 -threads 20

To resume analysis (mainly on a cluster)

./beast/bin/beast -seed 123456789 -threads 20 -resume brachyotum.xml

To combine log files from independent runs

./beast/bin/logcombiner -log ../run1/run1.log -log ../run2/run2.log -b 0 -o combined_run.log

Phylogenetic network inference

using PhyloNetworks
using PhyloPlots
using Distributed
addprocs(10)
@everywhere using PhyloNetworks
genetrees = readTrees2CF("starbeast3.trees.nw", writeTab = false, writeSummary = false)
msc = readTopology("starbeast3.tre.nw")
net0 = snaq!(msc, genetrees, hmax = 0, filename = "net0", seed = 1234)
net1 = snaq!(net0, genetrees, hmax = 1, filename = "net1", seed = 1234)
net2 = snaq!(net1, genetrees, hmax = 2, filename = "net2", seed = 1234)
net3 = snaq!(net2, genetrees, hmax = 3, filename = "net3", seed = 1234)
net4 = snaq!(net3, genetrees, hmax = 4, filename = "net4", seed = 1234)

Phylogenetic inference for the plastid data set

iqtree2 -p /ptd --prefix ptd -m MFP+MERGE -alrt 1000 -bb 1000 -bnni -nt AUTO

ASSESSMENT OF GENE TREE DISCORDANCE

Quartet sampling scores

~/quartetsampling-master/pysrc/quartet_sampling.py --tree concatenation.tre.nw --align nuclear_concat_brachyotum.phy --reps 10000 --threads 5 --lnlike 2 --engine iqtree --data-type nuc --engine-model GTR+G

Robinson-Foulds distances

Custom script to calculate and plot RF distances RFdistances.R

PLOTTING

To plot multiple values per branch on topologies (concatenation and plastid trees) tree_multiple_node.labels.R
Tanglegram: tanglegram_msc.R
Pseudolikelihood values from Heuristic search for Phylogenetic network: H_snaq.R