R语言处理DNA等位基因不平衡（一）

在生物信息学和基因组学研究中，等位基因不平衡分析是一种重要的方法，用于识别在特定生物过程或疾病状态中可能受到选择压力的基因或基因区域。等位基因不平衡（Allele Imbalance）指的是基因座上两个等位基因表达或存在的比例不等，这种不平衡可能是由于自然选择、遗传漂变或基因流等进化力量的作用。

1. 背景：

基因的等位基因表达不平衡可以揭示关于细胞类型特异性、基因调控机制、遗传变异对表型影响的重要信息。

例如，在肿瘤细胞中，特定基因的等位基因不平衡可能指示肿瘤发展过程中的选择性优势或遗传不稳定性。通过对来自不同组织、不同疾病状态或不同发展阶段的样本进行等位基因不平衡分析，研究人员可以识别出关键的遗传变异，这些变异可能在生物过程的调控或疾病的发生发展中起着关键作用。

本研究旨在通过对特定基因区域内的等位基因不平衡进行系统分析，来识别可能在特定生物过程或疾病中发挥作用的基因或基因区域。利用高通量测序技术产生的大量基因组数据，我们可以细致地分析个体或细胞群中的基因变异，并评估这些变异之间的不平衡关系。

2. demo

2.1 导入包和数据

读取计数文件（例如，某种类型的基因表达或变异数据），并对数据进行初步处理，比如过滤。

args = commandArgs(trailingOnly = TRUE)
filename.count= args[1]
filename.region= args[2]
filename.output= args[3]
print(args)# load libraries
suppressMessages(library(data.table))
suppressMessages(library(plyr))
suppressMessages(library(dplyr))
suppressMessages(library(rmutil))
suppressMessages(library(rtracklayer))
suppressMessages(library(pbapply))
suppressMessages(library(parallel))# edit counts
cnt= read.table(filename.count, header=T, skip = 87)

2.2 数据预处理

通过bedtools merge命令（在R中使用系统调用来执行），合并相邻的或重叠的区域，以减少数据中的冗余。

# select colnames
col.names= c("seqnames", "start", "ref.matches", "alt.matches", "ref", "alt")# edit counts
cnt= read.table(filename.count, header=T, skip = 87)colnames(cnt) = col.namescnt$N= rowSums(cnt[, c("ref.matches", "alt.matches")])# filter counts
cnt = subset(cnt, N >= 10)# add end position and name
cnt$end= cnt$start

2.3 等位基因不平衡分析

对每个区域内的变异进行等位基因不平衡测试，使用贝塔-二项式模型来估计每个位点的等位基因频率。
使用优化方法（如optimise和optimize）来找到使似然函数最大化的等位基因不平衡参数。

# function to estimate dispersion under assumption of no allele-imbalance
mtmp <- function(par, x, n){# likelihood given het sites p= dbetabinom(x, n, m= 0.5, s= par)# maximize likelihood for being het-sum(log(p))
}# estimate dispersion
print("estimating dispersion")m0<- optimise(mtmp, c(1e-05, 100), x= cnt$ref.matches, n= cnt$N)d<- m0$minimum# convert to granges
gr1= cnt %>% makeGRangesFromDataFrame(keep.extra.columns = T)# test regions
reg= import.bed(filename.region)reg$name= paste0(seqnames(reg), "_" , start(reg), "_" , end(reg))# function to estimate allele proportions
ll.new<- function(par, x, n, d){allelic.imbalance <- par# for first sitep1= dbetabinom(x[1], n[1], m= 0.5 + allelic.imbalance, s= d)# for subsequent siteslen = length(x)if(len > 1) {# precompute likelihoods of each subsequent SNP given 'in-phase' with first SNPsnp.phase1.like <- dbetabinom(x[2:len], n[2:len], m=0.5 + allelic.imbalance, s = d)# precompute likelihoods of each subsequent SNP given 'out-of-phase' with first SNPsnp.phase0.like <- dbetabinom(x[2:len], n[2:len], m=0.5 - allelic.imbalance, s = d)# create phase arrayphase1.like.array <- rep(NA, len)phase0.like.array <- rep(NA, len)# add likelihood for first sitephase1.like.array[1] <- p1phase0.like.array[1] <- p1for(i in 2:len) {# prior SNP was either in-phase or out-of-phase with first SNP, consider# both mutually exclusive possibilities when computing combined likelihood of# all possible combinations of phasesprev <- (0.5 * phase1.like.array[i-1]) + (0.5 * phase0.like.array[i-1])phase1.like.array[i] <- prev * snp.phase1.like[i-1]phase0.like.array[i] <- prev * snp.phase0.like[i-1]}# total likelihood is sum of last two elementsl = -log(0.5*phase1.like.array[len] + 0.5*phase0.like.array[len])} else {l = -sum(log(p1))}# returnreturn(l)
}# function to estimate allele imbalancefun = function(i){# test regiongr2= reg[i]# subset by overlapsdat= subsetByOverlaps(gr1, gr2, type = "any") %>% as.data.frame()if(nrow(dat)> 3){# sort by startdat=dat[order(dat$start),]rownames(dat)<- NULLm1<- optimize(ll.new, c(-0.49, 0.49), x= dat$ref.matches, n= dat$N, d= d)# resultres= data.frame(seqnames= seqnames(gr2),start= start(gr2), end= end(gr2), name= gr2$name, a =  m1$minimum ,nsites= nrow(dat))} else{res<- NULL}res
}

2.4 结果合成与输出

将每个区域的分析结果合并成一个列表，然后将该列表转换成数据框，最后写入CSV文件以供后续分析使用。

my.list=llply(1:length(reg), fun, .progress = progress_text(char="+"))# combine list
res= rbindlist(my.list) %>% as.data.frame()# output results
write.csv(res , filename.output, row.names = F)