給定一串數據序列，如何判斷是否存在（模糊的）重複子序列？有什麼相應的演算法？
只要求判斷是否存在重複的子序列就可以了。下面這張圖存在重複子序列，我們肉眼很容易判斷，但是如何讓計算機自動識別？每個子序列之間的數據並不是完全相等，但是從整體上看，他們是重複出現的。之前看過從一長串字元串中查找重複子序列的問題，例如「abcdabcdabcd」存在重複子序列「abcd」，那個比較簡單。
不知道自相關函數和傅里葉變換對於這個問題是否有用？
這個問題簡單，用不著上機器學習演算法，無非多一個序列比對。題主要是能完全明白我在說什麼，自己很快就能寫出代碼來做這件事。
先重新梳理下題干：
在一個長度為L的序列中，識別是否存在一個子序列，在全長序列中出現了至少兩次，同時要求容許一定程度的錯配（Edit distance = m）。關鍵詞：Kmer步驟：
1. 設Kmer=k，從第一個鹼基起逐個取全長序列的長度為k的子序列，最終得到L-k+1條子序列，以序列本身為鍵，出現次數為值建立哈希表。2. 取出哈希表所有鍵，兩兩之間進行序列比對（序列比對演算法就不詳述了，教科書上有），凡是Edit distance小於等於m的，將鍵和值合併到值更高的那個去。3. 凡是值大於1的所對應的鍵，其序列便是你要找的重複序列。簡單吧？其實最難的是Kmer的取值。Kmer值決定了你找的重複序列的最短長度，也就是說，長度短於Kmer的重複序列將不能被檢測到。另一方面，對於由ATCG組成的全長序列，Kmer值又不能低於log4(L)，否則概率上講即使全長序列是隨機序列，你也有很大幾率找到所謂重複序列，而這樣的重複序列往往就失去生物學意義了。另一個難點是將「出現至少兩次」作為重複序列的標準有時候太寬了，題主可以根據具體問題適當調高閾值。生物信息學裡有這個功能的需求，因為要在基因組裡面找重複元件。EMBOSS工具包裡面有幾個重複序列的程序，對應不同形式的重複序列。你可以找找看。
類似於生物信息學裡面motif finding的問題 所謂的模體識別就是從若干條基因序列裡面找出確定長度的保守的片段 這些片段可能有一些不同。EM演算法和Gibbs Samper之類的方法可以解決這個問題。我自己寫的R實現的Gibbs Samper在這裡= =非常慢
3、
實驗作業：
Write a program in
R or Python to sample motif occurrence from sequences in Lab_3_crp.fa given the motif probability model:
____________________________________________
Pos.
# A T
C G
____________________________________________
1 |
0.295 0.506 0.072 0.127
2 |
0.295 0.506 0.019 0.180
3
| 0.348 0.348 0.124 0.180
4 |
0.032 0.822 0.124 0.022
5 |
0.032 0.243 0.124 0.601
6 |
0.032 0.874 0.072 0.022
7 |
0.190 0.085 0.019 0.706
8 |
0.769 0.085 0.124 0.022
Assume that there
is one binding site in each sequence. Repeat the sampling step 1000 times, and
then plot a histogram of where the bind sites are for the first two sequences
(「cole1」 and 「ecoarabop」).
Note 1. You also need to consider the complementary strands in your sampling
program.
Note 2. For this problem R code will be much easier to write, use the R
function called sample.
library(seqinr)
seq&<-read.fasta("/home/guoyifan/下載/Lab_3_crp.fa")
a&<-0
t&<-0
c&<-0
g&<-0
for(s in seq){
for(b in s){
if(b=="a"){a&<-a+1}
if(b=="t"){t&<-t+1}
if(b=="c"){c&<-c+1}
if(b=="g"){g&<-g+1}
}
}
bg&<-c(a/(a+t+c+g),t/(a+t+c+g),c/(a+t+c+g),g/(a+t+c+g))
#計算背景基因頻率
seq2&<-seq
for(s in (1:18)){
for(b in (1:105)){
if(seq2[[s]][b]=="a"){seq2[[s]][b]&<-"t"}
else if(seq2[[s]][b]=="t"){seq2[[s]][b]&<-"a"}
else if(seq2[[s]][b]=="c"){seq2[[s]][b]&<-"g"}
else if(seq2[[s]][b]=="g"){seq2[[s]][b]&<-"c"}
}
}
seq&<-c(seq,seq2)
#添加互補鏈
pre&<-matrix(c(0.295,0.506,0.072,0.127,0.295,0.506,0.019,0.180,0.348,0.348,0.124,0.180,0.032,0.822,0.124,0.022,0.032,0.243,0.124,0.601,0.032,0.874,0.072,0.022,0.190,0.085,0.019,0.706,0.769,0.085,0.124,0.022),nrow=4,byrow=F)
IC1=0
for(m in 1:4){
for(n in 1:8){
if(pre[m,n]!=0){
IC1&<-IC1+pre[m,n]*log((pre[m,n]/bg[m]),2)
}
}
}
#計算初始矩陣的IC
sc&<-1
mi&<-1
a&<-c(0,0,0,0,0,0,0,0)
t&<-c(0,0,0,0,0,0,0,0)
c&<-c(0,0,0,0,0,0,0,0)
g&<-c(0,0,0,0,0,0,0,0)
for(s in (1:36)){
for(b in (1:(105-7))){
sc[b]&<-1
for (i in (0:7)){
if(seq[[s]][b+i]=="a"){sc[b]&<-sc[b]*pre[1,(i+1)]/bg[1]}
if(seq[[s]][b+i]=="t"){sc[b]&<-sc[b]*pre[2,(i+1)]/bg[2]}
if(seq[[s]][b+i]=="c"){sc[b]&<-sc[b]*pre[3,(i+1)]/bg[3]}
if(seq[[s]][b+i]=="g"){sc[b]&<-sc[b]*pre[4,(i+1)]/bg[4]}
}
}
for(j in 1:length(sc)){
if(sc[j]==max(sc))
mi[s]&<-j
} }
#初始矩陣選出起始點組
for(i in 1:36){
a&<-c(0,0,0,0,0,0,0,0)
t&<-c(0,0,0,0,0,0,0,0)
c&<-c(0,0,0,0,0,0,0,0)
g&<-c(0,0,0,0,0,0,0,0)
psfm&<-matrix(0,nrow=4,ncol=8)
for(j in 1:36){
if(j!=i){
for(k in (mi[j]:(mi[j]+7))){
if(seq[[j]][k]=="a"){a[k-mi[j]+1]&<-a[k-mi[j]+1]+1}
if(seq[[j]][k]=="t"){t[k-mi[j]+1]&<-t[k-mi[j]+1]+1}
if(seq[[j]][k]=="c"){c[k-mi[j]+1]&<-c[k-mi[j]+1]+1}
if(seq[[j]][k]=="g"){g[k-mi[j]+1]&<-g[k-mi[j]+1]+1}
}
}
}
for(i2 in (1:8)){
psfm[1:4,i2]&<-c(a[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),t[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),c[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),g[i2]/(a[i2]+t[i2]+c[i2]+g[i2]))
}
#刪除序列 重新計算psfm
sc[i]&<-0
for(b in (1:(105-7))){
sc[b]&<-1
for(i1 in (0:7)){
if(seq[[i]][b+i1]=="a"){sc[b]&<-sc[b]*psfm[1,(i1+1)]/bg[1]}
if(seq[[i]][b+i1]=="t"){sc[b]&<-sc[b]*psfm[2,(i1+1)]/bg[2]}
if(seq[[i]][b+i1]=="c"){sc[b]&<-sc[b]*psfm[3,(i1+1)]/bg[3]}
if(seq[[i]][b+i1]=="g"){sc[b]&<-sc[b]*psfm[4,(i1+1)]/bg[4]}
}
}
#用新的psfm給刪除的序列打分
cnt&<-0
while(TRUE){
a&<-c(0,0,0,0,0,0,0,0)
t&<-c(0,0,0,0,0,0,0,0)
c&<-c(0,0,0,0,0,0,0,0)
g&<-c(0,0,0,0,0,0,0,0)
psfm1&<-matrix(0,nrow=4,ncol=8)
si&<-sample((1:(105-7)),1,prob=sc)
#按照得分採樣
mi[i]&<-si
psfm1&<-matrix(0,nrow=4,ncol=8)
for(j in 1:36){
for(k in (mi[j]:(mi[j]+7))){
if(seq[[j]][k]=="a"){a[k-mi[j]+1]&<-a[k-mi[j]+1]+1}
if(seq[[j]][k]=="t"){t[k-mi[j]+1]&<-t[k-mi[j]+1]+1}
if(seq[[j]][k]=="c"){c[k-mi[j]+1]&<-c[k-mi[j]+1]+1}
if(seq[[j]][k]=="g"){g[k-mi[j]+1]&<-g[k-mi[j]+1]+1}
}
}
for(i2 in (1:8)){
psfm1[1:4,i2]&<-c(a[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),t[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),c[i2]/(a[i2]+t[i2]+c[i2]+g[i2]),g[i2]/(a[i2]+t[i2]+c[i2]+g[i2]))
}
IC2=0
for(m in 1:4){
for(n in 1:8){
if(psfm1[m,n]!=0){
IC2&<-IC2+psfm1[m,n]*log((psfm1[m,n]/bg[m]),2)
}
}
}
#更新psfm和IC
if(IC2&>IC1){
psfm&<-psfm1
IC1&<-IC2
break
}
cnt&<-cnt+1
if(cnt&>=1000){break}
}
}
mi
[1] 45 72 76 83 85 60 24 63 9 14 61 41 66 1 17 53 71 40 51 68 60 76 6 35 55 72 45 27 19 45 30 84 72 57 80 89
題主可以試試樣本熵sample entropy，能表徵相同模式信息出現的多少。
信號處理方面有類似的演算法，所用的就是自相關函數。
方法是設定一個窗長範圍，然後不斷的移動窗口，變換窗口大小，在窗口內求前後半段的自相關，所求的的最大自相關值對應的窗長的一半就是那個子序列長度了。
根據設置的自相關門限可以實現精確匹配或者模糊匹配。
建議去看一下melp語音編碼演算法中關於整數基音周期求取相關部分。
推薦閱讀：
※現在在網上上傳一張圖片，如果到一定程度就可自動被系統識別為違禁圖片，這些工作肯定是機器做的，我想問的是：這種圖像識別系統的原理是什麼？
※中國的算命能不能算是模式識別？特徵是生辰八字，輸出是特定的命運？
※如何向外行人解釋模式識別和機器學習中的 Generalization 機制？
※有沒有值得推薦的隨機森林 Random Forest教材？
※EM演算法怎麼用在聚類上？