手把手教你“破译”武汉新型冠状病毒（一天完成中文期刊《生物信息学》文章的工作）
手把手教你“破译”武汉新型冠状病毒（一天完成中文期刊《生物信息学》文章的工作）
写这篇“教程”的初衷是高教授发布的一篇文章，文章声称破译了2019新型冠状病毒。拜读后，深感科研文章用语之精妙：可变翻译，南开回文序列（Nankai complemented palindrome）、Nankai CDS、"分子功能与进化分析相结合的研究思想"，看起来都无比高大上。
本文就根据高教授在ResearchGate上公开的中文文章（以下简称“文章”），看一看此文做了什么工作，如何得到这些结果，并顺便根据此文之思想“发现”一些武汉新型冠状病毒基因序列中的信息。
注：本文实事上是一篇讽刺文，文章所写代码、步骤总用时不足一天，代码步骤都较为粗糙（也未写代码注释），但大致可还原高教授抢发的水文中的大部分实验结果。高教授在知乎文章中称“请大家支持中文期刊《生物信息学》”，我认为，支持中文期刊不是靠高声疾呼、高调灌水，不是大呼“给我上！”，而是实打实将高水平文章发在中文期刊上，以“跟我上！”之姿做表率，这才是真正的支持。
以下内容如您无兴趣皆可忽略，本文核心思想在上面的注释中都以表明。
1. 文章解读
图1 武汉 2019 冠状病毒的起源与可变翻译（CHEN Jiayuan, et al. 2020） 文章仅此一图。图1A展示的是不同病毒的“Nankai CDS”，作者称此区域包含一段22bp的互补回文序列，实际上可以看到，在部分碱基点突变后此区域并不“回文”；图1B展示的是不同病毒的CDS区序列去除图1A中8-11bp的“可变区”后根据序列比对结果构建的进化树，作者认为该方法“可以简化进化分析、提高信噪比，使分析结果更可靠，而且可以用于指导功能实验”（我感觉这就是一个CDS区编码蛋白的进化树……）；图1C展示的是Nankai CDS区可能的编码蛋白，称之为“可变翻译”，并认为可变翻译可以使病毒在不同宿主中表达不同的蛋白以适应环境（我认为简单的说就是病毒容易突变）。
2. 实验方法
2.1 寻找 Nankai回文序列
文章介绍，在前期工作中，他们在SARS冠状病毒中（DQ497008.1）发现了一段22 bp的互补回文序列（命名为Nankai complemented palindrome），并将其所在的CDS命名为Nankai CDS（见图1A）。
我们从GeneBank数据库中下载这一病毒的全基因组序列（DQ497008.1），保存为fasta文件；读取该文件并寻找最长的互补回文序列：
import re
fasta = {}
with open('DQ497008.1.fasta') as file:
sequence = ""
for line in file:
if line.startswith(">"):
name = line[1:].rstrip()
fasta[name] = ''
continue
fasta[name] += line.rstrip().upper()
def complement(s):
basecomplement = {
"A":"T",
"T":"A",
"G":"C",
"C":"G",
}
letters = list(s)
letters = [basecomplement[base] for base in letters]
return ''.join(letters)
maxlen = 0
for seqname,seq in fasta.items():
for i in range(len(seq)):
for j in range(10,50):
selectseq = seq[i:i+j]
reverse_complement_seq = complement(selectseq)[::-1]
if selectseq == reverse_complement_seq:
if len(selectseq) > maxlen:
maxlen = len(selectseq)
start = i+1
end = i+j+1
CPSeq = selectseq
print(start,end,CPSeq)
print(maxlen)
运行结果如下：
25962 25984 TCTTTAACAAGCTTGTTAAAGA
22
即该SARS病毒序列的25962-25984碱基位置是一段22bp长度的互补回文序列（TCTTTAACAAGCTTGTTAAAGA），与图1A一致。返回GeneBank，可以看到这段序列位于ORF 3b读码框，与文章描述一致。
2.2 构建进化树
作者在文章中是如此描述的：
进化树构建使用 13 条去除可变区的 Nankai CDS
简单的说就是使用各序列的ORF 3 CDS区进行局部比对，构建了一个进化树，并称之为：
当前大部分的冠状病毒基因组研究都是简单使用全基因组或某个病毒结构基因的序列， 后一种做法主观性较大。简单使用全基因组序列进行进化分析和重要功能或致病位点分析，不仅工作量巨大，而且大量与研究目的无关的基因组区域会对数据分析产生干扰。因此，不同于简单使用全基因组序列的其它研究，我们于2018年在国际上首次提出分子功能与进化分析相结合的研究思想， 并应用于冠状病毒基因组的研究。
这部分，我们可以使用NCBI的 blastn 工具，查询序列（Query Sequence）为 DQ497008.1，查询范围（Query subrange）是上述“Nankai CDS”，即25689-26153，程序优化选择 More dissimilar sequences；注意在算法参数（Algorithm parameters）中将显示结果由100改为20000，点击“BLAST”运行。网页展示运行结果后，进一步筛选出Beta冠状病毒（Betacoronavirus (taxid:694002)），结果如图2所示。
图2 部分检索结果，对号对应序列即为文章所用序列 这里略有疑惑的是，文章中表示
提取 NCBI GenBank 数据库中全部 beta 冠状病毒的完整基因组 （Complete genome）序列（以下简称病毒序列或序列）；去除与 Nankai CDS 相似度在 0.999 以上的序列后，留下 11 条病毒序列
而我在限定了物种为Beta冠状病毒（Betacoronavirus）并去除相似度在 0.999 以上的序列后，余下序列远不止11条。在这里，我从所有符合条件的序列中手动选出了文章中所说的10条序列（未添加此次的新型冠状病毒序列），使用blast提供的Blast Tree View查看进化树，结果如图3所示。与图1B基本一致。
此外，blast可以查看比对结果，逐个整理后便是图1A。
2.3 寻找“可变翻译”
恕我愚钝，并不能理解高教授在评论区一遍遍强调的“可变翻译”到底为何物。从文章中看，可变翻译指的是不同病毒的Nankai CDS区所预测翻译出的蛋白质。同样以SARS病毒（DQ497008.1）为例，导入fasta格式的序列进行分析：
import re
fasta = {}
with open('DQ497008.1.fasta') as file:
sequence = ""
for line in file:
if line.startswith(">"):
name = line[1:].rstrip()
fasta[name] = ''
continue
fasta[name] += line.rstrip().upper()
for seqname,seq in fasta.items():
useseq = seq[25688:26152] #选择CDS区序列
def trans(s):
transbase = {
'TCA': 'S', # Serine
'TCC': 'S', # Serine
'TCG': 'S', # Serine
'TCT': 'S', # Serine
'TTC': 'F', # Phenylalanine
'TTT': 'F', # Phenylalanine
'TTA': 'L', # Leucine
'TTG': 'L', # Leucine
'TAC': 'Y', # Tyrosine
'TAT': 'Y', # Tyrosine
'TAA': '*', # Stop
'TAG': '*', # Stop
'TGC': 'C', # Cysteine
'TGT': 'C', # Cysteine
'TGA': '*', # Stop
'TGG': 'W', # Tryptophan
'CTA': 'L', # Leucine
'CTC': 'L', # Leucine
'CTG': 'L', # Leucine
'CTT': 'L', # Leucine
'CCA': 'P', # Proline
'CCC': 'P', # Proline
'CCG': 'P', # Proline
'CCT': 'P', # Proline
'CAC': 'H', # Histidine
'CAT': 'H', # Histidine
'CAA': 'Q', # Glutamine
'CAG': 'Q', # Glutamine
'CGA': 'R', # Arginine
'CGC': 'R', # Arginine
'CGG': 'R', # Arginine
'CGT': 'R', # Arginine
'ATA': 'I', # Isoleucine
'ATC': 'I', # Isoleucine
'ATT': 'I', # Isoleucine
'ATG': 'M', # Methionine
'ACA': 'T', # Threonine
'ACC': 'T', # Threonine
'ACG': 'T', # Threonine
'ACT': 'T', # Threonine
'AAC': 'N', # Asparagine
'AAT': 'N', # Asparagine
'AAA': 'K', # Lysine
'AAG': 'K', # Lysine
'AGC': 'S', # Serine
'AGT': 'S', # Serine
'AGA': 'R', # Arginine
'AGG': 'R', # Arginine
'GTA': 'V', # Valine
'GTC': 'V', # Valine
'GTG': 'V', # Valine
'GTT': 'V', # Valine
'GCA': 'A', # Alanine
'GCC': 'A', # Alanine
'GCG': 'A', # Alanine
'GCT': 'A', # Alanine
'GAC': 'D', # Aspartic Acid
'GAT': 'D', # Aspartic Acid
'GAA': 'E', # Glutamic Acid
'GAG': 'E', # Glutamic Acid
'GGA': 'G', # Glycine
'GGC': 'G', # Glycine
'GGG': 'G', # Glycine
'GGT': 'G', # Glycine
}
for key,value in transbase.items():
s = s.replace(key,value)
return s
for seqname,seq in fasta.items():
useseq = seq[25688:26152]
transeq = []
dictrans = {}
for i in range(3):
for j in range(0,len(useseq),3):
threebases = useseq[i+j:i+j+3]
if len(threebases) == 3:
transeq.append(trans(threebases))
transeqname = str(i)+ '.' + seqname
dictrans[transeqname]= ''.join(transeq)
transeq.clear()
for i in range(3):
for j in range(0,len(useseq),3):
reverse_seq = useseq[::-1]
threebases = reverse_seq[i+j:i+j+3]
if len(threebases) == 3:
transeq.append(trans(threebases))
transeqname = '-' + str(i)+ '.' + seqname
dictrans[transeqname]= ''.join(transeq)
transeq.clear()
print(dictrans)
运行后可以看到DQ497008.1的“Nankai CDS”区六框翻译所得的蛋白质序列为
{'0.DQ497008.1 SARS coronavirus strain MA-15, complete genome': 'MMPTTLFAGTHITMTTVYHITVSQIQLSLLKVTAFQHQNSKKTTKLVVILRIGTQVLKTMSLYMAISPKFTTSLSLHKLLQTLVLKMLHSSSLTSLLKTHRMCKYTQSTALQELLIQQWIQFMMSRRRLLACLCKHKKVSTNLCTHSFRKKQVR', '1.DQ497008.1 SARS coronavirus strain MA-15, complete genome': '*CQLLCLLAHT*L*LLYTI*QCHRYNCRY*R*RHFNTKTQRRLPNWWLF*G*ALRC*RLCRCTWLFHRSLLPA*VYTNYYRHWY*KCYILHL*QAC*RPTECANTHNRRLFRSC*SSNGSNL**ADDDY*RAFVSTRK*VRTYVLIRFGRNRYV', '2.DQ497008.1 SARS coronavirus strain MA-15, complete genome': 'DANYFVCWHTHNYDYCIPYNSVTDTIVVTEGDGISTPKLKEDYQIGGYSEDRHSGVKDYVVVHGYFTEVYYQLESTQITTDTGIENATFFIFNKLVKDPPNVQIHTIDGSSGVANPAMDPIYDEPTTTTSVPL*AQESEYELMYSFVSEETGTL', '-0.DQ497008.1 SARS coronavirus strain MA-15, complete genome': 'IAWTKKALLTHVFKHE*KNTNVSVRSSAAAE*YLT*VTT*SLRTSRQLTHINV*ATQKLFEQFLLLTS*KLWSQTSLNTSEFDHHLKPLYRYMLLYQKLWTHG*ESYWWLNHQKKLKTTTLRQWKSLLLT*TL*QYTICHQYQYTHGRLFHQP*', '-1.DQ497008.1 SARS coronavirus strain MA-15, complete genome': 'LHGQRRLCLLMYSSMSERTRMFPCDHQQQPSSI*PR*RPNR*GLLGS*HT*TCKPPRNCSNNFYFLHRKSYGHRHH*THLSSTII*SHFIGTCCCIRNCGLTDRSLIGG*TIRRNSKPQLYGSGSHCC*HRHCDNIPYVISINTHTVVCFINRS', '-2.DQ497008.1 SARS coronavirus strain MA-15, complete genome': 'CMDKEGFAYSCIQA*VKEHECFRAIISSSRVVFNLGNDLIVEDFSAANTHKRVSHPEIVRTISTSYIVKVMVTDIIKHI*VRPSFEATLSVHVAVSEIVDSRIGVLLVVKPSEETQNHNFTAVEVIAVNIDTVTIYHMSSVSIHTRSFVSSTVV'}
进一步筛选
List = []
for seqname,seq in dictrans.items():
hit = re.finditer('M\w+\\*',seq)
for match in hit:
print(seqname)
print(match)
可以看到预测翻译所得的蛋白质及对应的翻译方法、翻译的起始位置。
2.DQ497008.1 SARS coronavirus strain MA-15, complete genome
<re.Match object; span=(117, 134), match='MDPIYDEPTTTTSVPL*'>
-0.DQ497008.1 SARS coronavirus strain MA-15, complete genome
<re.Match object; span=(92, 104), match='MLLYQKLWTHG*'>
-1.DQ497008.1 SARS coronavirus strain MA-15, complete genome
<re.Match object; span=(10, 34), match='MYSSMSERTRMFPCDHQQQPSSI*'>
-2.DQ497008.1 SARS coronavirus strain MA-15, complete genome
<re.Match object; span=(1, 15), match='MDKEGFAYSCIQA*'>
-2.DQ497008.1 SARS coronavirus strain MA-15, complete genome
<re.Match object; span=(70, 80), match='MVTDIIKHI*'>
输入11条序列的CDS区运行，整理后应该就是图1C。
3. 结语
至此，一篇中文核心期刊文章的大部分工作便介绍完了。此外，如果将2019新型冠状病毒（MN908947.3）的序列文件输入到本文2.1中的代码里运行，可以发现该病毒中存在一条20bp长的互补回文序列（ACACTGGTAATTACCAGTGT），位于5745-5765bp，我在这里非正式宣布其为Xiyang complemented palindrome，并将其所在的开放阅读框命名为Xiyang ORF。你还可以套用上述代码“破译”更多的病毒（最好是热点病毒），然后跟导师商量一下，挂个名就可以投中文核心试试了。（大部分情况下，你导师肯定会拒绝你！）
高教授的这篇文章槽点实在太多，不过此文的方法都很经典，比如我在这篇“教程”中用到的互补回文序列筛选、六框翻译等，倒是可以作为生物信息学课程的实践内容。此外，敲完这篇文章我也深深感到，即使灌水文章也不应太水，否则背黑锅担骂名的还是自己的导师。
本文可能有诸多不足及疏漏，若您有任何批评或疑问，欢迎在评论区留言。
后记
本教程里面出现了大量python自编脚本，如果你看不懂呀也不用着急哈，现在很多小学生都在学了，你可以向他们请教哦。反正现在疫情那么严重，大家都滞留在老家，充分沟通！