新冠病毒变异株核酸检测引物设计——代码实现1

     上一期我们对新冠变异株核酸检测的方法原理进行了简介，本期我们主要讲一下，如何利用python语言实现我们特异性检测的目标。

一、文件读取

    接着上一期我们下载新冠Alpha、Beta、Delta、Gamma株的全部基因组序列后，我们发现同一株型的新冠基因组，以“>”分开不同测序基因组，序列总共有上千条，为了后面程序简便，我们将同一株型的新冠基因组文件读取并储存数组中，如Alpha株，其数组形式为：Alpha[“序列1”，“序列2”，“序列3”，......，“序列n”]，其余毒株方法相同，代码如下：

Alpha = []
i = 0
a = ""
f = open("/home/lxh/Alpha.fasta","r")
for line in f.readlines():
  if ">" in line:
    i += 1
    Alpha.append(a)
    a = ""
  if ">" not in line:
    a += line.replace("\n","")

del Alpha[0]

f.close()

     简单说明下，上面代码意思是逐行读取，当遇到“>”建立新的数组元素，没有遇到的时候也就是将每一行核酸序列末尾的换行符“\n”去掉，加在一起变成一个字符串，这样一个字符串成为毒株数组中的一个元素，以此类推，完成所有基因组的读取和写入。

二、引物设计

      本次引物长度暂定为20bp，如要设计Alpha毒株的引物，那么引物必须满足能够同时出现Alpha数组的所有元素中，但不出现在其他毒株元素中。我们参考宏基因组测序kermer的分析方法。假设有30个序列，如果引物长度为20bp，那么会出现30-20+1=11条连续而且大概率不相同的引物。代码如下：

Alpha_kermer = []
for m in range(i-1):
  n = len(Alpha[m])
  for mm in range(n-19):
    Alpha_kermer.append(Alpha[m][mm:(20 + mm)])

       注意循环语句最后一句一定要从头空一行，表示循环语句从这里结束，新手一定注意，不懂的可以搜索WC3School，里面有各组编程语言的语法教程。

        以上代码将Alpha所有的基因组裁剪并形成20bp的数组，因为这些元素来自于同一株的新冠基因组，所以理论上就会出现大量相同序列的元素，下面的代码将数组中元素进行了频次统计，并进行元素出现频次高低的排序。新形成的bb成为一个无重复，频率由高到低出现的数组，通过检验bb的第一个元素bb[0]为一串polyA，果断放弃，第二个元素bb[1]为正常序列，出现频率为3879次，而我们上面写入的新冠基因组也有3879次，这说明此序列出现在每个基因组中。代码如下：

from collections import Counter
collection_words = Counter(Alpha_kermer)
i

bb = sorted(collection_words.items(),key = lambda x:x[1],reverse = True)

     出现3879次的序列远不止这一条，我们将出现频次为3879次的所有序列检索并形成新的数组，为了减少意外情况出现和重复计算，我们将结果导出保存，以便下次直接读取应用，代码如下：

f = open("Alpha_kermer.txt","a+")
fk = open("Alpha_kermer_all_in.txt","a+")
k_number = 0
for key,value in bb:
  f.write(key + "\n")
  if value == (i - 1):
    k_number += 1
    fk.write(key + "\n")

f.close()
fk.close()
k_number

        文件"Alpha_kermer.txt"为出现在Alpha中长度为20bp的所有引物序列，文件Alpha_kermer_all_in.txt为同时出现在所有Alpha株中长度为20bp的引物序列，两者为包含关系：“或”、“和”关系。

       其他株的代码，只需要查找和替代Alpha，换成其他株名就可以一键copy and run。