Plant, Polyploid & Non-Model Genomics

# 1. Count canonical 21-mers using KMC3
mkdir -p kmc_tmp
kmc -k21 -t16 -m32 -ci1 -cs10000 @hifi_reads.fofn kmer_counts kmc_tmp/
kmc_tools transform kmer_counts histogram kmer_k21.hist -cx10000

# 2. Run GenomeScope 2.0 (set -p 2 for diploid, -p 4 for tetraploid, -p 6 for hexaploid)
genomescope2 -i kmer_k21.hist -o genomescope_out -k 21 -p 2

# 3. Extract heterozygous k-mer pairs for Smudgeplot ploidy inference
L=$(smudgeplot.py cutoff kmer_k21.hist L)
U=$(smudgeplot.py cutoff kmer_k21.hist U)
kmc_tools transform kmer_counts -ci"$L" -cx"$U" dump -s kmer_L_U.dump
smudgeplot.py hetkmers -o kmer_pairs < kmer_L_U.dump
smudgeplot.py plot kmer_pairs_coverages.tsv -o plant_ploidy

# Haplotype-resolved assembly using PacBio HiFi + paired-end Hi-C reads
hifiasm -o cultivar_asm -t 48 \
  --h1 hic_R1.fastq.gz \
  --h2 hic_R2.fastq.gz \
  pacbio_hifi.fastq.gz

# Convert GFA assembly graphs to FASTA for Primary, Hap1, and Hap2
for hap in p_ctg hap1.p_ctg hap2.p_ctg; do
  awk '/^S/{print ">"$2"\n"$3}' cultivar_asm.hic.${hap}.gfa > cultivar_${hap}.fasta
done

# 1. Align Hi-C reads with BWA-MEM2 (-5SP flags preserve 5' ligation junctions)
bwa-mem2 index cultivar_p_ctg.fasta
bwa-mem2 mem -5SP -t 32 cultivar_p_ctg.fasta hic_R1.fastq.gz hic_R2.fastq.gz | \
  samblaster | samtools view -Sb -h -F 2316 - > hic_aligned.bam

# 2. Sort by read name and scaffold with YaHS
samtools sort -n -@ 16 -o hic_namesorted.bam hic_aligned.bam
samtools faidx cultivar_p_ctg.fasta
yahs cultivar_p_ctg.fasta hic_namesorted.bam -o yahs_out

# 1. BUSCO Completeness against Embryophyta (Land Plants) OrthoDB v10
busco -i yahs_out_scaffolds_final.fa -l embryophyta_odb10 -o busco_plant -m genome -c 16

# 2. Merqury Reference-Free Consensus Quality (QV) & K-mer Completeness
meryl k=21 count output read_db.meryl pacbio_hifi.fastq.gz
merqury.sh read_db.meryl yahs_out_scaffolds_final.fa merqury_qc

# 3. Calculate LTR Assembly Index (LAI) via LTR_retriever
LTR_retriever -genome yahs_out_scaffolds_final.fa -inharvest raw_ltr.scn -threads 16

# Run Extensive de-novo TE Annotator (EDTA) and soft-mask the genome
EDTA.pl \
  --genome yahs_out_scaffolds_final.fa \
  --species others \
  --step all \
  --sensitive 1 \
  --anno 1 \
  --threads 32

# Generate soft-masked genome for BRAKER3 (--xsmall converts repeats to lowercase a/t/g/c)
RepeatMasker -pa 32 -s -xsmall \
  -lib yahs_out_scaffolds_final.fa.mod.EDTA.TElib.fa \
  yahs_out_scaffolds_final.fa

# Run BRAKER3 via Singularity/Apptainer container with RNA-Seq BAMs and Viridiplantae proteins
singularity exec -B $PWD:$PWD docker://teambraker/braker3:latest braker.pl \
  --genome=yahs_out_scaffolds_final.fa.masked \
  --bam=leaf_rnaseq.bam,root_rnaseq.bam,fruit_rnaseq.bam \
  --prot_seq=Viridiplantae_OrthoDB11.fa \
  --threads=32 \
  --gff3 \
  --species=solanum_cultivar_v1 \
  --workingdir=braker3_out

