怎樣使用R語言利用vcf格式文件計算核苷酸多樣性

# 怎樣使用R語言利用vcf格式文件計算核苷酸多樣性

## 摘要  
核苷酸多樣性（π）是群體遺傳學中衡量遺傳變異的重要指標，通過分析VCF格式的基因組數據可以高效計算該參數。本文詳細介紹使用R語言從VCF文件處理到π值計算的完整流程，包括數據導入、質量過濾、等位基因頻率計算及統計檢驗方法，并提供可復用的代碼示例。

---

## 1. 背景知識

### 1.1 核苷酸多樣性的定義  
核苷酸多樣性（π）指在群體中隨機選取的兩個個體在特定基因組位點上出現不同核苷酸的概率，計算公式為：

\[ \pi = \frac{\sum_{i<j} 2p_i p_j}{n(n-1)/2} \]

其中：
- \( p_i, p_j \) 為等位基因頻率
- \( n \) 為樣本數

### 1.2 VCF文件格式  
Variant Call Format (VCF) 是存儲基因組變異的標準格式，包含：
- **元信息行**（##開頭）
- **標題行**（#CHROM POS...）
- **數據行**（每行代表一個變異位點）

關鍵字段：
- `CHROM/POS`: 染色體位置
- `REF/ALT`: 參考/替代堿基
- `FORMAT`: 樣本基因型編碼方式

---

## 2. 準備工作

### 2.1 軟件環境配置
```r
install.packages(c("vcfR", "adegenet", "pegas", "dplyr"))
library(vcfR)
library(adegenet)
library(pegas)
library(dplyr)

2.2 示例數據獲取

從1000 Genomes Project下載測試VCF文件：

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

3. 數據處理流程

3.1 VCF文件導入與解析

vcf <- read.vcfR("ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz")

# 查看基本信息
head(vcf@fix)  # 查看變異位點信息
vcf@gt[1:5, 1:5]  # 查看前5個樣本的基因型

3.2 數據質量控制

過濾低質量位點

# 根據GATK推薦標準過濾
vcf_filtered <- vcf[which(vcf@fix[,6] > 30 &  # QUAL > 30
                         grepl("PASS", vcf@fix[,7])), ] # FILTER=PASS

移除缺失數據

geno <- extract.gt(vcf_filtered)
geno[geno == "./."] <- NA
missing_rate <- apply(geno, 1, function(x) sum(is.na(x))/ncol(geno)
vcf_clean <- vcf_filtered[missing_rate < 0.2, ]  # 保留缺失率<20%的位點

4. 核苷酸多樣性計算

4.1 等位基因頻率計算

# 轉換為genind對象
genind_obj <- vcfR2genind(vcf_clean)

# 計算等位基因頻率
allele_freq <- adegenet::makefreq(genind_obj)

4.2 π值計算實現

單個位點計算

calc_pi_site <- function(freq){
  sum(freq * (1-freq)) * 2  # 二倍體校正
}

pi_per_site <- apply(allele_freq, 1, calc_pi_site)

滑動窗口統計（100kb窗口）

# 獲取位點坐標
positions <- as.numeric(vcf_clean@fix[,2])

# 創建窗口劃分
windows <- seq(1, max(positions), by=1e5)

# 計算窗口內π值
window_pi <- sapply(windows, function(w){
  idx <- which(positions >= w & positions < w+1e5)
  if(length(idx)>10){  # 至少包含10個有效位點
    mean(pi_per_site[idx], na.rm=T)
  } else { NA }
})

4.3 可視化分析

plot(windows/1e6, window_pi, type='l', 
     xlab="Position (Mb)", ylab="Nucleotide diversity (π)",
     main="Chromosome 22 Diversity Landscape")
abline(h=mean(window_pi, na.rm=T), col="red", lty=2)

5. 高級分析方法

5.1 群體間分化比較

# 按群體分組計算
pop_map <- read.table("population_labels.txt") # 需準備群體標簽文件
genind_obj$pop <- pop_map$V2

# 使用pegas計算群體間π
pi_by_pop <- pegas::nuc.div(genind_obj, pairwise=FALSE)

5.2 統計檢驗

# 檢驗群體間差異顯著性
pairwise.wilcox.test(pi_by_pop, pop_map$V2, p.adjust.method="fdr")

6. 常見問題解決

6.1 內存不足處理

對于大型VCF文件：

# 分塊讀取處理
vcf_chunk <- read.vcfR("large_file.vcf.gz", nrows=10000)

6.2 缺失數據處理建議

• 采用EM算法估計缺失基因型
• 使用SNPRelate等專業包處理

7. 完整代碼示例

# 完整流程腳本
library(vcfR)
vcf <- read.vcfR("input.vcf.gz")
geno <- extract.gt(vcf)
vcf_clean <- vcf[rowSums(is.na(geno)) < ncol(geno)*0.2, ]

library(adegenet)
genind_obj <- vcfR2genind(vcf_clean)
pi_values <- nuc.div(genind_obj)

write.csv(pi_values, "pi_results.csv")

參考文獻

1. Knaus BJ, Grünwald NJ (2017) vcfR: an R package to manipulate and visualize VCF format data. Mol Ecol Resour
2. Paradis E (2010) pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics

”`

該文檔包含： 1. 理論背景說明 2. 分步驟的代碼實現 3. 可視化示例 4. 常見問題解決方案 5. 完整可執行的代碼模板

實際應用時需根據具體數據調整參數，建議先使用小規模測試數據驗證流程。