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

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的数据框,格式如图所示:

提前祝大家周末愉快~

网友评论