R语言利用vcf文件计算等位基因频率和连锁不平衡（LD）R方

代码来源于论文
Assessment of linkage disequilibrium patterns between structural variants and single nucleotide polymorphisms in three commercial chicken populations
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-022-08418-7
vcf示例文件用之前介绍rMVP包做GWAS分析的那期推文的数据
首先使用beagle做基因型填充
beagle gt=smoove_filtered.vcf out=smoove.filtered.impute nthreads=2
读取vcf文件

library(data.table)
library(tidyverse)
dat.map<-fread(“smoove.filtered.impute.vcf.gz”,skip = “#CHR”)
把 0|1 这种基因型拆分成两列

gt<-data.table()

for(i in 10:ncol(dat.map)){
tmp<-str_split(dat.map[[i]],fixed(“|”),simplify = TRUE) %>%
as.data.frame()
colnames(tmp)<-paste0(colnames(dat.map)[i],c(“_1″,”_2”))

gt[,colnames(tmp):=tmp]
}
这里有一个:=这个符号，暂时没有搞明白这个写法是什么意思，可以一直把列添加到一个数据框里
以下代码把数据框转化成了一个列表

gt %>%
t() %>%
as.data.table() %>%
unclass() -> gt.list

class(gt.list)
计算等位基因频率

p<-list() n<-gt.list[[1]] %>% length()

for(i in 1:length(gt.list)){
p[[i]] <- table(gt.list[[1]])/n
}
自定义计算LD的函数

library(compiler)
calcLD <- cmpfun(function(x,pa,ht,p){ n<-length(x) ht_int <- lapply(ht,as.integer) R2 <- list() if(is.list(p)){ biv <- which(unlist(lapply(ht,function(x){length(levels(x))}))==2) if(length(biv) > 0){
pb <- unlist(lapply(p[biv],function(x){return(x[1])}))
pab <- unlist(lapply(ht_int[biv],function(y,x,n){sum(x==1 & y==1)/n},x=as.integer(x),n=n))
D <- pab – pa[1] * pb
R2 <- D^2/ ((pa[1] * pb) * ((1-pa[1]) * (1-pb)))
}
}else{
y <- ht
pb <- p

pab <- matrix(0,nlevels(x),nlevels(y))
rownames(pab) <- levels(x)
colnames(pab) <- levels(y)

if(sum(dim(pab))==4){
  pab <- sum(x == rownames(pab)[1] & y == colnames(pab)[1]) / n
  D <- pab - pa[1] * pb[1]
  R2[[1]] <- D^2/ prod(pa,pb) 

}else{
  for(i in 1:nrow(pab)) for(j in 1:ncol(pab)){
    pab[i,j] <- sum(x == rownames(pab)[i] & y == colnames(pab)[j]) / n 
  }
  D <- pab - pa %*% t(pb)
  R2[[1]] <- D^2/ ((pa %*% t(pb)) * ((1-pa) %*% (1-t(pb)))) # not correct yet? How to define multiallelic LD?
}

}
return(R2)
})
整个函数的逻辑还看不明白
这里自定义函数还用到了compiler这个R包，有什么作用暂时不太明白
函数是输入两个位点的等位基因和等位基因频率
calcLD(gt.list[[1]],p[[1]],gt.list[[3]],p[[3]])
gt.list 的格式

image.png
p的数据格式

image.png
以上是本期推文的内容
一个R语言的零散知识点：pivot_longer()函数把多列的数据转换成长格式
library(tidyverse)
df <- data.frame(
id = 1:3,
var1 = c(“A”, “B”, “C”),
var2 = c(“D”, “E”, “F”),
value1 = c(10, 20, 30),
value2 = c(100, 200, 300)
)

df %>%
pivot_longer(cols = c(var1,var2),
names_to = “ABCDE”) %>%
pivot_longer(cols = c(value1,value2),
values_to = “p”)
cols 参数的作用是 把向量里的两个列名单独生成一列
cols 里的列如果数据类型不一样是不能合并的
names_to 生成的是新生成的列的列名 values_to 也是指定列名
欢迎大家关注我的公众号
小明的数据分析笔记本