一、plink软件分析LDblock

plink --bfile 142samples_basefilter_only_snp_qc \
      --chr 5 \
      --from-bp 1477439038 \
      --to-bp 1477440735 \
      --r2 square \
      --out /home/184_Tg_RNA/142samples_Tg_rna/142samples_rna_vcf_result/LDblock_result_newstyle/cyp71_ld

二、利用R来进行组合图作图

# ==========================================================================
# Figure 3E: Proportional Locus Plot (GWAS + Gene Model + Uncompressed LD)
# Core Target: evm.TU.PTG002667L.28 (Cytochrome P450, CYP71)
# ==========================================================================

library(data.table)
library(dplyr)
library(ggplot2)
library(patchwork)
library(scales)

# --------------------------------------------------------------------------
# 0. 环境与目录设置
# --------------------------------------------------------------------------
setwd("E:/课题-根腐病GWAS/基因型数据分析/转录组关联分析结果/snp_gwas_gemma_cov/")

# --------------------------------------------------------------------------
# 1. 基础输入文件与核心配置
# --------------------------------------------------------------------------
gwas_file   <- "gwas_gemma_cov_trait4_assoc.assoc.txt"
ld_file     <- "NPC_ld.ld"

target_gene <- "evm.TU.PTG004844L.25"
target_chr  <- 1
lead_pos    <- 1317891688

gene_start  <- 1317891424
gene_end    <- 1317893405

region_start <- 1317891424
region_end   <- 1317893405

# 负链基因，按基因组坐标从小到大排序绘图
exon_df <- data.frame(
  exon  = c("Exon 3", "Exon 2", "Exon 1"),
  start = c(1317891424, 1317891833, 1317892301),
  end   = c(1317891737, 1317892009, 1317893405)
)



# --------------------------------------------------------------------------
# 2. 期刊指定专属配色系统 (全图色彩完美适配曼哈顿图)
# --------------------------------------------------------------------------
col_point <- "#5679A6"   # GWAS中性变异点 (莫兰迪蓝灰)
col_blue  <- "#2F5AA8"   # 暗示阈值线 (普鲁士蓝过渡色)
col_red   <- "#D94B3D"   # Lead 变异与 LD 强连锁核心高亮 (统一学术珊瑚红)
col_gene  <- "#1B365D"   # 候选基因骨架与外显子 (深邃普鲁士蓝)
col_text  <- "#222222"   
col_band  <- "grey75"    # 垂直微观定位共享对齐带

# --------------------------------------------------------------------------
# 3. 高保真读取并清洗 GWAS 数据
# --------------------------------------------------------------------------
gwas_all <- fread(gwas_file)

if ("p_wald" %in% colnames(gwas_all)) {
  p_col <- "p_wald"
} else if ("p_lrt" %in% colnames(gwas_all)) {
  p_col <- "p_lrt"
} else if ("p_score" %in% colnames(gwas_all)) {
  p_col <- "p_score"
} else if ("p" %in% colnames(gwas_all)) {
  p_col <- "p"
} else {
  stop("No p-value column found.")
}

gwas_all <- gwas_all %>%
  rename(Chr = chr, Pos = ps) %>%
  mutate(
    Chr  = as.numeric(gsub("Chr", "", Chr)),
    Pos  = as.numeric(Pos),
    P    = as.numeric(.data[[p_col]]),
    logP = -log10(P)
  ) %>%
  filter(!is.na(Chr), !is.na(Pos), !is.na(P), P > 0)

gwas_local <- gwas_all %>%
  filter(Chr == target_chr, Pos >= region_start, Pos <= region_end) %>%
  arrange(Pos)

if (nrow(gwas_local) == 0) stop("No GWAS markers found in the selected region.")

n_marker <- nrow(gwas_all)
genome_threshold     <- -log10(0.05 / n_marker)
suggestive_threshold <- 5

lead_df <- gwas_local %>%
  mutate(dist = abs(Pos - lead_pos)) %>%
  arrange(dist) %>%
  slice(1)

