GWAS-eQTL的共定位分析，听起来是不是很高大上。当时自己也是不太清楚，还到处找教程，脚本什么的。但是都无果，好不容易找到一个，折腾了半天安装了R包，只能做人类的！！！当场炸裂。但是其实了解了之后也就那么回事~

那就自己写一个脚本就行~

GWAS和eQTL都是对SNP进行分析，一个与表型进行关联分析，一个与基因的表达量进行关联分析；说到底都有一个共同的SNP信息。现在共定位分析只不过是将eQTL和GWAS的SNP计算的P-value值进行整合，换句话说就是在两个分析过程中有相同的SNP在两部分分析中都显示是超过阈值线的P值，那这个SNP是不是就是很重要的信息。

1.eQTL分析

可以看上一篇文章：https://www.jianshu.com/p/fba46efad6d8

2.GWAS分析

现在网上满天都是，稍微去查一下，找到适合自己的就可以做，在这里我推荐GAPIT,这个方便快捷。

3.共定位分析

文件准备

#整理eQTL结果为Chr3.eqtl.txt样式
awk 'NR > 1 {print $0}' modelANOVA.eqtl.txt | grep "Chr3" | awk '{print $1"\t"$4}' | awk -F "\t" 'BEGIN{OFS=FS}{$2=-log($2)/log(10)}1' > Chr3.eqtl.txt

$ head Chr3.eqtl.txt
rsid    pval
Chr3-35149904   42.5056
Chr3-8369353    41.4392
Chr3-35146484   39.5319
Chr3-5143404    32.4618
Chr3-35221782   31.0892
Chr3-3862546    28.8215
Chr3-4703118    27.3703
Chr3-4657827    27.3212
Chr3-4720052    27.2563
#整理GWAS结果为Chr3.eqtl.txt样式
awk '{print $1"\t"$4}' GWAS.csv | grep "Chr3" | awk -F "\t" 'BEGIN{OFS=FS}{$2=-log($2)/log(10)}1' > Chr3.gwas.txt

$ head Chr3.gwas.txt
rsid    pval
Chr3-6970       0.639774
Chr3-9305       1.15444
Chr3-19824      1.14868
Chr3-21749      0.357413
Chr3-22444      0.355889
Chr3-26997      1.0888
Chr3-29995      0.371762
Chr3-48017      0.665807
Chr3-48111      0.326685

awk -F "\t" 'BEGIN{OFS=FS}{$2=-log($2)/log(10)}1是将结果的pvalue值转化为-log10()的对数。注意文件的列名比较重要，也就是这个rsid,后面的文件要通过这个进行匹配。否则会发生错误。

共定位脚本

# 加载ggplot2包
library(ggplot2)

# 读取eQTL数据，确保文件路径正确
eqtl_data <- read.table(file="/you_path/Chr3.eqtl.txt", header=TRUE, stringsAsFactors = FALSE)
# 读取GWAS数据，确保文件路径正确
gwas_data <- read.table(file="/you_path/Chr3.gwas.txt", header=TRUE, stringsAsFactors = FALSE)

# 确保两个数据集中SNP标识符的格式相同
eqtl_data$rsid <- as.character(eqtl_data$rsid)
gwas_data$rsid <- as.character(gwas_data$rsid)

# 合并两个数据集，仅保留两个数据集中都有的SNP
merged_data <- merge(eqtl_data, gwas_data, by = "rsid")
names(merged_data) <- c("SNP", "eqtl_pval", "gwas_pval")

# 检查合并后的数据前几行
head(merged_data)

# 定义显著性阈值
gw_threshold <- 6   # GWAS的阈值
eq_threshold <- 15  # eQTL的阈值

# 创建一个新列来表示显著性状态
merged_data$significance <- with(merged_data, ifelse(eqtl_pval > eq_threshold & gwas_pval > gw_threshold, 'Red',
                                                     ifelse(gwas_pval > gw_threshold, 'Purple',
                                                     ifelse(eqtl_pval > eq_threshold, 'Green', 'Navyblue'))))

# 绘制散点图并增加显著性标记，添加居中标题和坐标轴标题
p <- ggplot(merged_data, aes(x = eqtl_pval, y = gwas_pval)) +
  geom_point(aes(color = significance), alpha = 0.5) +
  scale_color_manual(values = c("Navyblue" = "#000074", "Purple" = "#692A7D", "Green" = "#66BD99", "Red" = "#8A1C3D")) +
  geom_text(data = subset(merged_data, significance == 'Red'), aes(label = SNP), 
            hjust = 0, vjust = 0, check_overlap = TRUE, size = 3, color = "#8A1C3D") +
  geom_vline(xintercept = eq_threshold, color = "#53832F", linetype = "dashed") +
  geom_hline(yintercept = gw_threshold, color = "red", linetype = "dashed") +
  ggtitle("GWAS-eQTL co-localization analysis") +
  xlab(expression(paste("eQTL  -log"[10], "(p-value)"))) + # 设置x轴标题
  ylab(expression(paste("GWAS  -log"[10], "(p-value)"))) + # 设置y轴标题
  theme_minimal() +
  theme(legend.position = "bottom", 
        axis.text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
        axis.title.x = element_text(size = 16, face = "bold"), # 设置x轴标题样式
        axis.title.y = element_text(size = 16, face = "bold")) # 设置y轴标题样式

pdf("Chr3.pdf", width = 20, height = 12)
print(p)
dev.off()

4.结果展示

co-localization.png