Python应用——分析序列k-mer

什么是k-mer

在生物学里面，k-mer就是指一段生物学序列的长度为k的子串。例如：一段 ACGAGGTACGA 的DNA序列，其4-mer就是如下图所示：

4-mers
k-mer在生物学分析中有非常重要的应用，例如序列组装、序列相似性分析等。具体内容可以参见 wikipedia k-mer中关于k-mer的介绍。

以DNA序列为例，一般来说，一段由 L 个碱基组成的序列，会有 L-k+1 个k-mer，其可能的k-mer一共有 4^k 种。

分析步骤

在这里我们的分析步骤包括：

• 读取序列fasta文件
• 计算序列k-mer
• 对序列k-mer进行统计

读取序列fasta文件

首先需要了解fasta文件基本格式：fasta格式是一种非常简单的储存序列的格式，可以储存核酸序列（DNA/RNA）也可以储存蛋白质的氨基酸序列，主要分成2个部分。第一部分是以“>”为开始的一行主要储存的是序列的描述信息；剩下的是序列部分。例如下面就是一段氨基酸序列的fasta文件：

>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1 
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP

在这里我们读取fasta文件需要用到 strip 函数，其介绍为：
strip 用来去除头尾字符、空白符(包括\n、\r、\t、' '，即：换行、回车、制表符、空格)
lstrip 用来去除开头字符、空白符(包括\n、\r、\t、' '，即：换行、回车、制表符、空格)
rstrip 用来去除结尾字符、空白符(包括\n、\r、\t、' '，即：换行、回车、制表符、空格)
所以我们可以用下面的函数来读取fasta文件：

def load_fa(path):
    """a function to read fasta file from the path and store in a dict"""
    genes_seq = {}  #将序列存入字典
    with open(path,"r") as sequences:  #以读取方式打开文件
        lines = sequences.readlines()

    for line in lines:
        if line.startswith(">"):
            genename = line.split()[1]  #这个地方需要灵活调整
            genes_seq[genename] = ''  #序列为字符串
        else:
            genes_seq[genename] += line.strip()

    return genes_seq

计算序列k-mer

然后我们就需要对序列的k-mer进行计算，找出一段序列里面有哪些k-mer：

def build_kmers(seq, k_size):
    """a function to calculate kmers from seq"""
    kmers = []  # k-mer存储在列表中
    n_kmers = len(seq) - k_size + 1
    
    for i in range(n_kmers):
        kmer = seq[i:i + k_size]
        kmers.append(kmer)
        
    return kmers

对序列k-mer进行统计

在这一步，我们想统计出每段序列的各个k-mer的个数，需要使用到 collections 库中的 Counter 函数。

from collections import Counter
def summary_kmers(kmers):
    """a function to summarize the kmers"""
    kmers_stat = dict(Counter(kmers))
    return kmers_stat

实战

下面以这段fasta文件为例，分析其中基因的k-mer。

>chr18_4399328-4400910:+ ENSMUSG00000117579.2
CTGCAGTGTCAAGTCCCCCAATTATTTACTTTGCCTGTAAACGTGATGCACCCACCCACC
TTTACTTACTCACTGGCTCTAAAGGTGAAATTTGCTTTTTTCTTTTTCTTTTTAATGTCC
>chr18_4634878-4682869:+ ENSMUSG00000033960.7
GCGGTGGGCGGGACTGTGCGGGGCGGAGGGCGGGGCATGCGAGTGTGCTCCGAGCATGCT
CCACCCGGTAGTAGCAGGCTGGGTGGCTCGTGCTGGTCCCCGCCGGGCGAAGCGGCAGCG

genes_seq = load_fa(path="test.fa")
genes_kmers = {}
for gene in genes_seq.keys():
  genes_kmers[gene] = summary_kmers(build_kmers(seq=cluster1_genes_seq[gene], k_size=6))

通过上面这段代码，我们实际上得到了一个字典，字典的内容大概是这样的：

{ENSMUSG00000117579.2:{"ATCCGG":1,"ATCATT":3,...}}

一个字典种嵌套字典的结构，最外面字典的键为基因的ID，值为每个基因的k-mer信息，而每个基因的k-mer信息也是一个字典，它的键为k-mer序列，值为k-mer序列出现的频率。

当然，这肯定可读性不高，也不利于我们后续的数据处理，所以我们利用 pandas 库来对它进行变形：

import pandas as pd
pd.DataFrame(genes_kmers)

这样就会生成一个列为基因，行为k-mer的数据框，格式如图所示：

提前祝大家周末愉快~