167_植物基因组分析¶

# === 安装EDTA ===
conda install -c bioconda edta

# === 运行EDTA ===
EDTA.pl \
  --genome genome.fasta \
  --species others \
  --step all \
  --sensitive 1 \
  --anno 1 \
  --threads 32

# 参数:
# --species: Rice/Maize/others
# --sensitive 1: 高灵敏度模式
# --anno 1: 输出注释(mask文件等)
# --cds: 提供CDS文件排除基因(推荐)

# === 带CDS过滤(推荐) ===
EDTA.pl \
  --genome genome.fasta \
  --cds known_cds.fasta \
  --species others \
  --step all \
  --sensitive 1 \
  --anno 1 \
  --threads 32

# === 输出文件 ===
# genome.fasta.mod.EDTA.TElib.fa    : TE consensus库
# genome.fasta.mod.EDTA.TEanno.gff3 : TE注释GFF3
# genome.fasta.mod.EDTA.TEanno.sum  : 统计摘要
# genome.fasta.mod.MAKER.masked     : 硬mask基因组

# === 统计TE组成 ===
# EDTA输出包含分类统计:
# LTR retrotransposons: Copia, Gypsy
# DNA transposons: hAT, CACTA, Mutator, Harbinger
# MITE, Helitron, LINE, SINE

# === 安装BRAKER3 ===
# BRAKER3整合了GeneMark-ETP和Augustus
conda install -c bioconda braker3

# === 方法1: 仅用蛋白证据 ===
braker.pl \
  --genome=genome_masked.fasta \
  --prot_seq=proteins.fasta \
  --softmasking \
  --threads=32 \
  --species=my_plant

# === 方法2: RNA-seq + 蛋白证据(推荐) ===
# 先比对RNA-seq
hisat2-build genome_masked.fasta genome_idx
hisat2 -x genome_idx \
  -1 rnaseq_R1.fq.gz -2 rnaseq_R2.fq.gz \
  -p 32 --dta | samtools sort -@ 8 -o rnaseq.bam
samtools index rnaseq.bam

# 运行BRAKER3
braker.pl \
  --genome=genome_masked.fasta \
  --bam=rnaseq.bam \
  --prot_seq=OrthoDB_plants.fasta \
  --softmasking \
  --threads=32 \
  --species=my_plant \
  --gff3

# === 方法3: BRAKER结合Iso-Seq长读长 ===
braker.pl \
  --genome=genome_masked.fasta \
  --bam=rnaseq.bam \
  --prot_seq=proteins.fasta \
  --rnaseq_sets_ids=isoseq \
  --rnaseq_sets_dirs=isoseq_aligned.bam \
  --softmasking \
  --threads=32

# === 输出 ===
# braker/augustus.hints.gff3 : 基因注释
# braker/augustus.hints.aa   : 蛋白序列
# braker/augustus.hints.codingseq : CDS序列

# === 安装BUSCO ===
conda install -c bioconda busco

# === 基因组完整性评估 ===
busco -i genome.fasta \
  -l embryophyta_odb10 \
  -o busco_genome \
  -m genome \
  -c 16

# === 蛋白/转录本完整性评估 ===
busco -i predicted_proteins.fasta \
  -l embryophyta_odb10 \
  -o busco_proteins \
  -m protein \
  -c 16

# === 解读BUSCO结果 ===
# C:95.2%[S:92.1%,D:3.1%],F:2.3%,M:2.5%,n:1614
# C = Complete (完整)
# S = Single-copy (单拷贝完整)
# D = Duplicated (重复完整) - 多倍体中D高是正常的
# F = Fragmented (片段化)
# M = Missing (缺失)
# 质量标准: C>95%为优秀; C>90%为良好

# === 多倍体注意 ===
# 六倍体小麦的BUSCO Duplicated可能>50%(因为3套同源基因)
# 这不是错误，而是多倍体的特征

# === 共线性分析(MCScanX) ===
# 安装
git clone https://github.com/wyp1125/MCScanX.git
cd MCScanX && make

# 准备输入: 全蛋白BLASTp
makeblastdb -in all_proteins.fasta -dbtype prot -out blast_db
blastp -query all_proteins.fasta -db blast_db \
  -outfmt 6 -evalue 1e-10 -num_threads 16 \
  -out blast.out -max_target_seqs 5

# 准备GFF简化文件 (chr, gene, start, end)
awk '$3=="gene"' annotation.gff3 | \
  awk -F'\t' '{print $1"\t"$9"\t"$4"\t"$5}' > genes.gff

# 运行MCScanX
MCScanX genes  # 输入前缀

# 输出:
# genes.collinearity : 共线性块
# genes.tandem       : 串联重复基因
# genes.html         : 网页可视化

# === Ks分析检测WGD ===
# 使用wgd包或手动计算
library(ggplot2)

# 读取Ks值(从ParaAT+KaKs_Calculator获得)
ks_data <- read.delim("ks_values.tsv")

# Ks分布图
ggplot(ks_data, aes(x = Ks)) +
  geom_density(fill = "steelblue", alpha = 0.5) +
  xlim(0, 5) +
  geom_vline(xintercept = c(0.3, 0.8, 1.5), linetype = "dashed", color = "red") +
  annotate("text", x = c(0.3, 0.8, 1.5), y = rep(0.8, 3),
           label = c("Recent WGD", "Ancient WGD", "Older WGD")) +
  labs(x = "Ks (synonymous substitution rate)", y = "Density",
       title = "Ks Distribution - WGD Detection") +
  theme_minimal()
# 峰值对应全基因组复制事件时间点

# === 植物基因组注释完整流程 ===
# 1. 转座子注释
EDTA.pl --genome genome.fa --cds cds.fa --sensitive 1 --anno 1 --threads 32

# 2. 软mask基因组
bedtools maskfasta -fi genome.fa -bed EDTA_TE.bed -soft -fo genome_softmasked.fa

# 3. 基因预测
braker.pl --genome=genome_softmasked.fa --bam=rnaseq.bam \
  --prot_seq=proteins.fa --softmasking --threads=32 --gff3

# 4. BUSCO评估
busco -i braker/augustus.hints.aa -l embryophyta_odb10 -m protein -c 16

# 5. 功能注释
interproscan.sh -i predicted_proteins.fa -goterms -pa -f tsv