genome辅导、讲解R语言、R编程语言辅导、讲解MAP segmentation 辅导Web开发|解析R语言编程
- 首页 >> 其他 Introduction
We are interested in the segmentation of a sequence of DNA over a genome. A sequence of DNA is
composed of four bases “A”, “C”, “G” and “T”. The objective is to automatically segment the DNA
sequence into regions that are AT-rich, and regions that are CG-rich. The data come from the R
package seqinr.
#install.packages("seqinr") # if needed
library("seqinr")
ff <- system.file("sequences/ECOUNC.fsa", package = "seqinr")
seq <- read.fasta(ff)
dnaseq = seq$ECOUNC
dnaseq[1:50] # DNA sequence
## [1] "a" "a" "a" "g" "c" "a" "a" "a" "t" "a" "a" "a" "a" "t" "t" "t" "a"
## [18] "a" "t" "t" "t" "t" "t" "a" "t" "c" "a" "a" "a" "a" "a" "a" "a" "t"
## [35] "c" "a" "t" "a" "a" "a" "a" "a" "a" "t" "t" "g" "a" "c" "c" "g"
# Convert letters to numbers
y =c()
y[which(dnaseq=="a")] = 1
y[which(dnaseq=="c")] = 2
y[which(dnaseq=="g")] = 3
y[which(dnaseq=="t")] = 4
y[1:50]
## [1] 1 1 1 3 2 1 1 1 4 1 1 1 1 4 4 4 1 1 4 4 4 4 4 1 4 2 1 1 1 1 1 1 1 4 2
## [36] 1 4 1 1 1 1 1 1 4 4 3 1 2 2 3
n = length(y)
plot(y[1:1500],xlab="Genome Index", ylab="Base", , yaxt="n")
ytick<-seq(1, 4, by=1)
axis(side=2, at=ytick, labels = FALSE)
text(par("usr")[1], ytick,
labels = c('A','C', 'G', 'T'), srt = 45, pos = 2, xpd = TRUE)
0 500 1000 1500
Genome Index
We consider a hidden Markov model (X1, Y1,...,Xn, Yn) where (Y1,...,Yn) is the DNA sequence,
with Yt 2 {1, 2, 3, 4} (1 = A, 2 = C, 3 = G, 4 = T) and Xt 2 {1, 2} are the hidden states, taking value1 (AT-rich) or 2 (CG-rich). We define the state transition probabilities as
(1)
where p = 0.995 and uniform initial probability mass function.
The emission probabilities are given by Pr(Yt = j | Xt = i) = Bi,j where
B = 0.3 0.2 0.2 0.3
0.25 0.25 0.3 0.2. (2)
1 Exercise 1
(a) Use the viterbi R function provided in the lecture notes to calculate the maximum a posteriori
(MAP) estimate of the hidden state sequence
xb1:n = arg max x1:n2{1,2}n Pr(X1:n = x1:n|y1:n).
Check the numerical value of the message mt within that function. Discuss any numerical issue
arising. Propose and implement a modification of that R function to resolve this numerical issue,
and report the obtained MAP segmentation.
(b) Run the R function alpha recursion given in the lecture notes on the dataset and plot the filtering
probability mass function Pr(Xt = k|y1:t) over the genome index t for each class k = 1, 2. Again,
what numerical issue arise with this implementation?
(c) Create a novel function alpha2 recursion which takes the same inputs (y,mu,A,B) and returns the
filtering pmfs at time t = 1,...,n, using the following two-steps recursion, for t= 2,...,nt(xt) := p(xt, yt|y1:t1)= p(yt|xt)X
xt1
p(xt|xt1)p(xt1|y1:t1)
p(xt|y1:t) = t(xt)
with 1(x1) = p(y1|x1)p(x1) and p(x1|y1) = P1(x1)
x 1(x) .
(d) Run the function alpha2 recursion on the data, and plot the resulting filtering pmf over the
genome.
2 Exercise 2
(a) Obtain the MAP segmentation for di?erent values of p. Discuss how the segmentation relates
to the value of this parameter.
(b) Noting that
log(p(y1:n)) = log(p(y1)) +Xn
modify the function alpha2 recursion so that it also returns the loglikelihood.
(c) Using a grid-search approach, find the approximate value of p that maximizes the loglikelihood.
Briefly discuss any alternative approach that you could use to find the maximum likelihood
estimate of p (you do not need to implement this alternative approach).
We are interested in the segmentation of a sequence of DNA over a genome. A sequence of DNA is
composed of four bases “A”, “C”, “G” and “T”. The objective is to automatically segment the DNA
sequence into regions that are AT-rich, and regions that are CG-rich. The data come from the R
package seqinr.
#install.packages("seqinr") # if needed
library("seqinr")
ff <- system.file("sequences/ECOUNC.fsa", package = "seqinr")
seq <- read.fasta(ff)
dnaseq = seq$ECOUNC
dnaseq[1:50] # DNA sequence
## [1] "a" "a" "a" "g" "c" "a" "a" "a" "t" "a" "a" "a" "a" "t" "t" "t" "a"
## [18] "a" "t" "t" "t" "t" "t" "a" "t" "c" "a" "a" "a" "a" "a" "a" "a" "t"
## [35] "c" "a" "t" "a" "a" "a" "a" "a" "a" "t" "t" "g" "a" "c" "c" "g"
# Convert letters to numbers
y =c()
y[which(dnaseq=="a")] = 1
y[which(dnaseq=="c")] = 2
y[which(dnaseq=="g")] = 3
y[which(dnaseq=="t")] = 4
y[1:50]
## [1] 1 1 1 3 2 1 1 1 4 1 1 1 1 4 4 4 1 1 4 4 4 4 4 1 4 2 1 1 1 1 1 1 1 4 2
## [36] 1 4 1 1 1 1 1 1 4 4 3 1 2 2 3
n = length(y)
plot(y[1:1500],xlab="Genome Index", ylab="Base", , yaxt="n")
ytick<-seq(1, 4, by=1)
axis(side=2, at=ytick, labels = FALSE)
text(par("usr")[1], ytick,
labels = c('A','C', 'G', 'T'), srt = 45, pos = 2, xpd = TRUE)
0 500 1000 1500
Genome Index
We consider a hidden Markov model (X1, Y1,...,Xn, Yn) where (Y1,...,Yn) is the DNA sequence,
with Yt 2 {1, 2, 3, 4} (1 = A, 2 = C, 3 = G, 4 = T) and Xt 2 {1, 2} are the hidden states, taking value1 (AT-rich) or 2 (CG-rich). We define the state transition probabilities as
(1)
where p = 0.995 and uniform initial probability mass function.
The emission probabilities are given by Pr(Yt = j | Xt = i) = Bi,j where
B = 0.3 0.2 0.2 0.3
0.25 0.25 0.3 0.2. (2)
1 Exercise 1
(a) Use the viterbi R function provided in the lecture notes to calculate the maximum a posteriori
(MAP) estimate of the hidden state sequence
xb1:n = arg max x1:n2{1,2}n Pr(X1:n = x1:n|y1:n).
Check the numerical value of the message mt within that function. Discuss any numerical issue
arising. Propose and implement a modification of that R function to resolve this numerical issue,
and report the obtained MAP segmentation.
(b) Run the R function alpha recursion given in the lecture notes on the dataset and plot the filtering
probability mass function Pr(Xt = k|y1:t) over the genome index t for each class k = 1, 2. Again,
what numerical issue arise with this implementation?
(c) Create a novel function alpha2 recursion which takes the same inputs (y,mu,A,B) and returns the
filtering pmfs at time t = 1,...,n, using the following two-steps recursion, for t= 2,...,nt(xt) := p(xt, yt|y1:t1)= p(yt|xt)X
xt1
p(xt|xt1)p(xt1|y1:t1)
p(xt|y1:t) = t(xt)
with 1(x1) = p(y1|x1)p(x1) and p(x1|y1) = P1(x1)
x 1(x) .
(d) Run the function alpha2 recursion on the data, and plot the resulting filtering pmf over the
genome.
2 Exercise 2
(a) Obtain the MAP segmentation for di?erent values of p. Discuss how the segmentation relates
to the value of this parameter.
(b) Noting that
log(p(y1:n)) = log(p(y1)) +Xn
modify the function alpha2 recursion so that it also returns the loglikelihood.
(c) Using a grid-search approach, find the approximate value of p that maximizes the loglikelihood.
Briefly discuss any alternative approach that you could use to find the maximum likelihood
estimate of p (you do not need to implement this alternative approach).