613_比较基因组共线性分析¶

# === 从GitHub安装MCScanX ===
git clone https://github.com/wyp1125/MCScanX.git  # 克隆仓库
cd MCScanX  # 进入目录
make  # 编译

# 将MCScanX加入PATH
echo 'export PATH=$PATH:/path/to/MCScanX' >> ~/.bashrc  # 添加路径
source ~/.bashrc  # 刷新

# 验证安装
MCScanX -h  # 查看帮助

# JCVI是更现代的Python替代工具，比MCScanX快
pip install jcvi  # pip安装

# 依赖安装
conda install -c bioconda last  # LAST比对工具（比BLAST快）
pip install matplotlib  # 画图依赖

# 验证安装
python -m jcvi.compara.catalog --help  # 查看帮助

# === 准备GFF文件（基因位置文件）===
# MCScanX的GFF格式很简单：染色体 基因名 起始 终止
# 格式：chr\tgene_id\tstart\tend

# 从标准GFF3提取
awk '$3=="gene" {
    split($9, a, ";");     # 分割第9列
    split(a[1], b, "=");   # 获取ID
    gsub(/ID=/, "", b[2]); # 去掉前缀
    print $1"\t"b[2]"\t"$4"\t"$5  # 输出4列
}' genome_A.gff3 > species_A.gff  # 物种A的基因位置

awk '$3=="gene" {
    split($9, a, ";");
    split(a[1], b, "=");
    gsub(/ID=/, "", b[2]);
    print $1"\t"b[2]"\t"$4"\t"$5
}' genome_B.gff3 > species_B.gff  # 物种B的基因位置

# 合并两个物种的GFF
cat species_A.gff species_B.gff > input.gff  # 合并

# === 准备蛋白序列 ===
cat proteins_A.fasta proteins_B.fasta > all_proteins.fasta  # 合并蛋白序列

# === 运行BLAST ===
# 先建数据库
makeblastdb -in all_proteins.fasta \
    -dbtype prot \             # 蛋白数据库
    -out all_proteins_db       # 数据库名

# 全比全BLASTp
blastp \
    -query all_proteins.fasta \     # 查询序列
    -db all_proteins_db \           # 数据库
    -out input.blast \              # 输出文件
    -outfmt 6 \                     # 表格格式（MCScanX要求）
    -evalue 1e-10 \                 # E值阈值
    -num_threads 16 \               # 线程数
    -max_target_seqs 5 \            # 每个基因最多5个匹配
    -num_alignments 5               # 最多5个比对

# === 基本运行 ===
# 注意：input.gff和input.blast必须同名（前缀相同）
MCScanX input  # 运行共线性分析

# === 自定义参数 ===
MCScanX input \
    -a \          # 检测种间和种内共线性（默认只检测种内）
    -e 1e-10 \    # E值阈值
    -s 5 \        # 最少锚定基因对数（默认5）
    -m 25 \       # 最大间隔基因数（默认25）
    -b 2          # 匹配分值（默认50）

# 主要输出文件
# input.collinearity — 共线性块信息
# input.tandem — 串联重复基因
# input.html — HTML可视化

# 查看共线性块
head -30 input.collinearity

# 输出格式示例：
# ## Alignment 0: score=350 e_value=3.7e-20 N=7 ChrA&ChrB plus
# 0- 0: geneA_001 geneB_001   1e-50
# 0- 1: geneA_002 geneB_002   2e-45
# ...
# 含义：第0号共线性块，得分350，包含7对锚定基因，
# 位于ChrA和ChrB之间，方向相同（plus）

# === 1. 点阵图（Dot plot）===
java -cp MCScanX/downstream_analyses dot_plotter \
    -g input.gff \              # 基因位置
    -s input.collinearity \     # 共线性结果
    -c dot_plot_config.txt \    # 配置文件
    -o dot_plot.png             # 输出图片