# 1. Fast orthology-based functional annotation using EggNOG-mapper (DIAMOND mode)
emapper.py -i braker3_out/braker.aa \
  --itype proteins \
  -m diamond \
  --tax_scope 33090 \
  --cpu 24 \
  --output plant_proteome_func

# 2. Domain & family architecture scan via InterProScan
interproscan.sh -i braker3_out/braker.aa -f tsv,gff3 -dp -goterms -pa -cpu 24

# 1. All-vs-all DIAMOND BLASTP for collinearity detection
diamond makedb --in combined_proteomes.fa -d combined_db
diamond blastp -d combined_db -q combined_proteomes.fa -o all_vs_all.blast -e 1e-10 --max-target-seqs 5 -p 16

# 2. Identify syntenic collinear blocks with MCScanX
MCScanX species_comparison

# 3. Plot macrosynteny karyotype ribbons using JCVI
python -m jcvi.compara.catalog ortholog tomato_v1 potato_v1 --no_strip_names
python -m jcvi.graphics.karyotype seqids layout

# 1. Build an all-to-all chromosome pan-genome graph across 8 accessions using PGGB
pggb -i chr01_8_accessions.fa.gz \
  -p 95 \
  -s 10000 \
  -n 8 \
  -t 32 \
  -V 'ref_cultivar:#' \
  -o chr01_pangenome_graph

# 2. Call chromosomal inversions, translocations & PAVs between two assemblies with SyRI
minimap2 -ax asm5 --eqx -t 16 ref_genome.fa query_genome.fa | samtools sort -o asm_align.bam
syri -c asm_align.bam -r ref_genome.fa -q query_genome.fa -F B --prefix cultivar_vs_wild_
plotsr --sr cultivar_vs_wild_syri.out --genomes genomes.txt -o structural_variants_chr.pdf

# 1. Index an all-sites VCF (containing both variant and invariant/monomorphic sites)
tabix -p vcf crop_panel_allsites.vcf.gz

# 2. Compute unbiased sliding-window Pi, Fst, and dXY between Wild vs. Cultivated groups
pixy --stats pi fst dxy \
  --vcf crop_panel_allsites.vcf.gz \
  --populations pop_assignments.txt \
  --window_size 50000 \
  --n_cores 16 \
  --output_folder pixy_domestication_scan

# R Production Script: Multi-Locus GWAS using GAPIT 3 (BLINK & FarmCPU)
library(GAPIT)

myY <- read.table("phenotypes_yield_drought.txt", head = TRUE)
myG <- read.delim("hapmap_filtered_maf005.hmp.txt", head = FALSE)

myGAPIT <- GAPIT(
  Y = myY[, c(1, 2)],          # Taxa ID + Trait column
  G = myG,                     # HapMap genotype matrix
  PCA.total = 3,               # Top 3 PCs for population structure (Q)
  model = c("MLM", "FarmCPU", "BLINK"),
  Multiple_analysis = TRUE
)

# R Production Script: Sliding-Window Δ(SNP-index) and G' Analysis with QTLseqr
library(QTLseqr)

df <- importFromGATK(
  file = "bsa_bulks_gatk.table",
  highBulk = "Resistant_Bulk",
  lowBulk = "Susceptible_Bulk"
)

df_filt <- filterSNPs(
  SNPset = df,
  refAlleleFreq = 0.20,
  minTotalDepth = 40,
  maxTotalDepth = 250,
  minSampleDepth = 20
)

df_qtl <- runQTLseqAnalysis(df_filt, windowSize = 1e6, popStruc = "F2", bulkSize = c(40, 40), replications = 10000)
df_qtl <- runGprimeAnalysis(df_qtl, windowSize = 1e6, outlierFilter = "deltaSNP")
plotQTLStats(SNPset = df_qtl, var = "deltaSNP", plotIntervals = TRUE)

# Assemble circular Chloroplast genome directly from WGS Illumina reads
get_organelle_from_reads.py \
  -1 wgs_R1.fq.gz -2 wgs_R2.fq.gz \
  -o plastome_out \
  -R 15 -k 21,45,65,85,105 \
  -F embplant_pt \
  -t 16