生信自學筆記（十） 多序列區域性聯配演算法

區域性聯配，所以對兩段相似度較高的鹼基序列就不要求處於臨近的位置，通過比較雜湊索引的偏移，我們可以方便地實現這一點。

對這種迭代演算法乙個通常的比喻是：要均勻地把飯分到兩個碗裡面，先隨機分配，看到那個碗裡面飯比較多，就從中舀一些到另外乙個碗裡面去，直到看上去兩個碗裡面的飯完全一樣為止。

選擇乙個聯配寬度

基於聯配結果和寬度構建初始pssm矩陣

利用該pssm對每條序列進行逐位點掃瞄

計算每條序列各個位點與pssm的匹配概率

合計一條序列的所有匹配概率值，將上一步算出的匹配概率與這個總和相除。

對所有序列進行上述運算

基於上述步驟獲得的匹配概率值更新該pssm

重複步驟4~8一百次以上，直到pssm各列的鹼基頻率不再變化為止。

示意圖見下：

import numpy as np

t = 'agct'
def get_pssm(s, ss):
# s為整體的dna序列的列表，而ss為要計算的區域性序列矩陣，返回乙個pssm矩陣，每一行對應的鹼基分別為atcg
width = ss.shape[1]
pssm = np.zeros((4, width))
# 統計每一列鹼基出現的次數,構建位置頻度矩陣 (pfm)
for col in range(width):
temp = ''.join(list(ss[:,col]))
for row in range(4):
pssm[row, col] = temp.count(t[row])
# 計算鹼基頻率,構建位置概率矩陣 (ppm)
for col in range(width):
s = sum(pssm[:,col])
pssm[:, col] = pssm[:, col]/s
k = ''.join(s)
background = 
for row in range(4):
# 返回pssm矩陣和背景鹼基頻率列表
return pssm, background
def get_pmatric(pssm, b, s, width, offset): #獲得最佳定位概率矩陣
p = np.full((len(s), len(s[0]) - width + 1), fill_value=1.0)
for n in range(1):
for i in range(len(s)):
for j in range(len(s[i]) - width + 1):
for k in range(len(s[i])):
if k in range(j+offset, j + width):
p[i][j] = p[i][j] * pssm[t.index(s[i][k])][k - j-offset]
else:
p[i][j] = p[i][j] * b[t.index(s[i][k])]
for i in range(4):
p[i, :] = p[i, :] / sum(p[i, :])
return p
def update(pssm, p, s, offset): #更新pssm矩陣
for i in range(p.shape[0]):
for j in range(p.shape[1]):
pssm[t.index(s[i][j+offset])][offset] = pssm[t.index(s[i][j+offset])][offset] + p[i][j]
s = ['aatcgattc','atacaatga', 'taaggaata', 'tataatcga']
# 輸入四條待聯配的序列，其比對的目標子串行為aat，即聯配寬度設定為3
width = 3
ss = 
# 隨機生成乙個多序列聯配結果
for i in s:
r = np.random.randint(len(i)-width)
ss = np.vstack(ss)
# 構建初始pssm矩陣
pssm, b = get_pssm(s, ss)
for n in range(100):
for offset in range(width):
p = get_pmatric(pssm, b, s, width, offset)
update(pssm, p, s, offset)
print(pssm)
p = get_pmatric(pssm, b, s, width, 0)
print('原始序列\t\t 對應聯配')
for i in range(p.shape[0]):
t = p[i,:]
start = np.where(t == max(t))[0][0]
print(s[i], '\t', s[i][start:start+width])

原始序列		 對應聯配

aatcgattc 	 aat
atacaatga 	 aat
taaggaata 	 aat
tataatcga 	 aat

生信自學筆記（一） 概述

想涉獵生物資訊學已經很久了，只是苦於平時作業繁重，沒有連續的時間坐下來好好學習。好不容易挨到暑假，是時候補一波知識了。自學教材 生物資訊學 樊龍江主編 當初在圖書館找教材，找了好幾種，出版社不同，年代不同，裝幀風格也非常不一樣。迷茫的我就去請教了研究生師姐，然鵝她說課本這種東西都是大同小異的，於是我...

Java自學筆記（十）

要用到多型，一定是已經有子父類關係或者類實現介面等前提 格式 父類型別 變數名稱 new 子型別行 變數名稱.方法 具體體現 父子類，抽象類，介面 class fu class zi extends fu 類的多型使用 fu f new zi 這其實就是向上轉型 abstract class fu ...

python自學筆記（十）

1.生成單調的list 1,2,3,4,5,6 可用list range 1,11 2.生成有一定規律的list 1x1,2x2,3x3,10x10 l for x in range 1,11 x x for x in range 1,11 把要生成的表示式放在前面 for迴圈 判斷條件 根據情況選...