# === 2. 圈图（Circle plot）===
java -cp MCScanX/downstream_analyses circle_plotter \
    -g input.gff \
    -s input.collinearity \
    -c circle_config.txt \
    -o circle_plot.png

# === 3. 双基因组对比图（Dual synteny plot）===
java -cp MCScanX/downstream_analyses dual_synteny_plotter \
    -g input.gff \
    -s input.collinearity \
    -c dual_config.txt \
    -o dual_synteny.png

# === 4. 柱形对比图（Bar plot）===
java -cp MCScanX/downstream_analyses bar_plotter \
    -g input.gff \
    -s input.collinearity \
    -c bar_config.txt \
    -o bar_plot.png

# === 第1步：准备BED和CDS文件 ===
# 从GFF3和基因组提取BED格式和CDS序列
python -m jcvi.formats.gff bed \
    --type=mRNA \              # 提取mRNA特征
    --key=Name \               # 用Name作为基因ID
    genome_A.gff3 \            # GFF3文件
    -o species_A.bed           # 输出BED文件

python -m jcvi.formats.gff bed \
    --type=mRNA \
    --key=Name \
    genome_B.gff3 \
    -o species_B.bed

# 提取CDS序列
python -m jcvi.formats.fasta format \
    cds_A.fasta species_A.cds  # 格式化CDS文件
python -m jcvi.formats.fasta format \
    cds_B.fasta species_B.cds

# === 第2步：运行共线性分析 ===
python -m jcvi.compara.catalog ortholog \
    species_A species_B \      # 两个物种
    --no_strip_names           # 保留完整基因名

# 输出文件：
# species_A.species_B.anchors  — 锚定基因对
# species_A.species_B.lifted.anchors — 扩展的锚定对
# species_A.species_B.last — LAST比对结果

# === 第3步：画点阵图（Dot plot）===
python -m jcvi.graphics.dotplot \
    species_A.species_B.anchors  # 画点阵图

# === 第4步：画共线性对比图 ===
# 先创建布局文件（seqids和layout）

# seqids文件：每行列出要展示的染色体
echo "Chr1,Chr2,Chr3,Chr4,Chr5" > seqids   # 物种A的染色体
echo "chr1,chr2,chr3,chr4,chr5" >> seqids   # 物种B的染色体

# layout文件：定义图形布局
cat > layout << 'EOF'
# y, xstart, xend, rotation, color, label, va, bed
.6, .1, .8, 0, , Species_A, top, species_A.bed
.4, .1, .8, 0, , Species_B, top, species_B.bed
# edges
e, 0, 1, species_A.species_B.anchors.simple
EOF

# 画宏观共线性图
python -m jcvi.graphics.karyotype \
    seqids layout  # 输出 karyotype.pdf

# === 用JCVI计算Ks值 ===
python -m jcvi.compara.synteny mcscan \
    species_A.bed species_A.species_B.anchors \
    --iter=1  # 迭代次数

# === 画Ks分布图 ===
python -m jcvi.compara.synteny ksdistrib \
    species_A.species_B.anchors.ks  # Ks值文件

# === R语言画Ks密度图 ===
library(ggplot2)  # 加载ggplot2

# 读取Ks数据
ks_data <- read.table("ks_values.txt", header = TRUE)  # 读取Ks值

# 画Ks密度分布图
ggplot(ks_data, aes(x = Ks)) +
    geom_density(fill = "steelblue", alpha = 0.5) +  # 密度曲线
    xlim(0, 3) +  # x轴范围
    geom_vline(xintercept = c(0.5, 1.0, 1.5),  # 标记可能的WGD峰
               linetype = "dashed", color = "red") +
    labs(
        title = "Ks Distribution",
        x = "Synonymous substitution rate (Ks)",  # 同义替换率
        y = "Density"  # 密度
    ) +
    theme_minimal()

ggsave("ks_distribution.pdf", width = 8, height = 5)  # 保存

# Ks峰值对应WGD事件
# 如果有明显的峰，说明发生过全基因组复制