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
