当前位置：首页 > news >正文

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

news 2025/1/14 4:59:56

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

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

import pandas as pd
import numpy as npdef 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)
))