ymax_assoc <- max(max(gwas_local$logP, na.rm = TRUE) + 0.8, genome_threshold + 0.3, 6.5)

# --------------------------------------------------------------------------
# 4. 高级几何重组：构建标准不压缩 45 度倒三角连锁矩阵
# --------------------------------------------------------------------------
ld_mat <- as.matrix(fread(ld_file, header = FALSE))
if (nrow(ld_mat) != ncol(ld_mat)) stop("LD matrix must be square.")

n_ld <- nrow(ld_mat)
ld_pos_scaled <- seq(region_start, region_end, length.out = n_ld)

ld_index <- which(upper.tri(ld_mat, diag = TRUE), arr.ind = TRUE) %>% as.data.frame()
colnames(ld_index) <- c("i", "j")
ld_index$R2 <- ld_mat[cbind(ld_index$i, ld_index$j)]

# 核心修正：利用实际物理区间跨度精确计算单元格的几何缩放因子，避免高度变扁
total_span <- region_end - region_start
step_bp    <- total_span / (n_ld - 1)

ld_cells <- ld_index %>%
  mutate(
    pos_i = ld_pos_scaled[i],
    pos_j = ld_pos_scaled[j],
    # 计算标准旋转 45 度后的质心 X 坐标
    x_center = (pos_i + pos_j) / 2,
    # 🔥【高度不压缩关键】：Y 轴高度与其在矩阵中的步长完全等比映射，确保完美的45度斜角形态
    y_center = -(pos_j - pos_i) / 2
  )

# 多边形闭合位移量
half_dx <- step_bp / 2
half_dy <- step_bp / 2

ld_poly_list <- lapply(seq_len(nrow(ld_cells)), function(k) {
  data.frame(
    id = k, R2 = ld_cells$R2[k],
    x = c(ld_cells$x_center[k] - half_dx, ld_cells$x_center[k], ld_cells$x_center[k] + half_dx, ld_cells$x_center[k]),
    y = c(ld_cells$y_center[k], ld_cells$y_center[k] + half_dy, ld_cells$y_center[k], ld_cells$y_center[k] - half_dy)
  )
})
ld_poly <- do.call(rbind, ld_poly_list)

# 动态提取 LD 三角形物理高度的最大边界，用于画布高精度的比例锁定
max_ld_depth <- min(ld_cells$y_center)

# ==========================================================================
# 5. Panel A: 顶部 —— 局部关联曼哈顿图 (普鲁士蓝变异信号)
# ==========================================================================
p_assoc <- ggplot(gwas_local, aes(x = Pos, y = logP)) +
  annotate("rect", xmin = gene_start, xmax = gene_end, ymin = -Inf, ymax = Inf, fill = col_band, alpha = 0.12) +
  geom_point(color = col_point, size = 2.0, alpha = 0.9) +
  geom_hline(yintercept = suggestive_threshold, color = col_blue, linewidth = 0.45, linetype = "dashed") +
  geom_hline(yintercept = genome_threshold, color = col_red, linewidth = 0.45, linetype = "dashed") +
  geom_point(data = lead_df, aes(x = Pos, y = logP), inherit.aes = FALSE, shape = 23, size = 3.5, fill = col_red, color = "black", stroke = 0.4) +
  annotate("text", x = lead_df$Pos, y = lead_df$logP + 0.35, label = target_gene, fontface = "italic", size = 3.8, color = col_text, vjust = 0) +
  scale_x_continuous(limits = c(region_start, region_end), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, ymax_assoc), expand = expansion(mult = c(0, 0.03))) +
  labs(title = "Local association signal", x = NULL, y = expression(-log[10](italic(P)))) +
  theme_classic(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 11, hjust = 0),
    axis.text.x = element_blank(), axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold", size = 11),
    axis.text.y = element_text(color = "black", size = 9),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.7),
    plot.margin = margin(6, 8, 2, 8)
  )

