单倍型、候选基因关联分析

原创已于 2025-05-12 15:40:16 修改 · 845 阅读

7 ·

CC 4.0 BY-SA版权

文章标签：

#r语言

于 2025-01-12 22:20:00 首次发布

对于现成的VCF 这里补充一下单倍型整理的方法：

在直接有vcf文件的情况下，可以直接用R先获取这样的单倍型信息：

library(geneHapR)
setwd("")

vcf <- import_vcf("v1.vcf")
hapPrefix <- "H"        # 单倍型名称的前缀

# 从VCF开始单倍型鉴定
hapResult <- vcf2hap(vcf, hapPrefix = hapPrefix,
                     hetero_remove = TRUE, # 移除包含杂合位点的样本
                     na_drop = TRUE) # 移除包含基因型缺失的样本

# 对单倍型结果进行汇总整理
hapSummary <- hap_summary(hapResult, hapPrefix = hapPrefix)
write.hap(hapSummary, file = "nrt.hapSummary")

1.把单倍型文件处理成VCF文件

这里的单倍型信息形式可能非常多，要根据具体的数据类型来进行输入信息的处理vcf hapmap等皆可主要还是要包含材料位点染色体等需要结合分析的信息，这里随便选用一个基因与他的单倍型来进行测试。

import pandas as pd
import numpy as np


def create_vcf_input(genome_df, pheno_df):
    # 转置基因型数据并处理
    geno = genome_df.T

    # 创建VCF头部信息
    vcf_header = """##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=2>
##reference=B73_RefGen_v4
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{}""".format('\t'.join(pheno_df['Sample'].tolist()))

    # 处理每个SNP位点
    vcf_lines = []
    for idx, row in geno.iterrows():
        if 'REF' in idx:  # 确保是SNP行
            # 解析位置和等位基因信息
            pos = idx.split('(')[0]
            ref_alt = idx.split('(REF:')[1].split('/ALT:')
            ref = ref_alt[0]
            alt = ref_alt[1].strip(')')

            # 创建VCF行
            genotypes = []
            for gt in row[1:]:  # 跳过第一列（Unnamed: 0）
                if gt == ref:
                    genotypes.append('0/0')
                elif gt == alt:
                    genotypes.append('1/1')
                else:
                    genotypes.append('./.')

            vcf_line = f"2\t{pos}\t.\t{ref}\t{alt}\t.\tPASS\t.\tGT\t{'\t'.join(genotypes)}"
            vcf_lines.append(vcf_line)

    # 合并所有内容
    vcf_content = vcf_header + '\n' + '\n'.join(vcf_lines)

    # 保存为文件
    with open('Zm00001d002439.vcf', 'w') as f:
        f.write(vcf_content)


# 读取文件（使用你的原始路径）
f_g = pd.read_csv(r"Zm00001d002439_DI_SX\Zm00001d002439.genome.csv")
f_pht = pd.read_csv(r"Zm00001d002439_DI_SX\Zm00001d002439.Phe.csv")

# 创建VCF文件
create_vcf_input(f_g, f_pht)

2.表型文件的准备

对应的也就是Tassel中的Phenotype文件

这里的格式要求非常诡异试了好久

<Trait>  trait1
773-2G  56.55270655
169 31.31313131
391 68.88888889
706 16.23931624
Fu746Mu 73.33333333
Ji63    23.14814815
.....

3.Tassel

接下来的操作就比较简单了：主要是gui版的Tassel

读入两个文件——>

vcf+phe——>

选中两个文件

——>直接用GLM分析就好了

——>

但是Tassel的图片太粗糙了，所以使用R代码自己绘图

在这个曼哈顿图中，90% quantile = 0.676 的含义是：
在这个曼哈顿图中，90% quantile = 0.676 的含义是：

1. 统计学解释：
- 在所有的 -log10(p) 值中，90%的点都小于0.676
- 换句话说，只有10%的SNP位点的 -log10(p) 值高于0.676

2. 生物学意义：
- 这个阈值可以帮助我们识别相对更显著的关联信号
- 高于这个阈值的位点（约10%的SNP）可能代表与性状有更强关联的候选区域
- 在您的数据中，有几个位点的显著性明显高于这个阈值，特别是接近 -log10(p) = 0.9 的那些位点

3. 实际应用：
- 这不是严格的显著性阈值（像GWAS中常用的 p < 0.05 或更严格的阈值）
- 而是一个相对的参考线，帮助识别在这个候选区域中信号相对较强的位点
- 对于候选基因分析，这个阈值可以帮助我们关注信号较强的区域进行后续功能验证

需要注意的是，因为这是候选基因分析而不是全基因组关联分析，这个阈值主要用于在已知的候选区域内进行信号的相对比较。

# 加载必要的包
library(ggplot2)

# 读取数据
assoc_data <- read.table("Zm00001d002439_DI_SX.txt", header=TRUE, sep="\t")

# 计算90%分位数阈值
threshold <- quantile(-log10(assoc_data$p), 0.9)

# 找出最显著的SNPs（前三个）
top_snps <- assoc_data[order(-log10(assoc_data$p), decreasing=TRUE)[1:3], ]

# 创建Manhattan图
p <- ggplot(assoc_data, aes(x=Pos, y=-log10(p))) +
  # 添加点
  geom_point(aes(color=-log10(p)), size=3, alpha=0.8) +
  # 使用从浅蓝到红色的渐变
  scale_color_gradientn(
    colors = c("#B8D7FF", "#FF9E9E", "#FF0000"),
    values = c(0, 0.5, 1),
    name = "-log10(p)"
  ) +
  # 添加90%分位数阈值线
  geom_hline(yintercept=threshold,
             linetype="dashed",
             color="red",
             alpha=0.5) +
  # 添加阈值标签
  annotate("text",
           x=min(assoc_data$Pos),
           y=threshold,
           label=paste0("90% quantile = ", round(threshold, 3)),
           hjust=-0.1,
           vjust=-0.5,
           color="red") +
  # 为最显著的SNPs添加位置标签
  geom_text(data=top_snps,
            aes(label=format(Pos, scientific=FALSE)),
            hjust=-0.1,
            vjust=0.5,
            size=3) +
  # 自定义主题
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size=12),
    axis.title = element_text(size=14),
    plot.title = element_text(hjust=0.5, size=16),
    plot.subtitle = element_text(hjust=0.5, size=12),
    legend.position = "right"
  ) +
  # 添加标签
  labs(
    x = "Position on Chromosome 2 (Mb)",
    y = "-log10(p-value)",
    title = "Local Association Analysis"
  ) +
  # 设置x轴刻度格式
  scale_x_continuous(labels = function(x) format(x/1e6, digits=3))

# 保存图片
ggsave("manhattan_plot.svg", p, width=12, height=6, dpi=300)

# 输出统计信息
cat("\n关键统计信息：\n")
cat("分析的SNP数量：", nrow(assoc_data), "\n")
cat("最显著的SNP：", top_snps$Marker[1], "\n")
cat("90%分位数阈值：", round(threshold, 3), "\n")

# 输出最显著的SNPs
cat("\n最显著的3个SNP：\n")
print(data.frame(
  Marker = top_snps$Marker,
  Position = top_snps$Pos,
  P_value = signif(top_snps$p, 3),
  Negative_log10P = signif(-log10(top_snps$p), 3)
))