如何使用Python在指定目录中计算并筛选所有DNA序列的Kmer？

本文共计337个文字，预计阅读时间需要2分钟。

背景：Kmer是基因组组装算法中经常接触到的概念，简单来说，Kmer就是长度为k的核苷酸序列。

内容：一种长为m的reads可以分解成m-k+1个Kmer。Kmer的长度和阈值直接影响组装的效果。+Denovo组装流程。

背景

Kmer是基因组组装算法中经常接触到的概念，简单来说，Kmer就是长度为k的核苷酸序列。一般长短为m的reads可以分成m-k+1个Kmer。Kmer的长度和阈值直接影响到组装的效果。

Denovo组装流程：​​原始数据——数据过滤——纠错——kmer分析——denovo组装​​。

组装测序策略：根据基因组大小和具体情况选择个大概的k值，构建contig所需的数据量以及所需的构建的文库数量。对于植物基因组一般考虑的是大kmer（>31），动物一般在27左右，具体根据基因组情况调整。需要在短片段数据量达到20X左右的时候进行kmer分析。Kmer分析正常后，继续加测数据以达到最后期望的数据量。

编码

importos
importsys

#convertcommandlineargumentstovariables
kmer_size=int(sys.argv[1])
count_cutoff=int(sys.argv[2])

#definethefunctiontosplitdna
defsplit_dna(dna,kmer_size):
kmers=[]
forstartinrange(0,len(dna)-(kmer_size-1),1):
kmer=dna[start:start+kmer_size]
kmers.append(kmer)
returnkmers

#createanemptydictionarytoholdthecounts
kmer_counts={}

#processeachfilewiththerightname
forfile_nameinos.listdir("."):
iffile_name.endswith(".dna"):
dna_file=open(file_name)

#processeachDNAsequenceinafile
forlineindna_file:
dna=line.rstrip("\n")

#increasethecountforeachk-merthatwefind
forkmerinsplit_dna(dna,kmer_size):
current_count=kmer_counts.get(kmer,0)
new_count=current_count+1
kmer_counts[kmer]=new_count

#printk-merswhosecountsareabovethecutoff
forkmer,countinkmer_counts.items():
ifcount>count_cutoff:
print(kmer+":"+str(count))

作者：Bioinfarmer

若要及时了解动态信息，请关注同名微信公众号：Bioinfarmer。