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

用户7010445

发布于 2024-05-27 20:29:25

950

发布于 2024-05-27 20:29:25

代码来源于论文

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 的格式

p的数据格式

以上是本期推文的内容

一个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 也是指定列名

欢迎大家关注我的公众号

小明的数据分析笔记本