Comparative Genomics & Synteny

# Prepare: one protein FASTA per species in a directory
mkdir -p orthofinder_input
cp species1_proteins.fa orthofinder_input/Species1.fa
cp species2_proteins.fa orthofinder_input/Species2.fa
cp species3_proteins.fa orthofinder_input/Species3.fa

# Run OrthoFinder
orthofinder -f orthofinder_input/ -t 16 -a 8

# Key output files:
# Orthogroups/Orthogroups.tsv     — gene families
# Orthogroups/Orthogroups.GeneCount.tsv
# Comparative_Genomics/Orthogroups_SpeciesOverlaps.tsv
# Species_Tree/SpeciesTree_rooted.txt
# Orthologues/Species1__v__Species2.tsv

# Step 1: Prepare input files
# species.gff — simplified gene positions: chr  gene  start  end
awk -F'\t' '$3=="gene" {
  match($9, /ID=([^;]+)/, id);
  print $1"\t"id[1]"\t"$4"\t"$5
}' species1.gff3 species2.gff3 > combined.gff

# Step 2: All-vs-all BLAST of protein sequences
cat species1_proteins.fa species2_proteins.fa > combined_proteins.fa
makeblastdb -in combined_proteins.fa -dbtype prot
blastp -query combined_proteins.fa -db combined_proteins.fa \
  -evalue 1e-10 -outfmt 6 -num_threads 16 \
  -max_target_seqs 5 -out combined.blast

# Step 3: Run MCScanX
MCScanX combined
# Outputs: combined.collinearity, combined.html, combined.tandem

# Step 4: Classify duplications
duplicate_gene_classifier combined
# Outputs gene pairs as: WGD/segmental, tandem, proximal, dispersed

# Install
pip install jcvi

# Step 1: Prepare BED and CDS files from GFF
python -m jcvi.formats.gff bed --type=mRNA --key=Name species1.gff3 -o species1.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name species2.gff3 -o species2.bed

# Step 2: Pairwise synteny search
python -m jcvi.compara.catalog ortholog species1 species2 --no_strip_names

# Step 3: Dotplot
python -m jcvi.graphics.dotplot species1.species2.anchors

# Step 4: Macro-synteny (karyotype view)
python -m jcvi.compara.synteny screen --minspan=30 \
  species1.species2.anchors species1.species2.anchors.new

# Create layout file for macro-synteny
echo "# y, xstart, xend, rotation, color, label, va, bed
.6, .1, .8, 0, , Species 1, top, species1.bed
.4, .1, .8, 0, , Species 2, bottom, species2.bed" > layout.csv

echo "# anchor file
species1.species2.anchors.simple" > seqids

python -m jcvi.graphics.karyotype seqids layout.csv

# Step 5: Micro-synteny (zoom into a region)
python -m jcvi.graphics.synteny species1.species2.anchors.simple \
  species1.bed species2.bed --genenames

# Using circlize in R (easier than Perl Circos)
library(circlize)

# Read synteny data (from MCScanX collinearity)
synteny <- read.delim("synteny_links.tsv", header=TRUE)
# Columns: chr1, start1, end1, chr2, start2, end2

# Chromosome lengths
karyotype <- read.delim("karyotype.tsv", header=TRUE)
# Columns: chr, start, end

# Initialize circular plot
circos.clear()
circos.par(gap.degree = 2, cell.padding = c(0, 0, 0, 0))
circos.genomicInitialize(karyotype)

# Add chromosome ideogram
circos.track(ylim = c(0, 1), bg.col = rainbow(nrow(karyotype), alpha=0.3),
             panel.fun = function(x, y) {
               circos.text(CELL_META$xcenter, 0.5, CELL_META$sector.index,
                           cex = 0.6, facing = "clockwise", niceFacing = TRUE)
             })

# Add synteny links
for(i in 1:nrow(synteny)) {
  circos.link(synteny$chr1[i], c(synteny$start1[i], synteny$end1[i]),
              synteny$chr2[i], c(synteny$start2[i], synteny$end2[i]),
              col = adjustcolor("dodgerblue", alpha=0.3), border=NA)
}

# Using wgd tool
pip install wgd

# Calculate Ks for paralog pairs
wgd dmd proteins.fasta
wgd ksd cds.fasta proteins.fasta.mcl

# Plot Ks distribution
wgd viz -o ks_plot.pdf

# Using ParaAT + KaKs_Calculator (alternative)
ParaAT.pl -h homologs.txt -n cds.fa -a protein.fa -p proc_num -o paraAT_output
KaKs_Calculator -i paraAT_output/aligned.axt -o kaks_results.txt -m YN