文章目录
一、简介二、文献来源三、文献工作流程梳理以及本文复现目标四、文献复现4.1 数据下载4.2 kmer的生成与过滤4.3 多序列比对4.4 序列编码与建模4.5 模型验证 五、复现效果评估与优化方向一、简介
在转录组学数据挖掘中，转录因子结合位点预测算是比较重要的一个研究方向。在这个专栏中，我们将从的一篇文献开始，逐渐了解在转录因子结合位点预测方向的一些工作。
本文基于在Nucleic Acids Research上发表的文献《基于信念传播的DNA模体解析》，通过下载相应的数据集进行文献复现，体会隐马尔可夫模型（HMM）在生物信息学中的作用。
值得注意的是，由于文章没有提供源码，并且文献不会告诉你全部的细节。于是，本文仅仅下载了一个数据集并且对文章的部分流程进行复现。不过，通过这种方式依然能够对该文献的工作有一定的了解。
二、文献来源
Ka-Chun Wong and others, DNA motif elucidation using belief propagation,Nucleic Acids Research, Volume 41, Issue 16, 1 September , Page e153./10.1093/nar/gkt574
三、文献工作流程梳理以及本文复现目标
文献已经提供了Workflow，如下图所示：
首先，对数据源进行解读：
文献是基于PBM数据来分析的。PBM即Protein Binding Microarray，是一种类似DNA microarray的实验技术，与DNA microarray不同，PBM是研究蛋白和核酸探针的结合情况的，从而能够进一步探讨蛋白与序列结合的情况。如果使用的蛋白是转录因子，那么这个实验就是研究转录因子结合位点（一段序列）的，这一段序列就是“基序(Motif)”的一种。文献首先基于PBM数据，该数据的组成有两列：第一列是核酸探针（包含引物序列）的序列情况，第二列是该探针所表现的信号强度（归一化后）。下图展示了一个数据集的情况：
接下来，对Workflow进行一些解读：
文献第一步，以长度为 k k k的窗口依次读取探针序列，如果探针的总长度为 L L L，那么对于每一个探针来说，就会生成 L − k + 1 L-k+1 L−k+1个“窗口”，即kmer。位于同一条探针上的kmer对应的信号强度应当相同。
文献第二步，由于不同的探针有可能产生相同的kmer，于是将相同kmer的不同信号强度放在一起。
文献第三步，在上一步汇总完成后，对于每一个kmer的一组信号强度，取其中位数的信号强度作为该kmer的信号强度。（第三步图中未体现，正文中提及）
文献第四步，在进行比对之前，先对kmer进行过滤。设样本的信号强度序列 { I 1 , I 2 , . . , I m } \{I_1,I_2,..,I_m\} {I1,I2,..,Im}的中位数为 m e d med med，中位数绝对偏差(MAD)为 m a d mad mad，那么对于每个kmer的信号强度 I k m e r I_{kmer} Ikmer，我们只保留这样的kmer：
I k m e r > m e d + 4 ∗ m a d / 0.6745 I_{kmer}>med+4*mad/0.6745 Ikmer>med+4∗mad/0.6745
并且称这样的kmer为“正kmer”（Positive kmer）。（第四步图中未体现，正文中提及）
文献第五步，对过滤得到的“正kmer”进行多重比对。文献第六步，将比对后的结果作为HMM模型的输入，训练出模型的参数 θ \theta θ。文献使用了50个隐变量。文献第七步(a)，将训练好的模型进行验证。验证方法是，对新数据的每一条探针，生成 L − k + 1 L-k+1 L−k+1个长度为 k k k的kmer，使用模型预测生成每个kmer的概率，取最大的概率作为该条探针对转录因子的偏好性。最后，检验探针的偏好性与其信号强度的m a d / 0.6745 mad/0.6745 mad/0.6745其实就是对信号强度序列的标准差的估计值。
Spearman相关系数，就可以反映模型的有效性。文献第七步(b)，使用N-Max-Product算法找出最有可能的 N N N条隐状态路径，并输出对应的观测序列，作为对转录因子结合位点基序的预测。本文的复现目标：基于转录因子Zif268的PBM实验数据，使用Python实现上述第1~7(a)步的任务。
至于为什么不复现7(b)步，是因为发现复现效果与文献有较大出入，所以在此也不做展示。
四、文献复现
4.1 数据下载
本文主要基于转录因子Zif268的PBM实验数据，可在Uniprobe上下载。
Zif268数据介绍(点击访问)
Zif268数据下载界面(点击访问)
下载得到一个zip文件，解压后找到里面的Zif268_deBruijn_v1.txt和Zif268_deBruijn_v2.txt，这两个文件就是我们需要处理的文件了，分别将其命名为Zif268_1.txt和Zif268_2.txt。我们将会以1号文件为训练集进行模型参数训练，并使用2号文件作为模型的验证数据集。
首先，我们需要导入相关的包：
import pandas as pdimport numpy as npfrom hmmlearn import hmmimport matplotlib.pyplot as plt
接着，使用pandas包读入文件：
def read_file(path: str):"""读取文件生成Dataframe并预处理path: PBM文件路径"""f = pd.read_csv(path, sep="\t")f.columns = ['Seq', 'Intensity']f['Seq'] = f['Seq'].str.strip()f['Intensity'] = pd.to_numeric(f['Intensity'], errors='coerce')return ff1 = read_file("Zif268_1.txt")f2 = read_file("Zif268_2.txt")
这就完成了数据准备工作。
4.2 kmer的生成与过滤
接下来，我们需要对数据f1进行kmer的提取和过滤。
首先，我们定义探针长度为35（建模时不需要引物序列），kmer长度为8：
kmer_len = 8probe_len = 35
之后依次读取kmer，并收集其对应的信号强度。收集完毕后，取中值作为该kmer的信号强度：
def probe2kmer(data: pd.DataFrame):"""将探针数据转换成kmerdata: 原始数据, 使用pandas数据框读取k: kmer的长度probe_len: 探针的长度"""# 生成kmer字典kmer_dict = {}for row in range(len(data)):seq = data['Seq'][row]val = data['Intensity'][row]if val != float("nan"):j = kmer_lenwhile j <= probe_len:kmer = seq[j-kmer_len:j]if kmer not in kmer_dict.keys():kmer_dict[kmer] = []kmer_dict[kmer].append(float(val))j += 1# 取信号中值for key, value in kmer_dict.items():kmer_dict[key] = np.median(value)return kmer_dictkmer_dict = probe2kmer(f1)
值得注意的是，原数据里面有部分探针的信号强度为空，在上一步载入数据时，已经将其转化为nan了。遇到这种探针，本文的处理是将其抛弃，只保留有数值的kmer。
上述程序运行完毕后，我们就得到了kmer字典，其key为kmer序列，value为对应的信号强度。这一步获得的kmer数量是65536个。
接下来，我们需要过滤获得“正kmer”。在此之前，我们设计一个函数，用于获取数据框中某一列的中位数和MAD值：
def get_param(data: pd.DataFrame, name: str):"""返回数据框某一列的中位数及其中位数绝对偏差data: 原始数据的数据框name: 列名"""median_val = data[name].median()mad_val = (data[name] - data[name].median()).abs().median()return median_val, mad_val
这样就可以根据筛选规则来筛选“正kmer”了：
def pos_kmer(kmer_dict: dict, med: float, mad: float):"""kmer过滤为正kmerkmer_dict: kmer字典med: 所有信号的中位数mad: 所有信号的中位数绝对偏差"""kmer_pos = []for key, value in kmer_dict.items():if value > med+4*mad/0.6745:kmer_pos.append(key)return kmer_posmed, mad = get_param(f1, "Intensity")kmer_pos = pos_kmer(kmer_dict, med, mad)
最终我们获得了一个“正kmer”的列表，共计444个。
4.3 多序列比对
获得“正kmer”后，就需要进行多序列比对了。这里笔者并不打算使用Python接口去调用比对程序，而是将kmer列表输出后，手动去网站比对：
def output(data: list):"""将过滤的kmer列表输出为fasta"""f = open("./kmer.fasta", "w")n = 1for i in data:f.write(">Seq%d\n%s\n"%(n, i))n += 1f.close()output(kmer_pos)
获得输出的fasta文件后，笔者使用ClustalOmega网站进行了多序列比对。
Clustal Omega (点击访问)
在比对时选择DNA序列，传入生成的fasta文件，其他参数保持默认，提交。获得比对文件后下载，命名为alignment.txt：
对于文件后面的“8”，笔者认为可能是我提交的kmer的长度；对于上述比对过程，可以使用biopython库里的接口完成，这里不阐述。4.4 序列编码与建模
获得多重比对后，我们读取保留中间的序列部分，准备开始建立模型了：
align = open("alignment.txt").readlines()align = [i.split("\t")[0].split(" ")[-1] for i in align]
读取时注意去掉标题和空行，只保留比对数据本身。
由于本次模型是离散型的隐马尔可夫模型，所以首先要对序列中出现的字符进行编码：
def labelencode(kmer_list: list):"""将kmer序列编码kmer_list: 由kmer组成的列表"""label = {'-': 0, 'A': 1, 'G': 2, 'C': 3, 'T': 4}obseq = []for i in kmer_list:seq = [label[j] for j in list(i)]obseq.append(seq)return np.array(obseq)x = labelencode(align)
对于每一个kmer，我们生成对应的编码，并且需要转化为np.array的数据类型。
接下来就可以开始建模了。对于离散型的HMM模型，我们使用hmmlearn中的CategoricalHMM()函数进行建模，设定隐状态数量为50个（与文献一致），迭代次数为20次。
# 训练model = hmm.CategoricalHMM(n_components=50, n_iter=20)model.fit(x)
由于隐状态数量较多，所以训练耗时在13秒左右。
4.5 模型验证
已经获得模型了，接下来就使用f2数据集进行验证。
首先，我们需要计算每一条序列的偏好性：
def prefer_test(data: pd.DataFrame, model: hmm.CategoricalHMM):"""准备偏好测试data: 用于测试的原始数据model: 拟合的模型"""prob_list = []val_list = []n = 0for i in list(data['Seq']):kmer_list = [i[j:j+kmer_len] for j in range(probe_len-kmer_len+1)]temp_prob = [np.exp(model.score(labelencode([i]))) for i in kmer_list]temp_prob = [i for i in temp_prob if not np.isnan(i)]temp_prob.append(0)prob_list.append(max(temp_prob))if not np.isnan(data['Intensity'][n]):val_list.append(data['Intensity'][n])else:val_list.append(0)n += 1return prob_list, val_listprob_list, val_list = prefer_test(f2, model)
代码中kmer_list用于存储每一条探针的所有kmer，并使用model.score()方法对每一条kmer打分。由于可能存在nan，所以我们只保留非nan的部分。同样地，对于不存在信号强度值的点，我们插入值0。
为了保证一条探针的偏好性存在，笔者特意为
temp_prob变量添加了0元素，防止temp_prob被过滤后出现空集的现象，否则无法计算max。
值得注意的是，由于数据量较大，上述代码会运行6分钟左右不等，这也是值得优化的地方。
获得每一条探针的偏好性列表prob_list和对应的信号强度列表val_list后，我们就可以计算其Spearman相关系数并画图了。相关系数的计算直接采用数据框的corr方法，可视化以散点图为载体：
def corr_plot(x: list, y: list):"""计算相关性并画图"""# 相关性result = pd.DataFrame({'Probability':x,'Intensity':y})result = result.corr(method='spearman')coef = result['Probability'][1]# 散点图fig, ax = plt.subplots()ax.spines['right'].set_visible(False)ax.spines['top'].set_visible(False)plt.scatter(x, y, c=y, cmap='plasma', s=10)plt.xlabel("Binding Preference")plt.ylabel("Probe Intensity")plt.title("Scatter Plot")plt.text(2.5e-9, 60000, "Spearman Correlation = %.3f"%coef)plt.show()corr_plot(prob_list, val_list)
可以看到，偏好性值与信号强度的Spearman相关性值为0.255，而文献中给出的相关性值为0.336，复现效果还是将就的。
五、复现效果评估与优化方向
复现效果评估：本文对验证集的复现效果还是不错的，在没有得到文章实现细节的情况下，仍然能够达到一个可观的效果；同时，这也反映了隐马尔可夫模型在处理多模态序列识别方面是具有优势的。
优化方向：
HMM模型的训练次数可以增加，并且可以增加交叉验证的内容；缺乏了N-Max-Product算法的实现；划分单个kmer的方法会忽略序列上下文信息，这是该模型的限制。后续如果能够提高k值，并且借助ChIP-seq等NGS数据进行分析，也许会有更进一步地效果。
