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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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.names

cnt$N= rowSums(cnt[, c("ref.matches", "alt.matches")])

# filter counts
cnt = subset(cnt, N >= 10)

# add end position and name
cnt$end= cnt$start
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

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 site
  p1= dbetabinom(x[1], n[1], m= 0.5 + allelic.imbalance, s= d)
  
  # for subsequent sites
  len = length(x)
  
  if(len > 1) {
    
    # precompute likelihoods of each subsequent SNP given 'in-phase' with first SNP
    snp.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 SNP
    snp.phase0.like <- dbetabinom(x[2:len], n[2:len], m=0.5 - allelic.imbalance, s = d)
    
    # create phase array
    phase1.like.array <- rep(NA, len)
    phase0.like.array <- rep(NA, len)
    
    # add likelihood for first site
    phase1.like.array[1] <- p1
    phase0.like.array[1] <- p1
    
    for(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 phases
      
      prev <- (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 elements
    l = -log(0.5*phase1.like.array[len] + 0.5*phase0.like.array[len])
    
  } else {
    
    l = -sum(log(p1))
    
  }
  
  # return
  return(l)
}

# function to estimate allele imbalance

fun = function(i){
  
  # test region
  gr2= reg[i]
  
  # subset by overlaps
  dat= subsetByOverlaps(gr1, gr2, type = "any") %>% as.data.frame()
  
  if(nrow(dat)> 3){
    
    # sort by start
    dat=dat[order(dat$start),]
    
    rownames(dat)<- NULL
    
    m1<- optimize(ll.new, c(-0.49, 0.49), x= dat$ref.matches, n= dat$N, d= d)
    
    # result
    res= data.frame(seqnames= seqnames(gr2),start= start(gr2), end= end(gr2), name= gr2$name, a =  m1$minimum ,nsites= nrow(dat))
    
  } else{
    
    res<- NULL
    
  }
  res
}
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

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)
1
2
3
4
5
6
7

这part的demo特针对于DNA数据，代码已debug，根据各自搭建的环境运行，记得提前加载R包。

下期出RNA的处理～