8000份大豆重测序数据文章的流程

8000份大豆重测序数据文章的流程

目录

  1. 分析流程总览
  2. 数据分析具体流程
  3. 图表制作

分析流程总览

整体流程架构

1
2
3
4
5
6
7
8
9
10
11
原始数据 (Raw Data)
    ↓
数据质控 (Quality Control)
    ↓
比对到参考基因组 (Alignment)
    ↓
变异检测 (Variant Detection)
    ↓
变异注释 (Variant Annotation)
    ↓
进化和群体遗传学分析 (Evolutionary and Population Genetics Analysis)

数据分析具体流程

一、数据质控

1. FastQC - 质量评估

软件功能：对测序数据进行质量评估

基本用法：

1
fastqc -t 8 -o /path/to/output /path/to/seq_data/*.fastq

参数说明：

  • -t 8：使用8个线程
  • -o：指定输出目录
  • 输入文件：fastq格式的测序数据

输出结果：

  • HTML格式的质量报告文件
  • 包含测序质量分布、GC含量、序列长度分布、接头含量等信息

2. fastp - 质控与过滤

软件功能：

  • 质量剪切（quality trimming）
  • 过滤低质量reads
  • 去接头
  • UMI（Unique Molecular Identifier）处理
  • 过滤polyG（Illumina NovaSeq特征）
  • 生成HTML和JSON格式的质控报告

基本用法：

1
2
3
fastp -i sample_R1.fq.gz -I sample_R2.fq.gz \
      -o clean_R1.fq.gz -O clean_R2.fq.gz \
      -w 8

常用参数详解：

参数 说明 默认值
-i 输入的R1（前向）reads文件 -
-I 输入的R2（反向）reads文件 -
-o 输出的清洁R1 reads文件 -
-O 输出的清洁R2 reads文件 -
-w 线程数 2
-q 质量值阈值，低于此值的碱基将被剪切 15
-u 低质量reads百分比阈值，超过此比例的reads将被过滤 30
-n reads长度阈值，短于此长度的reads将被过滤 0
-l 保留的最小reads长度 15
--detect_adapter_for_pe 自动检测接头（推荐） -
--trim_poly_g 去除polyG尾（NovaSeq数据必需） -
--poly_g_min_len 触发polyG剪切的 最小长度 10
--cut_front 从前端剪切指定数量的碱基 -
--cut_tail 从尾端剪切指定数量的碱基 -
--cut_mean_quality 按平均质量值剪切 -
--length_required reads最小长度要求 15
--json 输出JSON格式质控报告 fastp.json
--html 输出HTML格式质控报告 fastp.html
--report_title 报告标题 fastp report

二、序列比对

BWA-MEM 比对

软件功能：将清洁后的reads比对到参考基因组

完整流程：

2.1 建立参考基因组索引
1
bwa index reference.fasta

说明：此步骤会生成多个索引文件（.amb, .ann, .bwt, .pac, .sa）

2.2 序列比对
1
2
3
bwa mem -t 8 -M reference.fasta \
    clean_R1.fq.gz clean_R2.fq.gz \
    > sample.sam

参数说明：

  • -t 8：使用8个线程进行比对，提高速度
  • -M：将shorter split hits标记为次优，以便兼容Picard工具
  • reference.fasta：参考基因组序列文件
  • 输入：清洁后的R1和R2 reads（gzip压缩格式）
  • 输出：SAM格式比对结果
2.3 SAM转BAM并排序
1
2
3
4
5
6
7
8
# SAM转BAM（压缩格式）
samtools view -Sb sample.sam > sample.bam

# 对BAM文件进行排序
samtools sort sample.bam -o sample.sorted.bam

# 建立索引，便于后续快速访问
samtools index sample.sorted.bam
2.4 标记PCR重复
1
2
3
4
5
6
7
8
# 使用Picard标记重复
picard MarkDuplicates \
    I=sample.sorted.bam \
    O=sample.dedup.bam \
    M=metrics.txt

# 建立索引
samtools index sample.dedup.bam

说明：

  • I=：输入文件
  • O=：输出文件
  • M=：输出重复统计信息
2.5 添加/read groups信息（可选但推荐）
1
2
3
4
5
6
7
8
gatk AddOrReplaceReadGroups \
    -I sample.dedup.bam \
    -O sample.rg.bam \
    -RGID id1 \
    -RGLB lib1 \
    -RGPL illumina \
    -RGPU unit1 \
    -RGSM sample1

参数说明：

  • -RGID：Read Group ID
  • -RGLB：Library名称
  • -RGPL：测序平台（illumina, iontorrent等）
  • -RGPU：测序单元
  • -RGSM：样本名称

三、变异检测

GATK Best Practices 流程

GATK（Genome Analysis Toolkit）是目前最广泛使用的变异检测软件之一。

3.1 HaplotypeCaller - 单样本变异检测

基本用法（GVCF模式）：

1
2
3
4
5
gatk HaplotypeCaller \
    -R reference.fasta \
    -I sample.dedup.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF

参数说明：

  • -R：参考基因组fasta文件
  • -I：输入的BAM文件
  • -O：输出的gVCF文件（gzip压缩）
  • -ERC GVCF：产生GVCF格式，用于后续联合分析

说明：GVCF模式会记录所有位点的信息，包括变异位点和非变异位点，便于后续多样本联合分析。

3.2 CombineGVCFs - 合并多个gVCF
1
2
3
4
5
6
7
gatk CombineGVCFs \
    -R reference.fasta \
    --variant sample1.g.vcf.gz \
    --variant sample2.g.vcf.gz \
    --variant sample3.g.vcf.gz \
    ... \
    -O cohort.g.vcf.gz

说明：将多个样本的gVCF文件合并成一个联合gVCF文件

批处理示例：

1
2
3
4
5
6
7
8
9
# 创建variant列表文件
ls *.g.vcf.gz > variant_list.txt

# 合并
gatk CombineGVCFs \
    -R reference.fasta \
    --variant sample1.g.vcf.gz \
    --variant sample2.g.vcf.gz \
    -O cohort.g.vcf.gz
3.3 GenotypeGVCFs - 联合基因分型
1
2
3
4
gatk GenotypeGVCFs \
    -R reference.fasta \
    -V cohort.g.vcf.gz \
    -O cohort.vcf.gz

说明：对合并后的gVCF进行联合基因分型，生成最终的VCF文件

3.4 VariantFiltration - 变异过滤

SNP过滤标准：

1
2
3
4
5
6
7
8
9
gatk VariantFiltration \
    -R reference.fasta \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-name "QD_filter" --filter-expression "QD < 2.0" \
    --filter-name "FS_filter" --filter-expression "FS > 60.0" \
    --filter-name "MQ_filter" --filter-expression "MQ < 40.0" \
    --filter-name "MQRankSum_filter" --filter-expression "MQRankSum < -12.5" \
    --filter-name "ReadPosRankSum_filter" --filter-expression "ReadPosRankSum < -8.0"

InDel过滤标准：

1
2
3
4
5
6
7
gatk VariantFiltration \
    -R reference.fasta \
    -V cohort.vcf.gz \
    -O cohort.indel.filtered.vcf.gz \
    --filter-name "QD_filter" --filter-expression "QD < 2.0" \
    --filter-name "FS_filter" --filter-expression "FS > 200.0" \
    --filter-name "ReadPosRankSum_filter" --filter-expression "ReadPosRankSum < -20.0"

过滤指标说明：

指标 全称 说明 SNP阈值 InDel阈值
QD QualByDepth 每深度质量分数 < 2.0 < 2.0
FS FisherStrand 链偏向性 > 60.0 > 200.0
MQ RMSMappingQuality 根均方映射质量 < 40.0 -
MQRankSum MappingQualityRankSumTest 映射质量秩和检验 < -12.5 -
ReadPosRankSum ReadPositionRankSumTest Read位置秩和检验 < -8.0 < -20.0
SOR StrandOddsRatio 链比值 > 3.0 > 10.0
3.5 SelectVariants - 提取高质量变异

提取通过所有过滤的变异：

1
2
3
4
5
gatk SelectVariants \
    -R reference.fasta \
    -V cohort.filtered.vcf.gz \
    --exclude-filtered \
    -O cohort.high_quality.vcf.gz
3.6 分别提取SNP和InDel

提取SNP：

1
2
3
4
gatk SelectVariants \
    -V cohort.high_quality.vcf.gz \
    -select-type SNP \
    -O cohort.snp.vcf.gz

提取InDel：

1
2
3
4
gatk SelectVariants \
    -V cohort.high_quality.vcf.gz \
    -select-type INDEL \
    -O cohort.indel.vcf.gz

四、变异注释

ANNOVAR 注释流程

软件功能：对检测到的变异进行功能注释，预测变异对基因功能的影响

4.1 准备工作

下载注释数据库：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 建立数据库目录
mkdir humandb/

# 下载RefSeq基因数据库
annotate_variation.pl -buildver dm6 \
    -downdb -webfrom annovar \
    refGene humandb/

# 下载dbNSFP数据库（包含多个功能预测工具的结果）
annotate_variation.pl -buildver dm6 \
    -downdb -webfrom annovar \
    dbnsfp42a humandb/

# 其他常用数据库
annotate_variation.pl -buildver dm6 \
    -downdb -webfrom annovar \
    gnomad_genome humandb/
4.2 转换VCF格式为ANNOVAR输入格式
1
2
3
4
5
# SNP转换
convert2annovar.pl -format v4 cohort.snp.vcf.gz > cohort.snp.avinput

# InDel转换
convert2annovar.pl -format v4 cohort.indel.vcf.gz > cohort.indel.avinput
4.3 执行注释

SNP注释：

1
2
3
4
5
6
7
8
table_annovar.pl cohort.snp.avinput humandb/ \
    -buildver dm6 \
    -out cohort.snp_annotated \
    -remove \
    -protocol refGene,dbnsfp42a,gnomad_genome \
    -operation g,f,f \
    -nastring . \
    -vcfinput

InDel注释：

1
2
3
4
5
6
7
8
table_annovar.pl cohort.indel.avinput humandb/ \
    -buildver dm6 \
    -out cohort.indel_annotated \
    -remove \
    -protocol refGene \
    -operation g \
    -nastring . \
    -vcfinput

参数说明：

  • -buildver：基因组版本（如dm6, hg38等）
  • -protocol：指定要使用的注释数据库
  • -operation：
    • g：基于基因的注释
    • f：基于过滤的注释
  • -nastring：缺失值的表示字符串
  • -vcfinput：输入为VCF格式

注释结果解读：

  • exonic：外显子区的变异
  • splicing：剪接位点附近的变异（±2bp）
  • intronic：内含子区的变异
  • upstream/downstream：基因上游/下游区域
  • intergenic：基因间区
  • UTR5/UTR3：5’/3’非翻译区

功能影响预测：

  • synonymous SNV：同义突变（不改变氨基酸）
  • nonsynonymous SNV：非同义突变
    • missense：错义突变
    • nonsense：无义突变（产生终止密码子）
    • frameshift：移码突变

图表制作

一、群体遗传学分析

1.1 PCA分析（主成分分析）

分析目的：

  • 评估群体结构
  • 识别样本聚类模式
  • 检测离群样本
  • 揭示样本间的遗传关系

常用软件：

  • PLINK：命令行工具，快速高效
  • GCTA：支持混合线性模型
  • SNPRelate（R包）：灵活，便于可视化

PLINK进行PCA分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# Step 1: 将VCF转换为PLINK格式
plink --vcf cohort.vcf.gz --make-bed --out population

# Step 2: 数据过滤（可选）
plink --bfile population \
      --maf 0.05 \
      --geno 0.1 \
      --hwe 1e-6 \
      --make-bed \
      --out population_filtered

# Step 3: 进行PCA分析
plink --bfile population_filtered \
      --pca 20 \
      --out pca_result

# 输出文件：
# - pca_result.eigenval：特征值（解释方差比例）
# - pca_result.eigenvec：特征向量（样本的PC坐标）

R语言可视化：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 读取PCA结果
pca_data <- read.table("pca_result.eigenvec", header=FALSE)
colnames(pca_data) <- c("FID", "IID", paste0("PC", 1:20))

# 读取群体信息（如果有的话）
group_info <- read.table("group_info.txt", header=TRUE)
pca_data <- merge(pca_data, group_info, by="IID")

# 绘制PC1 vs PC2
library(ggplot2)
ggplot(pca_data, aes(x=PC1, y=PC2, color=Group)) +
    geom_point(size=3, alpha=0.7) +
    theme_bw() +
    labs(x="PC1", y="PC2",
         title="PCA Analysis",
         subtitle=paste0("PC1: ",
         round(summary(prcomp(pca_data[,3:22]))$importance[2,1]*100, 2),
         "% variance"))

# 绘制PC1 vs PC3
ggplot(pca_data, aes(x=PC1, y=PC3, color=Group)) +
    geom_point(size=3, alpha=0.7) +
    theme_bw() +
    labs(x="PC1", y="PC3")

# 绘制碎石图（Scree plot）
eigenvals <- read.table("pca_result.eigenval")
var_explained <- eigenvals$V1 / sum(eigenvals$V1) * 100
plot(1:20, var_explained[1:20], type="b",
     xlab="Principal Component",
     ylab="Variance Explained (%)",
     main="Scree Plot")

1.2 系统发育树构建

分析目的：

  • 展示样本间的进化关系
  • 识别单系群
  • 推断群体历史

常用方法：

方法 特点 适用场景 软件
邻接法(NJ) 速度快，适用于大数据 初步分析 MEGA, PHYLIP
最大似然法(ML) 准确性高，计算量大 精确分析 RAxML, IQ-TREE
贝叶斯法 最准确，计算量极大 小数据集精细分析 MrBayes

完整流程：

Step 1: 提取变异位点并格式转换

1
2
3
4
5
# 使用vcf2phylip（推荐）
vcf2phylip.py --input cohort.vcf.gz --output phylip_file

# 或使用seqtk
seqtk subseq reference.fasta region_list.txt > alignment.fasta

Step 2: 构建系统发育树（IQ-TREE）

1
2
3
4
5
6
7
8
9
10
11
12
13
# 基本分析
iqtree -s phylip_file -m GTR+G -bb 1000 -nt AUTO

# 参数说明：
# -s: 输入文件（PHYLIP或FASTA格式）
# -m: 替代模型（GTR+G常用，或使用MFP让软件自动选择）
# -bb: 自举重复次数（1000为标准）
# -nt: 线程数（AUTO自动检测）

# 输出文件：
# - *.treefile: 最佳系统发育树（Newick格式）
# - *.log: 运行日志
# - *.iqtree: 详细结果

Step 3: 使用RAxML（替代方案）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
raxmlHPC -s phylip_file \
         -n output \
         -m GTRGAMMA \
         -p 12345 \
         -# 100 \
         -f a

# 参数说明：
# -s: 输入文件
# -n: 输出文件后缀
# -m: 替代模型
# -p: 随机种子
# -#: 自举重复次数
# -f a: 快速自举分析和搜索最佳树

Step 4: 可视化

在线工具：

本地软件：

  • FigTree：图形界面软件
  • ggtree（R包）：可编程，灵活

R语言ggtree示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
library(ggtree)
library(treeio)

# 读取树文件
tree <- read.tree("output.treefile")

# 基本绘图
ggtree(tree) +
    geom_tiplab() +
    geom_tippoint(aes(color=group), size=3) +
    theme_tree2()

# 添加进化距离标尺
ggtree(tree) +
    geom_treescale() +
    geom_tiplab(size=3) +
    geom_nodepoint(aes(label=label), color="red", size=3)

# 添加自展值
ggtree(tree) +
    geom_nodelab(aes(label=bootstrap), size=3, color="blue")

1.3 群体结构分析

分析目的：

  • 评估群体混合程度
  • 估计祖先群体数量
  • 检测基因流事件

ADMIXTURE分析：

Step 1: 数据准备

1
2
3
4
5
6
7
8
9
10
11
12
# 转换为PLINK BED格式
plink --vcf cohort.vcf.gz --make-bed --out population

# 过滤相关样本（可选）
plink --bfile population \
      --indep-pairwise 50 5 0.2 \
      --out prune

plink --bfile population \
      --extract prune.prune.in \
      --make-bed \
      --out population.pruned

Step 2: 运行ADMIXTURE

1
2
3
4
5
6
7
8
# 测试不同的K值（祖先群体数量）
for K in {1..10}; do
    admixture --cv population.pruned.bed $K | tee log${K}.out
done

# 输出文件：
# - population.pruned.${K}.Q: 个体 ancestry proportions
# - population.pruned.${K}.P: SNP allele frequencies

Step 3: 选择最优K值

1
2
3
4
# 查看交叉验证误差
grep -h CV log*.out | sort -t':' -k2 -n

# 选择CV误差最小的K值

Step 4: 可视化（R语言）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
library(ggplot2)

# 读取Q矩阵（ancestry proportions）
K <- 3  # 假设最优K=3
Q <- read.table(paste0("population.pruned.", K, ".Q"))

# 设置群体顺序（按第一祖先成分排序）
order_idx <- order(Q[,1])
Q_ordered <- Q[order_idx,]

# 准备绘图数据
plot_data <- as.data.frame(Q_ordered)
colnames(plot_data) <- paste0("Ancestry", 1:K)
plot_data$Sample <- factor(1:nrow(Q_ordered))

# 转换为长格式
library(reshape2)
plot_data_long <- melt(plot_data, id.vars="Sample",
                       variable.name="Ancestry",
                       value.name="Proportion")

# 绘制堆叠条形图
ggplot(plot_data_long, aes(x=Sample, y=Proportion, fill=Ancestry)) +
    geom_bar(stat="identity", width=0.8) +
    scale_fill_brewer(palette="Set3") +
    theme_minimal() +
    theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position="right"
    ) +
    labs(x="Sample", y="Ancestry Proportion",
         title=paste0("ADMIXTURE Analysis (K=", K, ")"))

STRUCTURE分析（替代软件）：

1
2
3
4
5
# 运行STRUCTURE（需要图形界面或命令行）
structure -K 3 -i population_input.txt -o structure_output

# 使用STRUCTURE HARVESTER选择最优K
# 使用CLUMPAK或pophelper可视化

二、连锁不平衡分析

2.1 LD衰减分析

分析目的：

  • 评估连锁不平衡随着距离增加而衰减的程度
  • 估算历史重组率
  • 指导关联分析中的标记密度选择

关键指标：

指标 说明 取值范围
r² 两个位点间的相关系数 0-1
D’ 标准化连锁不平衡系数 0-1

LD衰减曲线：当r²降至其最大值的一半时，对应的物理距离即为LD衰减距离

软件1: PopLDdecay（推荐）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 基本用法
PopLDdecay -InVCF cohort.vcf.gz \
           -OutStat ld_decay_stat

# 参数说明：
# -InVCF: 输入VCF文件
# -OutStat: 输出统计文件

# 高级参数
PopLDdecay -InVCF cohort.vcf.gz \
           -OutStat ld_decay_stat \
           -MaxDist 1000000 \
           -MAF 0.05 \
           -Miss 0.2 \
           -Het 1

# -MaxDist: 最大距离（bp）
# -MAF: 最小等位基因频率
# -Miss: 最大缺失率
# -Het: 包含杂合位点

# 绘制LD衰减曲线
PopLDdecay -InVCF cohort.vcf.gz \
           -OutStat ld_decay_stat \
           -OutPlot ld_decay_plot

# 输出文件：
# - ld_decay_stat.stat: LD统计结果
# - ld_decay_stat.pdf: LD衰减曲线图

软件2: PLINK

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 计算配对LD
plink --vcf cohort.vcf.gz \
      --r2 \
      --ld-window 99999 \
      --ld-window-kb 100 \
      --ld-window-r2 0 \
      --out ld_pairwise

# 参数说明：
# --r2: 计算r²值
# --ld-window: 窗口内的最大SNP数量
# --ld-window-kb: 窗口大小（kb）
# --ld-window-r2: 只报告r²高于此值的配对

# 滑动窗口LD
plink --vcf cohort.vcf.gz \
      --ld-window 100 \
      --ld-window-kb 100 \
      --ld-window-r2 0.2 \
      --r2 --ld-window-r2 0 \
      --out ld_window

R语言可视化：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(ggplot2)

# 读取PopLDdecay结果
ld_data <- read.table("ld_decay_stat.stat", header=TRUE)
colnames(ld_data) <- c("Distance", "R2")

# 按距离分组求平均值
ld_avg <- aggregate(R2 ~ Distance, data=ld_data, FUN=mean)

# 绘制LD衰减曲线
ggplot(ld_avg, aes(x=Distance/1000, y=R2)) +
    geom_line(size=1, color="steelblue") +
    geom_smooth(method="loess", se=FALSE, color="red", linetype="dashed") +
    labs(x="Distance (kb)", y="Linkage Disequilibrium (r²)",
         title="LD Decay Analysis") +
    theme_bw() +
    annotate("segment", x=0, xend=Inf, y=0.1, yend=0.1,
             linetype="dotted", color="gray50")

# 计算半衰距离（r²降至0.1时的距离）
half_decay_distance <- approx(ld_avg$R2, ld_avg$Distance, xout=0.1)$y
print(paste0("Half-decay distance: ", round(half_decay_distance/1000, 2), " kb"))

2.2 LD热图

目的：展示特定染色体区域的LD模式

R语言实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
library(LDheatmap)
library(SNPRelate)

# 将VCF转换为GDS格式
snpgdsVCF2GDS("cohort.vcf.gz", "cohort.gds")

# 打开GDS文件
genofile <- snpgdsOpen("cohort.gds")

# 提取特定区域的SNP
snp_id <- snpgdsSelectSNP(genofile, chromosome=1,
                         start.position=1000000,
                         end.position=2000000)

# 计算LD矩阵
ld_matrix <- snpgdsLDMat(genofile, snp.id=snp_id,
                         method="r", slide=100,
                         LD.threshold=0.2)

# 绘制LD热图
png("LD_heatmap.png", width=2000, height=2000, res=150)
LDheatmap(ld_matrix$snpdist, ld_matrix$LDmatrix,
          flip=TRUE,
          title="LD Heatmap (Chr1: 1-2 Mb)",
          geneMap="no genes")
dev.off()

# 关闭GDS文件
snpgdsClose(genofile)

三、核苷酸多样性分析

3.1 π (Pi) 计算

定义：
π（核苷酸多样性）是指从群体中随机抽取的两个染色体的任意位点上，核苷酸差异的平均比例。

计算公式：

1
π = Σ πij / [n(n-1)/2]

其中，πij是第i个和第j个序列之间的核苷酸差异，n是样本数量

生物学意义：

  • 反映群体内的遗传多样性水平
  • 高π值：高遗传多样性
  • 低π值：低遗传多样性（可能经历了选择或瓶颈效应）

使用VCFtools计算：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 全基因组π（单点估计）
vcftools --vcf cohort.vcf.gz \
         --site-pi \
         --out pi_per_site

# 滑动窗口π
vcftools --vcf cohort.vcf.gz \
         --window-pi 100000 \
         --window-pi-step 50000 \
         --out pi_window

# 参数说明：
# --window-pi: 窗口大小（bp）
# --window-pi-step: 滑动步长（bp）

# 输出文件：
# - pi_per_site.sites.pi: 每个位点的π值
# - pi_window.windowed.pi: 窗口π值

使用PopGenome（R包）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
library(PopGenome)

# 读取VCF文件
data <- readVCF("cohort.vcf.gz", include.unknown=TRUE)

# 计算全基因组核苷酸多样性
data <- diversity.stats(data)

# 提取π值
pi_genome <- data@nuc.diversity.within
print(paste0("Genome-wide π: ", pi_genome))

# 滑动窗口分析
data <- sliding.window.transform(data, width=100000, step=50000)
data <- diversity.stats(data)

# 提取窗口π值
window_pi <- data@nuc.diversity.within
window_positions <- data@sliding.window.positions

# 可视化
plot(window_positions/1000000, window_pi,
     type="l", lwd=2,
     xlab="Genome Position (Mb)",
     ylab="Nucleotide Diversity (π)",
     main="Sliding Window π Analysis")

群体间π比较：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 群体1
vcftools --vcf cohort.vcf.gz \
         --keep population1.txt \
         --window-pi 100000 \
         --out pi_pop1

# 群体2
vcftools --vcf cohort.vcf.gz \
         --keep population2.txt \
         --window-pi 100000 \
         --out pi_pop2

# 群体3
vcftools --vcf cohort.vcf.gz \
         --keep population3.txt \
         --window-pi 100000 \
         --out pi_pop3

3.2 Watterson’s θ (Theta)

定义：
基于分离位点（segregating sites）数目的核苷酸多样性估计量

计算公式：

1
θ = S / an

其中：

  • S：分离位点数目
  • an：调和级数和，an = Σ[1/(i-1)], i从2到n

生物学意义：

  • 与π类似，也是衡量遗传多样性的指标
  • 在中性进化条件下，θ ≈ π
  • π < θ：说明群体经历了扩张
  • π > θ：说明存在平衡选择或群体收缩

计算方法：

1
2
3
4
5
6
7
8
9
10
11
12
13
library(PopGenome)

# 读取数据
data <- readVCF("cohort.vcf.gz", include.unknown=TRUE)

# 计算Watterson's θ
data <- neutrality.stats(data)

# 提取结果
theta <- data@Watterson.theta
theta <- data@Tajimas.D

print(paste0("Watterson's θ: ", theta))

Tajima’s D检验：

1
2
3
4
5
6
7
# 使用VCFtools
vcftools --vcf cohort.vcf.gz \
         --TajimaD 100000 \
         --out tajimaD

# 输出文件：
# - tajimaD.Tajima.D: 窗口Tajima's D值

四、固定指数 (FST)

4.1 FST计算

定义：
FST（Fixation Index）是衡量群体间遗传分化程度的指标

取值范围及意义：

FST值 分化程度 解释
0 - 0.05 很小 群体间几乎无分化
0.05 - 0.15 中等 群体间有一定分化
0.15 - 0.25 较大 群体间分化显著
> 0.25 极大 群体间分化极大

使用VCFtools计算：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 准备群体文件（每行一个样本ID）
# population1.txt
Sample001
Sample002
...

# population2.txt
Sample101
Sample102
...

# 计算两群体间FST（单点）
vcftools --vcf cohort.vcf.gz \
         --weir-fst-pop population1.txt \
         --weir-fst-pop population2.txt \
         --out fst_pop1_pop2

# 滑动窗口FST
vcftools --vcf cohort.vcf.gz \
         --weir-fst-pop population1.txt \
         --weir-fst-pop population2.txt \
         --fst-window-size 100000 \
         --fst-window-step 50000 \
         --out fst_window

# 参数说明：
# --weir-fst-pop: 指定一个群体的样本文件
# --fst-window-size: 窗口大小（bp）
# --fst-window-step: 滑动步长（bp）

使用PopGenome（R包）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
library(PopGenome)

# 读取数据
data <- readVCF("cohort.vcf.gz", include.unknown=TRUE,
                gffpath="annotation.gff")

# 定义群体
populations <- list(
    pop1 = read.table("population1.txt")$V1,
    pop2 = read.table("population2.txt")$V1,
    pop3 = read.table("population3.txt")$V1
)

# 设置群体
data <- set.populations(data, populations)

# 计算FST
data <- F_ST.stats(data, mode="nucleotide")

# 提取结果
fst_matrix <- data@F_ST.matrix
print(fst_matrix)

# 成对群体FST
for (i in 1:(nrow(fst_matrix)-1)) {
    for (j in (i+1):nrow(fst_matrix)) {
        cat(sprintf("FST(%s, %s) = %.4f\n",
                    rownames(fst_matrix)[i],
                    colnames(fst_matrix)[j],
                    fst_matrix[i,j]))
    }
}

可视化FST分布：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(ggplot2)

# 读取窗口FST结果
fst_data <- read.table("fst_window.windowed.weir.fst", header=TRUE)

# 绘制曼哈顿式FST图
ggplot(fst_data, aes(x=POS/1000000, y=WEIR_AND_COCKERHAM_FST)) +
    geom_point(size=0.5, alpha=0.6, color="steelblue") +
    geom_hline(yintercept=0.15, linetype="dashed", color="red") +
    facet_wrap(~CHROM, scales="free_x", ncol=4) +
    labs(x="Position (Mb)", y="FST",
         title="Genome-wide FST between Population 1 and 2") +
    theme_bw() +
    theme(strip.text = element_text(size=8))

FST热图：

1
2
3
4
5
6
7
8
9
library(pheatmap)

# 绘制FST矩阵热图
pheatmap(fst_matrix,
         display_numbers=TRUE,
         number_format="%.3f",
         main="Pairwise FST between Populations",
         color=colorRampPalette(c("white", "orange", "red"))(50),
         breaks=seq(0, max(fst_matrix), length.out=51))

五、选择清除分析

5.1 核苷酸多样性比值 (π-ratio)

原理：
受选择区域的核苷酸多样性会显著降低，通过计算不同群体间π的比值可以检测选择信号

计算公式：

1
π-ratio = π_selected / π_control

低π-ratio值：该区域可能经历了正选择

分析流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Step 1: 计算各群体π
# 群体1（假设为受选择群体）
vcftools --vcf cohort.vcf.gz \
         --keep population1.txt \
         --window-pi 100000 \
         --window-pi-step 50000 \
         --out pi_pop1

# 群体2（对照群体）
vcftools --vcf cohort.vcf.gz \
         --keep population2.txt \
         --window-pi 100000 \
         --window-pi-step 50000 \
         --out pi_pop2

# Step 2: 计算π-ratio（使用Python或R）
# Python示例
import pandas as pd

pi1 = pd.read_csv("pi_pop1.windowed.pi", sep="\t")
pi2 = pd.read_csv("pi_pop2.windowed.pi", sep="\t")

pi_ratio = pi1.copy()
pi_ratio["PI_RATIO"] = pi1["PI"] / pi2["PI"]
pi_ratio.to_csv("pi_ratio.txt", sep="\t", index=False)

R语言可视化：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
library(ggplot2)

# 读取数据
pi1 <- read.table("pi_pop1.windowed.pi", header=TRUE)
pi2 <- read.table("pi_pop2.windowed.pi", header=TRUE)

# 计算π-ratio
pi_ratio <- data.frame(
    CHROM = pi1$CHROM,
    BIN_START = pi1$BIN_START,
    BIN_END = pi1$BIN_END,
    PI_RATIO = pi1$PI / pi2$PI
)

# 识别选择信号（阈值：5%分位数）
threshold <- quantile(pi_ratio$PI_RATIO, 0.05)
candidates <- pi_ratio[pi_ratio$PI_RATIO < threshold, ]

# 可视化
ggplot(pi_ratio, aes(x=BIN_START/1000000, y=PI_RATIO)) +
    geom_line(color="steelblue", size=0.5) +
    geom_hline(yintercept=threshold, linetype="dashed",
               color="red", size=1) +
    geom_point(data=candidates, aes(x=BIN_START/1000000, y=PI_RATIO),
               color="red", size=2, alpha=0.7) +
    facet_wrap(~CHROM, scales="free_x", ncol=4) +
    labs(x="Position (Mb)", y="π Ratio (Pop1/Pop2)",
         title="Nucleotide Diversity Ratio Scan") +
    theme_bw()

# 输出候选区域
write.table(candidates, "selection_candidates.txt",
            sep="\t", row.names=FALSE, quote=FALSE)

5.2 XP-CLR (跨群体复合似然比检验)

原理：
检测由于选择引起的群体间等位基因频率差异，考虑了连锁不平衡和遗传重组

特点：

  • 不需要单体型相位信息
  • 适用于大数据集
  • 对古老选择信号敏感

分析流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Step 1: 准备输入文件
# XP-CLR要求特定的格式：Chrom  Position  Likelihood

# 将VCF转换为XP-CLR格式
# 可以使用自定义脚本或vcftools

# 提取多等位基因信息
vcftools --vcf cohort.vcf.gz \
         --freq \
         --out allele_freq

# Step 2: 运行XP-CLR
# 安装XP-CLR后
xpclr -xpclr cohort.vcf.gz \
       -w1 1 200 2000 \
       -w2 1 100 1000 \
       -p0 0.001 \
       -maxSNP 1000 \
       -out xpclr_result

# 参数说明：
# -w1: 窗口参数（大小，步长，最大SNP数）
# -w2: 第二个窗口参数
# -p0: 置信度阈值
# -maxSNP: 每个窗口的最大SNP数

结果可视化：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
library(ggplot2)

# 读取XP-CLR结果
xpclr <- read.table("xpclr_result", header=FALSE)
colnames(xpclr) <- c("CHROM", "Position", "XP_CLR_Score",
                     "Ref_Freq", "Alt_Freq")

# 计算滑动窗口平均
window_size <- 100000  # 100 kb
xpclr$Window <- floor(xpclr$Position / window_size) * window_size
xpclr_avg <- aggregate(XP_CLR_Score ~ Window + CHROM,
                       data=xpclr, FUN=mean)

# 绘制XP-CLR扫描图
ggplot(xpclr_avg, aes(x=Window/1000000, y=XP_CLR_Score)) +
    geom_line(color="steelblue", size=1) +
    facet_wrap(~CHROM, scales="free_x", ncol=4) +
    labs(x="Position (Mb)", y="XP-CLR Score",
         title="XP-CLR Scan for Selection Signals") +
    theme_bw()

# 识别显著峰（阈值：top 1%）
threshold <- quantile(xpclr_avg$XP_CLR_Score, 0.99)
peaks <- xpclr_avg[xpclr_avg$XP_CLR_Score > threshold, ]

# 标注显著峰
ggplot(xpclr_avg, aes(x=Window/1000000, y=XP_CLR_Score)) +
    geom_line(color="gray") +
    geom_point(data=peaks, aes(x=Window/1000000, y=XP_CLR_Score),
               color="red", size=2) +
    geom_hline(yintercept=threshold, linetype="dashed",
               color="red") +
    facet_wrap(~CHROM, scales="free_x", ncol=4) +
    labs(x="Position (Mb)", y="XP-CLR Score",
         title="Significant Selection Peaks") +
    theme_bw()

5.3 XP-EHH (交叉群体扩展单体型纯合性)

原理：
检测一个群体中已经完成的正选择（单体型频率高且扩展范围大），而另一群体中未完成的选择

特点：

  • 适用于检测近期的选择事件
  • 需要单体型相位信息

分析流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Step 1: 准备PLINK格式数据
plink --vcf cohort.vcf.gz --make-bed --out cohort

# Step 2: 相位推断（使用SHAPEIT或BEAGLE）
shapeit -B cohort \
        -M genetic_map.txt \
        -O cohort.phased \
        --thread 8

# 或使用beagle
beagle gt=cohort.bed \
       out=cohort.phased \
       nthreads=8

# Step 3: 运行XP-EHH（使用selscan）
selscan --xpehh \
        --vcf cohort.phased.vcf.gz \
        --map genetic_map.txt \
        --pop1 population1.txt \
        --pop2 population2.txt \
        --out xpehh_result

# 参数说明：
# --xpehh: 运行XP-EHH模式
# --pop1: 群体1的样本列表
# --pop2: 群体2的样本列表

# Step 4: 标准化XP-EHH得分
norm --xpehh --files xpehh_result

结果可视化：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(ggplot2)

# 读取标准化后的XP-EHH结果
xpehh <- read.table("xpehh_result.norm.xpehh", header=TRUE)

# 绘制XP-EHH扫描图
ggplot(xpehh, aes(x=POS/1000000, y=normXPEHH)) +
    geom_point(size=0.5, alpha=0.6, color="steelblue") +
    geom_hline(yintercept=0, linetype="dashed", color="gray50") +
    geom_hline(yintercept=2, linetype="dashed", color="red") +
    geom_hline(yintercept=-2, linetype="dashed", color="blue") +
    facet_wrap(~CHR, scales="free_x", ncol=4) +
    labs(x="Position (Mb)", y="Standardized XP-EHH",
         title="XP-EHH Scan") +
    theme_bw()

# 识别显著信号（|XP-EHH| > 2）
candidates_pos <- xpehh[xpehh$normXPEHH > 2, ]
candidates_neg <- xpehh[xpehh$normXPEHH < -2, ]

5.4 Tajima’s D

定义：
基于等位基因频谱的中性检验统计量

计算公式：

1
D = (π - θ) / SD(π - θ)

取值及生物学意义：

D值 意义 进化解释
D >> 0 显著为正 群体收缩或平衡选择
D ≈ 0 接近0 符合中性进化
D << 0 显著为负 群体扩张或正选择

使用VCFtools计算：

1
2
3
4
5
6
7
8
9
10
# 滑动窗口Tajima's D
vcftools --vcf cohort.vcf.gz \
         --TajimaD 100000 \
         --out tajimaD

# 参数说明：
# --TajimaD: 窗口大小（bp）

# 输出文件：
# - tajimaD.Tajima.D: 窗口Tajima's D值

使用PopGenome（R包）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
library(PopGenome)

# 读取数据
data <- readVCF("cohort.vcf.gz", include.unknown=TRUE)

# 滑动窗口分析
data <- sliding.window.transform(data, width=100000, step=50000)

# 计算Tajima's D
data <- neutrality.stats(data)

# 提取结果
tajima_d <- data@Tajimas.D
window_positions <- data@sliding.window.positions

# 可视化
plot(window_positions/1000000, tajima_d,
     type="l", lwd=2, col="steelblue",
     xlab="Genome Position (Mb)",
     ylab="Tajima's D",
     main="Tajima's D Scan")

# 添加阈值线（显著性水平）
abline(h=2, col="red", lty=2)  # 显著正D值
abline(h=-2, col="blue", lty=2)  # 显著负D值

# 识别显著窗口
sig_pos <- which(tajima_d > 2)
sig_neg <- which(tajima_d < -2)

六、全基因组关联分析 (GWAS)

6.1 分析流程总览

完整流程：

1
表型数据收集 → 基因型数据质控 → 群体结构控制 → 关联分析 → 多重检验校正 → 结果可视化 → 候选基因鉴定

6.2 数据准备

表型文件：

1
2
3
4
5
Sample_ID    Trait1    Trait2    Trait3
Sample001    10.5      25.3      1.2
Sample002    12.3      28.1      1.5
Sample003    9.8       24.7      1.1
...

基因型文件：

  • VCF格式
  • PLINK BED格式
  • HapMap格式

群体结构文件（可选但推荐）：

  • Q矩阵：从ADMIXTURE或STRUCTURE获得
  • PC矩阵：从PCA获得

6.3 基因型数据质控

质控标准：

参数 阈值 说明
MAF > 0.05 最小等位基因频率
Missing rate < 0.1 缺失率（位点）
HWE P > 1e-6 哈代-温伯格平衡p值
Geno missing < 0.1 缺失率（个体）

使用PLINK进行质控：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# Step 1: 转换格式
plink --vcf cohort.vcf.gz --make-bed --out gwas_raw

# Step 2: 基本质控
plink --bfile gwas_raw \
      --maf 0.05 \
      --geno 0.1 \
      --hwe 1e-6 \
      --mind 0.1 \
      --make-bed \
      --out gwas_cleaned

# Step 3: 去除相关样本（可选）
plink --bfile gwas_cleaned \
      --indep-pairwise 50 5 0.2 \
      --out prune

plink --bfile gwas_cleaned \
      --extract prune.prune.in \
      --make-bed \
      --out gwas_pruned

# 参数说明：
# --maf: 最小等位基因频率
# --geno: 位点最大缺失率
# --hwe: 哈代-温伯格平衡p值阈值
# --mind: 个体最大缺失率
# --indep-pairwise: LD剪裁（窗口大小，步长，r²阈值）

6.4 关联分析

常用软件对比：

软件 模型 特点 适用场景
TASSEL GLM/MLM 图形界面，易于使用 小数据集
GAPIT 多种模型 R包，灵活 中等数据集
GEMMA LMM 快速，准确 大数据集
EMMAX LMM 高效 大数据集

使用GEMMA进行关联分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# Step 1: 计算关系矩阵（Kinship Matrix）
gemma -bfile gwas_cleaned \
      -gk 1 \
      -o kinship_matrix

# -gk 1: 中心化关系矩阵（VanRaden method）
# 输出：kinship_matrix.cXX.txt

# Step 2: 准备表型文件
# 格式：Sample_ID  Phenotype
# （需要与PLINK的.fam文件顺序一致）

# Step 3: 关联分析（线性混合模型）
gemma -bfile gwas_cleaned \
      -k output/kinship_matrix.cXX.txt \
      -lmm 4 \
      -phenotype phenotype.txt \
      -o gwas_result

# 参数说明：
# -k: 关系矩阵文件
# -lmm: 混合模型方法（1=exact LRT, 4= Wald test）
# -phenotype: 表型文件

# 输出文件：
# - gwas_result.assoc.txt: 关联结果

使用GAPIT（R包）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
library(GAPIT3)

# 读取数据
myG <- read.table("gwas_cleaned.raw", header=TRUE)
myY <- read.table("phenotype.txt", header=TRUE)

# PCA分析（作为协变量）
myPCA <- GAPIT.PCA(myG)

# 运行GWAS
myGWAS <- GAPIT(
    Y = myY,
    G = myG,
    PCA.total = 3,  # 使用前3个PC作为协变量
    model = "MLM",  # 混合线性模型
    SNP.MAF = 0.05
)

# 多模型比较
myGWAS_multi <- GAPIT(
    Y = myY,
    G = myG,
    PCA.total = 3,
    model = c("GLM", "MLM", "FarmCPU", "Blink")
)

# 输出：
# - myGWAS: Manhattan图和QQ图
# - GAPIT.GWAS.Result.txt: 关联结果

6.5 结果解读

关联结果文件格式：

列名 说明
CHR 染色体编号
SNP SNP标识符（位置或ID）
BP 物理位置
P p值
BETA 效应大小
SE 标准误
R2 解释方差比例

多重检验校正：

方法 阈值 说明
Bonferroni 0.05/N 保守，适合严格标准
FDR 0.05 允许5%假阳性
Suggestive 1e-5 建议显著性阈值

基因组显著性阈值：

1
2
3
4
5
Bonferroni correction: 0.05 / N
其中 N = 总SNP数量

例如，N = 1,000,000：
阈值 = 0.05 / 1,000,000 = 5e-8

6.6 结果可视化

曼哈顿图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
library(qqman)

# 读取结果
gwas <- read.table("gwas_result.assoc.txt", header=TRUE)

# 基本曼哈顿图
manhattan(gwas,
          chr = "CHR",
          bp = "BP",
          p = "P",
          genomewideline = -log10(5e-8),
          suggestiveline = -log10(1e-5),
          highlight = NULL,
          main = "GWAS Manhattan Plot")

# 高亮特定SNP
manhattan(gwas,
          chr = "CHR",
          bp = "BP",
          p = "P",
          highlight = gwas$SNP[gwas$P < 5e-8],
          main = "Significant SNPs")

# 添加基因名称
manhattan(gwas,
          chr = "CHR",
          bp = "BP",
          p = "P",
          annotatePval = 1e-6,  # 标注p值小于1e-6的SNP
          annotateTop = TRUE)   # 只标注最显著的SNP

使用ggplot2绘制曼哈顿图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
library(ggplot2)
library(dplyr)

# 准备数据
gwas_plot <- gwas %>%
    mutate(
        neglog10p = -log10(P),
        # 为每个染色体分配颜色
        chr_color = ifelse(as.numeric(CHR) %% 2 == 0, "even", "odd"),
        # 转换染色体为因子
        CHR = factor(CHR, levels=unique(CHR))
    )

# 计算阈值
threshold <- 5e-8

# 绘图
ggplot(gwas_plot, aes(x=BP/1000000, y=neglog10p)) +
    geom_point(aes(color=chr_color), size=0.5, alpha=0.7) +
    facet_wrap(~CHR, scales="free_x", ncol=4) +
    scale_color_manual(values=c("odd"="steelblue", "even"="orange")) +
    geom_hline(yintercept=-log10(threshold),
               linetype="dashed", color="red", size=1) +
    labs(x="Position (Mb)", y="-log10(P)",
         title="Genome-wide Association Study",
         subtitle=paste0("Significance threshold: ", threshold)) +
    theme_bw() +
    theme(
        legend.position="none",
        strip.text = element_text(size=8),
        axis.text.x = element_text(angle=45, hjust=1)
    )

QQ图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 使用qqman包
qq(gwas$P, main="Q-Q Plot")

# 使用ggplot2
library(ggplot2)

# 计算期望p值
expected <- -log10((1:nrow(gwas)) / (nrow(gwas) + 1))
observed <- -log10(sort(gwas$P))

qq_data <- data.frame(Expected=expected, Observed=observed)

ggplot(qq_data, aes(x=Expected, y=Observed)) +
    geom_point(size=0.5, alpha=0.6) +
    geom_abline(intercept=0, slope=1, linetype="dashed",
                color="red", size=1) +
    labs(x="Expected -log10(P)",
         y="Observed -log10(P)",
         title="Q-Q Plot") +
    theme_bw()

# 计算膨胀因子
lambda <- median(qchisq(1 - gwas$P, 1)) / qchisq(0.5, 1)
cat(paste0("Genomic inflation factor (λ): ", round(lambda, 3)))

火山图：

1
2
3
4
5
6
7
8
9
ggplot(gwas, aes(x=BETA, y=-log10(P))) +
    geom_point(aes(color=-log10(P) > -log10(5e-8)),
               alpha=0.6, size=0.5) +
    scale_color_manual(values=c("FALSE"="gray", "TRUE"="red")) +
    geom_hline(yintercept=-log10(5e-8), linetype="dashed") +
    labs(x="Effect Size (BETA)",
         y="-log10(P)",
         title="Volcano Plot") +
    theme_bw()

区域关联图（Locus Zoom）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# 选择显著信号区域
sig_snp <- gwas %>%
    filter(P < 5e-8) %>%
    arrange(P) %>%
    slice(1)  # 最显著的SNP

chr <- as.character(sig_snp$CHR)
center_pos <- sig_snp$BP
window <- 500000  # 上下各500kb

# 提取区域数据
region_data <- gwas %>%
    filter(CHR == chr & BP >= (center_pos - window) &
           BP <= (center_pos + window))

# 绘制区域图
ggplot(region_data, aes(x=BP/1000000, y=-log10(P))) +
    geom_point(color="steelblue", size=1.5) +
    geom_vline(xintercept=center_pos/1000000,
               linetype="dashed", color="red") +
    labs(x=paste0("Position on Chr", chr, " (Mb)"),
         y="-log10(P)",
         title=paste0("Regional Association Plot (", round(center_pos/1000000, 2), " Mb)")) +
    theme_bw()

# 使用LocusZoom在线工具
# https://locuszoom.org/

6.7 候选基因鉴定

提取显著SNP附近的基因：

1
2
3
4
5
6
7
8
9
# 提取显著SNP
awk 'NR>1 && $9 < 5e-8' gwas_result.assoc.txt > significant_snps.txt

# 根据基因组注释文件提取附近基因
# 假设有GFF3格式的注释文件
bedtools closest -a significant_snps.bed \
                -b genes.gff3 \
                -d 100000 \
                > candidate_genes.txt

使用BEDTools：

1
2
3
4
5
# 将SNP转换为BED格式
awk 'BEGIN{OFS="\t"} {print $1, $2-1, $2, $3, $4}' significant_snps.txt > snps.bed

# 提取SNP附近的基因（上下各100kb）
bedtools window -a snps.bed -b genes.bed -w 100000 > snp_genes.txt

功能富集分析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 使用clusterProfiler进行GO和KEGG富集分析
library(clusterProfiler)
library(org.At.tair.db)  # 拟南芥注释（大豆需相应数据库）

# 读取候选基因列表
genes <- read.table("candidate_genes.txt", header=FALSE)$V1

# GO富集分析
go_enrich <- enrichGO(
    gene = genes,
    OrgDb = org.At.tair.db,
    keyType = "TAIR",
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.10
)

# 绘图
dotplot(go_enrich, showCategory=20) + ggtitle("GO Enrichment")
barplot(go_enrich, showCategory=20) + ggtitle("GO Enrichment")
emapplot(go_enrich) + ggtitle("GO Enrichment")

# KEGG通路分析
kegg_enrich <- enrichKEGG(
    gene = genes,
    organism = 'ath',  # 拟南芥
    pvalueCutoff = 0.05
)

dotplot(kegg_enrich) + ggtitle("KEGG Pathway Analysis")

七、其他分析

7.1 基因家族分析

目的：鉴定特定基因家族的成员、进化历史和表达模式

分析流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# Step 1: 下载目标基因家族的已知蛋白序列
# 从Pfam、InterPro或文献获取

# Step 2: 在大豆基因组中搜索同源基因
# 使用BLASTP或HMMER

blastp -query known_family_proteins.fa \
       -db soybean_proteome.fa \
       -evalue 1e-5 \
       -outfmt 6 \
       -out blast_results.txt

# 或使用HMMER（更准确）
hmmbuild gene_family.hmm known_family_alignment.sto
hmmsearch --domtblout hmm_results.txt gene_family.hmm soybean_proteome.fa

# Step 3: 提取候选基因序列
samtools faidx soybean_proteome.fa $(awk '{print $1}' blast_results.txt) > candidate_genes.fa

# Step 4: 多序列比对
mafft --auto candidate_genes.fa > candidate_genes_aligned.fa

# Step 5: 构建系统发育树
iqtree -s candidate_genes_aligned.fa -m MFP -bb 1000 -nt AUTO

# Step 6: 计算Ka/Ks比值
# 使用KaKs_Calculator
KaKs_Calculator -i coding_sequences.fa -o ka_ks_results.txt

# Step 7: 染色体定位
# 使用基因GFF文件提取位置信息
grep -f candidate_genes_list.txt soybean_annotation.gff3 > gene_locations.gff3

可视化基因家族：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(ggtree)
library(ggplot2)

# 读取树文件
tree <- read.tree("candidate_genes_aligned.fa.treefile")

# 绘制系统发育树
ggtree(tree) +
    geom_tiplab(size=2) +
    geom_tippoint(aes(color=group), size=2) +
    theme_tree2() +
    labs(title="Gene Family Phylogenetic Tree")

# 添加热图（表达数据）
if (file.exists("expression_data.txt")) {
    expr_data <- read.table("expression_data.txt", header=TRUE, row.names=1)
    gheatmap(tree, expr_data, offset=0.5, width=0.5,
             colnames_angle=-45, colnames_offset_y=5)
}

7.2 选择压力分析

Ka/Ks比值分析：

定义：

  • Ka：非同义替换率
  • Ks：同义替换率
  • Ka/Ks：选择压力指标

比值解读：

Ka/Ks 选择类型 说明
Ka/Ks < 1 纯化选择 功能约束，去除有害突变
Ka/Ks = 1 中性进化 无选择压力
Ka/Ks > 1 正选择 适应性进化

使用KaKs_Calculator：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 准备输入文件（FASTA格式的编码序列对）
# 格式：>gene1_seq1
#       ATG...
#       >gene1_seq2
#       ATG...

# 运行KaKs_Calculator
KaKs_Calculator -i gene_pairs.fa -o ka_ks_results.txt -m MYN

# 参数说明：
# -i: 输入文件
# -o: 输出文件
# -m: 方法（MYN, NG, LWL, MLP等）

# 输出包含：
# - Ka, Ks, Ka/Ks
# - 标准误和置信区间

使用PAML (codeml)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 准备必要的文件：
# 1. 序列比对文件（PHYLIP格式）
# 2. 树文件（Newick格式）
# 3. 控制文件（codeml.ctl）

# codeml.ctl 示例
seqfile = gene_alignment.phy
treefile = gene_tree.nwk
outfile = codeml_results.txt
runmode = 0
model = 0
NSsites = 0
icode = 0

# 运行codeml
codeml codeml.ctl

可视化Ka/Ks分布：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(ggplot2)

# 读取结果
ka_ks <- read.table("ka_ks_results.txt", header=TRUE)

# 绘制Ka/Ks分布直方图
ggplot(ka_ks, aes(x=Ka_Ks)) +
    geom_histogram(binwidth=0.1, fill="steelblue", color="black") +
    geom_vline(xintercept=1, linetype="dashed", color="red", size=1) +
    labs(x="Ka/Ks Ratio", y="Count",
         title="Distribution of Ka/Ks Ratios") +
    theme_bw()

# Ka vs Ks散点图
ggplot(ka_ks, aes(x=Ks, y=Ka)) +
    geom_point(alpha=0.6) +
    geom_abline(intercept=0, slope=1, linetype="dashed",
                color="red", size=1) +
    labs(x="Ks (synonymous)",
         y="Ka (nonsynonymous)",
         title="Ka vs Ks Plot") +
    theme_bw()

7.3 单倍型分析

目的：分析基因区域的单倍型组成和频率分布

使用Haploview：

1
2
3
4
5
6
7
8
# 导入数据
# File > Import Data > PLINK format

# 自动定义单倍型块
# Blocks > Define Blocks

# 输出LD图和单倍型频率
# File > Export Plots

使用PLINK：

1
2
3
4
5
6
7
8
9
10
11
12
13
# 计算单倍型频率
plink --bfile gwas_cleaned \
      --blocks no-pheno-req \
      --hap-freq \
      --out haplotype_freq

# 参数说明：
# --blocks: 定义单倍型块
# --hap-freq: 输出单倍型频率

# 输出文件：
# - haplotype_freq.frq: 单倍型频率
# - haplotype_freq.blocks: 单倍型块定义

使用VCFtools：

1
2
3
4
5
6
7
8
9
# 提取特定区域
vcftools --vcf cohort.vcf.gz \
         --chr 1 \
         --from-bp 1000000 \
         --to-bp 2000000 \
         --recode \
         --out region

# 使用Haploview或其他软件分析单倍型

7.4 受选择基因功能富集

目的：确定受选择基因显著富集的功能类别

分析流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# 提取受选择区域的基因
selected_genes <- read.table("selected_genes.txt", header=FALSE)$V1

# GO富集分析
library(clusterProfiler)
library(org.At.tair.db)

go_bp <- enrichGO(
    gene = selected_genes,
    OrgDb = org.At.tair.db,
    keyType = "TAIR",
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.01,
    qvalueCutoff = 0.05
)

# 可视化
dotplot(go_bp, showCategory=25) +
    ggtitle("GO Biological Process Enrichment")

cnetplot(go_bp, showCategory=15) +
    ggtitle("Gene Concept Network")

goplot(go_bp) +
    ggtitle("GO Relationship")

# KEGG富集
kegg <- enrichKEGG(
    gene = selected_genes,
    organism = 'atha',
    pvalueCutoff = 0.05
)

dotplot(kegg, showCategory=20) +
    ggtitle("KEGG Pathway Enrichment")

# 比较不同群体受选择基因
# 群体1受选择基因
genes_pop1 <- read.table("selected_genes_pop1.txt")$V1

# 群体2受选择基因
genes_pop2 <- read.table("selected_genes_pop2.txt")$V1

go_pop1 <- enrichGO(gene = genes_pop1, OrgDb = org.At.tair.db,
                    ont = "BP", pAdjustMethod = "BH")
go_pop2 <- enrichGO(gene = genes_pop2, OrgDb = org.At.tair.db,
                    ont = "BP", pAdjustMethod = "BH")

# 比较富集结果
compare <- compareCluster(gene = list(Population1 = genes_pop1,
                                      Population2 = genes_pop2),
                          fun = "enrichGO",
                          OrgDb = org.At.tair.db,
                          ont = "BP")

dotplot(compare, showCategory=20) +
    ggtitle("GO Enrichment Comparison between Populations")

总结

本流程文档详细描述了8000份大豆重测序数据的完整分析pipeline，涵盖了从原始数据到最终发表图表的所有关键步骤：

核心流程

  1. 数据质控：FastQC和fastp确保数据质量
  2. 序列比对：BWA-MEM高效比对到参考基因组
  3. 变异检测：GATK Best Practices流程检测SNP和InDel
  4. 变异注释：ANNOVAR进行功能注释和效应预测

高级分析

  1. 群体遗传学分析

    • PCA：评估群体结构
    • 系统发育树：重建进化历史
    • ADMIXTURE：分析祖先组成

  2. 选择信号分析

    • 核苷酸多样性(π)：遗传多样性评估
    • FST：群体分化分析
    • XP-CLR/XPEHH：选择信号扫描
    • Tajima’s D：中性检验

  3. 全基因组关联分析(GWAS)

    • 线性混合模型
    • 多重检验校正
    • 候选基因鉴定
    • 功能富集分析

数据可视化

  • 曼哈顿图和QQ图（GWAS结果）
  • LD衰减曲线
  • PCA图和系统发育树
  • ADMIXTURE条形图
  • FST和π分布图
  • 选择信号扫描图

所有分析方法均提供了详细的命令行示例、参数说明和R语言可视化代码，可直接应用于实际数据分析。该流程适用于大豆及其他作物的重测序数据分析，为遗传学研究和分子育种提供了完整的技术方案。

生物信息学
#生物信息学
8000份大豆重测序数据文章的流程
https://lixiang117423.github.io/article/soybean8000/
作者
李详【Xiang LI】
发布于
2026年1月19日
许可协议