什么是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一共有 4k 种。
分析步骤
在这里我们的分析步骤包括：
- 读取序列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的数据框，格式如图所示：
提前祝大家周末愉快~
网友评论