新发变异株奥密克戎（Omicron）变异位点的获得及其核酸检测的引物设计

        

      继2021年11月9号首次在南非报道了第一例Omicron病例，短短3个月内，南非豪登省出现该病例激增，且南非所有省份的Omicron变异毒株病例“似乎都在增加”。为应对可能到来的新发变异毒株，核酸检测是必不可少的工具，特异性的检测依赖特异性的引物；氨基酸变异位点是疫苗设计和抗体检测的依据，因此本章我们有两个任务，一个是设计特异性检测Omicron的引物，一个是找出Omicron突变位点。

一、核酸检测引物设计

         引物获得方法如前，由于NCBI上Omicron的基因组数据量暂时很少，只有18条，高质量拼接的基因组只有1条。为了获得较多的分析数据，我们把含有特殊符号“N‘’的基因组也全部用上了，但“N”可能代表"A/G/C/T"里面的任何一位，是无法分析的，因此我们要去掉含有“N‘’的kermer。

Omicron_kermer = []
for m in range(len(Omicron) - 1):
  n = len(Omicron[m])
  for mm in range(n-20):
    a = "".join(Omicron[m][mm:(20 + mm)])
      if "N" not in a:
        Omicron0_kermer.append(Omicron0[0][mm:(20 + mm)])

           下图为我们最终获得的Omicron的kermer结果，过程代码我们将不再描述，感兴趣的可以查阅上一期文章新冠病毒变异株核酸检测引物设计——代码实现2“”。

        同样为了特异性检测Omicro，我们必须将出现在Omicro中的kermer在其他新冠变异株的kermer中比对一下，结果如下：

           我们调整了Omicro的kermer出现频率为大于90%，其他株出现频率为0%，代码如下：

f9_1 = open("Omicron_primer.txt","a+")
f9_1.write("O_kermer" + "\t" + "per_O_kermer" + "\t" + "per_A_kermer" + "\t" + "per_B_kermer" + "\t" + "per_G_kermer" + "\t" + "per_D_kermer" + "\n")
for i in range(len(Omicron) - 1):
  Omicron_little = Omicron[i].split("\t")
  if (float(Omicron_little[1].replace("%","")) > 90) & (float(Omicron_little[2].replace("%","")) < 5) & (float(Omicron_little[3].replace("%","")) < 5) & (float(Omicron_little[4].replace("%","")) < 5):
    if float(Omicron_little[5].replace("%","")) < 5:
      f9_1.write(Omicron[i] + "\n")

f9_1.close()

         结果如下：

          由于Omicro中出现连续A/G/C/T较多，所以我们放宽了筛选条件，只要求一条序列满足Tm值以及无连续的A/G/C/T，其可以作为前或者后qpcr引物。

f_primer = open("Omicron_primer.txt","a+")
NN = len(Omicron_num)
for e in range(NN - 1):
  mgb = ''.join(Omicron0_kermer[Omicron_num[e]])
  A = mgb.count("A")
  G = mgb.count("G")
  C = mgb.count("C")
  T = mgb.count("T")
  Tm1 = 4 * (G + C) + 2 * (A + T)
  if (Tm1 < 70) & (Tm1 > 50) & ("AAA" not in mgb) & ("GGG" not in mgb) & ("CCC" not in mgb) & ("TTT" not in mgb) & ((mgb[0] == "A") or (mgb[0] == "T")) & ((mgb[19] == "C") or (mgb[19] == "G")):
    f_primer.write("<" + Omicron0_kermer[Omicron_num[e]] + ">" + "\n" + str(number) + "\n")

f_primer.close()

  结果如下，我们得到了11条候选的序列，后面可以根据primer5观察其键能、自我配对情况、以及与自设计的另外一条引物的互补配对情况，筛选最合适的引物序列，在这里我们不做过多的介绍，后面有机会可以写写这方面的知识。

<AAGTACTTAATGAGAGGTGC>
<ACTTAATGAGAGGTGCTCTG>
<TTAATGAGAGGTGCTCTGCC>
<ACCAATGTGCTATGAGGCAC>
<ACTCAGACTAAGTCTCATCG>
<TCAGACTAAGTCTCATCGGC>
<AGACTAAGTCTCATCGGCGG>
<TAGCCATGGCAGGTTCCAAC>
<AGCCATGGCAGGTTCCAACG>
<ATGGCAGGTTCCAACGGTAC>
<TTACTAATTATTATGCGGAC>

         到这里，Omicro的核酸检测引物我们算基本完成了。