# ==========================================================================
# 6. Panel B: 正中间 —— 候选基因外显子结构
# ==========================================================================
p_gene <- ggplot() +
  annotate("rect", xmin = gene_start, xmax = gene_end, ymin = -Inf, ymax = Inf, fill = col_band, alpha = 0.12) +
  geom_segment(aes(x = gene_start, xend = gene_end, y = 1, yend = 1), color = col_gene, linewidth = 1.0) +
  geom_rect(data = exon_df, aes(xmin = start, xmax = end, ymin = 0.78, ymax = 1.22), fill = col_gene, color = col_gene) +
  geom_text(data = exon_df, aes(x = (start + end) / 2, y = 1.35, label = exon), size = 3.0, color = col_gene, fontface = "bold") +
  annotate("text", x = gene_end + (region_end - region_start) * 0.02, y = 1.0, 
           label = paste0(target_gene, "\ncytochrome P450, family 71"), 
           hjust = 0, vjust = 0.5, size = 3.8, color = col_text, fontface = "italic") +
  scale_x_continuous(limits = c(region_start, region_end), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0.4, 1.6)) +
  labs(title = "Gene structure", x = NULL, y = NULL) +
  theme_classic(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 11, hjust = 0),
    axis.text.x = element_blank(), axis.ticks.x = element_blank(),
    axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.7),
    plot.margin = margin(2, 8, 2, 8)
  )

# ==========================================================================
# 7. Panel C: 最底部 —— 饱满不压缩的珊瑚红标准 LD 三角热图
# ==========================================================================
p_ld <- ggplot(ld_poly, aes(x = x, y = y, group = id, fill = R2)) +
  annotate("rect", xmin = gene_start, xmax = gene_end, ymin = -Inf, ymax = Inf, fill = col_band, alpha = 0.08) +
  geom_polygon(color = "grey93", linewidth = 0.15) +
  # 🔥【色彩搭配微调】：采用和曼哈顿图完全一体化的珊瑚红渐变，摒弃刺眼高饱和纯红
  scale_fill_gradient(low = "#FFFFFF", high = col_red, limits = c(0, 1), name = expression(italic(r)^2)) +
  scale_x_continuous(limits = c(region_start, region_end), labels = comma, expand = c(0, 0)) +
  # 🔥【精准控高机制】：动态绑定 Y 轴的真实映射深度，防止网格因拼图产生扁平化形变
  scale_y_continuous(limits = c(max_ld_depth - half_dy, 0 + half_dy), expand = c(0, 0)) +
  labs(title = "LD pattern", x = paste0("Genomic Position on Chromosome ", target_chr, " (bp)"), y = NULL) +
  theme_void(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 11, hjust = 0, color = "black"),
    axis.title.x = element_text(face = "bold", size = 11, margin = margin(t = 6)),
    axis.text.x = element_text(color = "black", size = 9, face = "bold"),
    legend.position = "right",
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 8),
    plot.margin = margin(2, 8, 6, 8)
  )

# --------------------------------------------------------------------------
# 8. 顶级多维组图组装控制中心 (Patchwork 协同比例微调)
# --------------------------------------------------------------------------
# 🔥 通过调整 heights 里的第三项参数（从 2.2 增加到 3.5），
# 为底部的 LD 热图提供更充裕的垂直像素空间，配合等比纵向质心映射，让倒三角挺拔饱满！
p_final <- p_assoc / p_gene / p_ld +
  plot_layout(heights = c(2.5, 0.9, 3.5)) + 
  plot_annotation(
    title = paste0("Representative candidate locus: Chr", target_chr, " - ", target_gene),
    theme = theme(
      plot.title = element_text(face = "bold", size = 14, hjust = 0.5, color = "black")
    )
  )

# --------------------------------------------------------------------------
# 9. 导出高DPI无损定稿图
# --------------------------------------------------------------------------
ggsave("Fig3E_PTG004844L25_uncompressed_perfect.pdf", p_final, width = 8.5, height = 7.8)
ggsave("Fig3E_PTG004844L25_uncompressed_perfect.png", p_final, width = 8.5, height = 7.8, dpi = 600)

message("🎉 调整完毕！饱满不压缩且色彩深度对齐的高级关联扫描图已安全导出。")