二、突变核苷酸查找

      首先基于相同kermer（20bp）片段的频率统计，决定了任意一个突变位点都会出现至少20个连续索引号的含有突变位点的kermer，如果两个或者两个以上突变位点相距小于等于20个核苷酸，那么出现连续索引号的kermer将大于20个，我们可以预见相邻两个小于20bp的突变位点至少产生20+n个连续索引号的kermer，其中n为第一个突变位点到最后一个突变位点（连续相邻突变位点小于等于20）的距离。

          基于以上事实，我们统计单独在Omicron中高频率出现，而在其他毒株中低频率出现的kermer的索引号，我们以Omicron的参考基因组作为索引序号来源。

f9_1 = open("Omicron_little_in_Alpha_Beta_Gamma_Delta.txt","r")
Omicron_only = []
ii = 0
for i in f9_1.readlines():
  ii += 1
  if ii > 1:
    Omicron_little = i.replace("\n","").split("\t")
    Omicron_only.append(Omicron_little[0])

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

del Omicron0[0]
f.close()
Omicron0_kermer = []
for mm in range(len(Omicron0[0]) - 20):
  Omicron0_kermer.append(Omicron0[0][mm:(20 + mm)])

index_Omicron = {} 
for i in range(AA - 1):
  index = Omicron0_kermer.index(Omicron_only[i])
  index_Omicron[Omicron_only[i]] = index

index_sort_A = sorted(index_Omicron.items(),key = lambda x:x[1],reverse = False)

Omicron_num = []
for x,y in index_sort_A:
  Omicron_num.append(y)

  上图即为Omicron对参考基因组的索引号，并做了排序，结果我们可以看出，连续数字较多，符合我们上述推断，我们对连续的索引数字进行整理排序，即可推断出突变位点在什么地方以及突变后的核苷酸是什么。下面为整理Omicron连续索引号的代码。

f = open("/home/lxh/Omicron_mutation_index.fasta","a+")
f.write("continuous_index" + "\n" + str(Omicron_num[0]) + "-")
for i in range(len(Omicron_num) - 1):
  if Omicron_num[i] + 1 != Omicron_num[i + 1]:
    f.write(str(Omicron_num[i]) + ":" + Omicron0_kermer[Omicron_num[i]] + "\n" + str(Omicron_num[i + 1]) + "-")

f.write(str(Omicron_num[- 1]) + ":" + Omicron0_kermer[Omicron_num[- 1]])
f.close()

         找到突变序列后，我们需要更进一步找到突变位点的位置和突变的核苷酸。下面以Alpha突变位点为例，连续kermer的最后一个kermer突变位点在第一位，因此减去第一位核苷酸后的序列在其他毒株中应该是高频出现，并且其序列排列方式如同这个kermer一样，也是除第一位不同外，其他序列完全相同，基于这一思路，完成这个操作的代码如下：

f9 = open("Alpha_mutation_index.txt","r")
i = 0
Alpha_mutation = []
for line in f9.readlines():
  i += 1
  if i > 1:
    Alpha_mutation.append(line.replace("\n",""))

f9.close()
f10 = open("Beta_kermer.txt","r")
Beta_kermer = []
mm = 0
for line in f10.readlines():
  mm += 1
  if mm > 1:
    Beta_kermer.append(line.replace("\n",""))

f10.close()
f11 = open("Alpha_mutation.txt","a+")
for line1 in range(len(Alpha_mutation) - 1):
  b = "".join(Alpha_mutation[line1]).split(":")
  a = "".join(b[0]).split("-")
  for line in range(len(Beta_kermer) - 1):
    aa = "".join(Beta_kermer[line]).split("\t")
    if (aa[0][1:19] == b[1][1:19]) & (float("".join(aa[5]).replace("%","")) > 80):
      f11.write(Alpha_mutation[line1] + ":" + aa[0][0] + str(int(a[1]) + 1) + b[1][0] + "\n")

f11.close()

  结果如下，对下图做一番简单的解释，889-908表示以长度为20bp的kermer连续索引号从889开始908结束，并且相隔正好20个核苷酸，这表明索引号889的最后一位核苷酸正好是索引号为908的第一个核苷酸，并且这个核苷酸即为突变位点。第一个冒号（：）后面的20bp的核苷酸为索引号908的kermer序列。第二个冒号（：）后面的C909T表示909位核苷酸由原来的胞嘧啶（C）突变为现在的胸腺嘧啶（T）。

   找到突变核苷酸位点后，进一步我们需要找到其突变氨基酸，首先根据20bp的kermer我们很容易锁定突变发生在基因组上的位置，其次，我们根据突变位点的索引号我们可以判断其突变位点位于三联密码子的第几位，如第一行结果C909T，突变位点为第909个核苷酸，除以3，余数为0，这说明这个核苷酸位于三联密码子的最三位，如果余数为1，表明这个核苷酸位于三联密码子的第一位，如果余数为2，表明这个核苷酸位于三联密码子的最二位。只需要简单对照密码子表即可知道突变前后氨基酸的变化（将表中U直接当成T就行）。