Case 3. mRNA 疫苗设计¶
1. 获取序列¶
SARS-CoV-2 Spike NCBI Refseq and CDS fasta.
In [5]:
import requests
from pathlib import Path
PROXY = "http://127.0.0.1:7897"
NCBI_EFETCH = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
proxies = {
"http": PROXY,
"https": PROXY,
}
def download_fasta(params, output_file):
response = requests.get(
NCBI_EFETCH,
params=params,
proxies=proxies,
timeout=60,
)
response.raise_for_status()
fasta = response.text.strip()
if not fasta.startswith(">"):
raise RuntimeError(f"NCBI 返回内容异常：\n{fasta[:500]}")
Path(output_file).write_text(fasta + "\n", encoding="utf-8")
sequence = "".join(
line.strip()
for line in fasta.splitlines()
if not line.startswith(">")
)
print(f"文件：{output_file}")
print(f"序列长度：{len(sequence)} bp")
# SARS-CoV-2 Wuhan-Hu-1 RefSeq 完整基因组
download_fasta(
params={
"db": "nuccore",
"id": "NC_045512.2",
"rettype": "fasta",
"retmode": "text",
},
output_file="SARS-CoV-2_NC_045512.2_RefSeq.fasta",
)
# Spike CDS
# NC_045512.2 中 Spike CDS 坐标为 21563..25384
download_fasta(
params={
"db": "nuccore",
"id": "NC_045512.2",
"seq_start": 21563,
"seq_stop": 25384,
"rettype": "fasta",
"retmode": "text",
},
output_file="SARS-CoV-2_Spike_CDS.fasta",
)
# Spike 蛋白序列
download_fasta(
params={
"db": "protein",
"id": "YP_009724390.1",
"rettype": "fasta",
"retmode": "text",
},
output_file="SARS-CoV-2_Spike_protein.fasta",
)
# 打印说明行与长度
text = Path("SARS-CoV-2_Spike_protein.fasta").read_text(encoding="utf-8").strip()
header, _, seq = text.partition("\n")
seq = seq.replace("\n", "")
print(f"说明：{header}")
print(f"长度：{len(seq)} aa")
文件：SARS-CoV-2_NC_045512.2_RefSeq.fasta 序列长度：29903 bp 文件：SARS-CoV-2_Spike_CDS.fasta 序列长度：3822 bp 文件：SARS-CoV-2_Spike_protein.fasta 序列长度：1273 bp 说明：>YP_009724390.1 surface glycoprotein [Severe acute respiratory syndrome coronavirus 2] 长度：1273 aa
2. 预测Refseq DNA的蛋白产物¶
In [8]:
from Bio import SeqIO
# 本节为纯本地计算：读取已下载的 fasta，不联网、不使用代理
# Spike CDS 在 NC_045512.2 中的坐标（1-based，含端点）
CDS_START, CDS_STOP = 21563, 25384
genome = SeqIO.read("SARS-CoV-2_NC_045512.2_RefSeq.fasta", "fasta")
# 从 RefSeq DNA 中切出 CDS（坐标转 0-based 左闭右开）
cds = genome.seq[CDS_START - 1 : CDS_STOP]
# 按标准遗传密码表(table=1)翻译
protein = cds.translate(table=1, to_stop=False)
# 去掉末尾终止密码子对应的 "*"
protein = protein[:-1] if protein.endswith("*") else protein
print(f"CDS 长度：{len(cds)} bp")
print(f"预测蛋白长度：{len(protein)} aa")
print(f"蛋白序列：{protein}")
# 与 NCBI 注释蛋白 YP_009724390.1 比较
ref = SeqIO.read("SARS-CoV-2_Spike_protein.fasta", "fasta")
print(f"\n与 {ref.id} 一致：{str(protein) == str(ref.seq).rstrip('*')}")
# 保存预测的蛋白产物
protein_record = SeqIO.SeqRecord(
protein,
id=f"NC_045512.2:{CDS_START}-{CDS_STOP}",
description="predicted surface glycoprotein (translated from RefSeq CDS)",
)
_ = SeqIO.write(protein_record, "SARS-CoV-2_Spike_predicted_protein.fasta", "fasta")
print("已保存：SARS-CoV-2_Spike_predicted_protein.fasta")
CDS 长度：3822 bp 预测蛋白长度：1273 aa 蛋白序列：MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFSNVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIVNNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLEGKQGNFKNLREFVFKNIDGYFKIYSKHTPINLVRDLPQGFSALEPLVDLPIGINITRFQTLLALHRSYLTPGDSSSGWTAGAAAYYVGYLQPRTFLLKYNENGTITDAVDCALDPLSETKCTLKSFTVEKGIYQTSNFRVQPTESIVRFPNITNLCPFGEVFNATRFASVYAWNRKRISNCVADYSVLYNSASFSTFKCYGVSPTKLNDLCFTNVYADSFVIRGDEVRQIAPGQTGKIADYNYKLPDDFTGCVIAWNSNNLDSKVGGNYNYLYRLFRKSNLKPFERDISTEIYQAGSTPCNGVEGFNCYFPLQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVCGPKKSTNLVKNKCVNFNFNGLTGTGVLTESNKKFLPFQQFGRDIADTTDAVRDPQTLEILDITPCSFGGVSVITPGTNTSNQVAVLYQDVNCTEVPVAIHADQLTPTWRVYSTGSNVFQTRAGCLIGAEHVNNSYECDIPIGAGICASYQTQTNSPRRARSVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSMTKTSVDCTMYICGDSTECSNLLLQYGSFCTQLNRALTGIAVEQDKNTQEVFAQVKQIYKTPPIKDFGGFNFSQILPDPSKPSKRSFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVLPPLLTDEMIAQYTSALLAGTITSGWTFGAGAALQIPFAMQMAYRFNGIGVTQNVLYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVEAEVQIDRLITGRLQSLQTYVTQQLIRAAEIRASANLAATKMSECVLGQSKRVDFCGKGYHLMSFPQSAPHGVVFLHVTYVPAQEKNFTTAPAICHDGKAHFPREGVFVSNGTHWFVTQRNFYEPQIITTDNTFVSGNCDVVIGIVNNTVYDPLQPELDSFKEELDKYFKNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIKWPWYIWLGFIAGLIAIVMVTIMLCCMTSCCSCLKGCCSCGSCCKFDEDDSEPVLKGVKLHYT 与 YP_009724390.1 一致：True 已保存：SARS-CoV-2_Spike_predicted_protein.fasta
3. Coden Usage Analysis¶
分别计算 Human codon usage， SARS-CoV-2 codon usage，以及mRNA vaccines clean.txt中的mRNA疫苗序列 计算RSCU， CAI, GC%, GC3等指标
3.1 数据载入与 CDS 提取¶
读取 mRNA vaccines clean.txt 与 SARS-CoV-2 注释 CDS，自动识别每条核酸序列的编码区（CDS）：疫苗 mRNA 取最长完整 ORF，病毒基因组序列按 frame 0 判断是否有内部终止密码子。
In [9]:
# ================== 3.1 载入序列并提取编码区（CDS） ==================
import re
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import requests
from Bio import SeqIO
import python_codon_tables as pct
sns.set_theme(style="whitegrid", context="notebook")
plt.rcParams["figure.dpi"] = 120
plt.rcParams["font.sans-serif"] = ["Noto Sans CJK JP", "WenQuanYi Micro Hei", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False
VACCINE_FILE = "mRNA vaccines clean.txt" # 疫苗 mRNA 与 S 基因序列
VIRAL_CDS_FILE = "SARS-CoV-2_NC_045512.2_CDS_all.fasta" # 病毒全部注释 CDS
BASES = "TCAG"
CODONS = [a + b + c for a in BASES for b in BASES for c in BASES]
STOP_CODONS = {"TAA", "TAG", "TGA"}
def read_fasta(path):
"""轻量 FASTA 解析，统一为大写 DNA（U → T）"""
records, name = {}, None
for line in Path(path).read_text().splitlines():
line = line.strip()
if not line:
continue
if line.startswith(">"):
name = line[1:].split()[0]
records[name] = ""
else:
records[name] += line
return {k: v.upper().replace("U", "T") for k, v in records.items()}
def frame0_stops(seq):
"""返回第 1 读框（frame 0）中终止密码子的起始位置"""
return [i for i in range(0, len(seq) - len(seq) % 3, 3) if seq[i:i + 3] in STOP_CODONS]
def find_orf(seq, min_aa=100):
"""扫描所有 ATG，返回最长完整 ORF 的 (start, end)，end 含终止密码子"""
best = None
for i in range(len(seq) - 3):
if seq[i:i + 3] != "ATG":
continue
j = i
while j + 3 <= len(seq) and seq[j:j + 3] not in STOP_CODONS:
j += 3
if j + 3 > len(seq): # 读到最后也没有终止密码子
continue
if (j - i) // 3 >= min_aa and (best is None or j - i > best[1] - best[0]):
best = (i, j + 3)
return best
def extract_cds(seq, min_aa=100):
"""提取编码区：优先认定 frame 0（无内部终止密码子），否则搜索最长完整 ORF"""
if len(seq) % 3 == 0:
stops = frame0_stops(seq)
if all(i == len(seq) - 3 for i in stops): # [] 或只有末尾一个终止密码子
return seq, "frame 0 编码区（无内部终止密码子）"
orf = find_orf(seq, min_aa)
if orf:
return seq[orf[0]:orf[1]], f"最长 ORF {orf[0]}..{orf[1]}"
return seq, "未找到完整 ORF，按 frame 0 处理"
# ---------- 1) mRNA 疫苗 / S 基因序列 ----------
vaccine_raw = read_fasta(VACCINE_FILE)
vaccine_cds, vaccine_note = {}, {}
print(f"【{VACCINE_FILE}】共 {len(vaccine_raw)} 条序列")
print(f"{'序列名':<32s}{'全长(nt)':>10s}{'CDS(nt)':>10s}{'蛋白(aa)':>10s} 说明")
for name, seq in vaccine_raw.items():
cds_seq, note = extract_cds(seq)
vaccine_cds[name] = cds_seq
vaccine_note[name] = note
print(f"{name:<32s}{len(seq):>10d}{len(cds_seq):>10d}{(len(cds_seq) - 3) // 3:>10d} {note}")
# ---------- 2) SARS-CoV-2 全部注释 CDS（本地缺失时直连 NCBI 下载，不使用代理） ----------
if not Path(VIRAL_CDS_FILE).exists():
url = ("https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
"?db=nuccore&id=NC_045512.2&rettype=fasta_cds_na&retmode=text")
resp = requests.get(url, timeout=60) # 直连，不设置 proxies
resp.raise_for_status()
Path(VIRAL_CDS_FILE).write_text(resp.text, encoding="utf-8")
print(f"\n已下载：{VIRAL_CDS_FILE}")
# 同一基因可能有多条重叠记录（如 ORF1ab 与 ORF1a），每个基因只保留最长的一条
viral_cds = {}
for rec in SeqIO.parse(VIRAL_CDS_FILE, "fasta"):
m = re.search(r"\[gene=([^\]]+)\]", rec.description)
gene = m.group(1) if m else rec.id
seq = str(rec.seq).upper()
if gene not in viral_cds or len(seq) > len(viral_cds[gene]):
viral_cds[gene] = seq
print(f"\n【{VIRAL_CDS_FILE}】去重后共 {len(viral_cds)} 个基因")
print(f"{'基因':<10s}{'长度(nt)':>10s}{'蛋白(aa)':>10s}")
for gene, seq in viral_cds.items():
print(f"{gene:<10s}{len(seq):>10d}{(len(seq) - 3) // 3:>10d}")
print(f"{'合计':<10s}{sum(len(s) for s in viral_cds.values()):>10d}")
【mRNA vaccines clean.txt】共 5 条序列 序列名 全长(nt) CDS(nt) 蛋白(aa) 说明 WHO_BNT162b2 4284 3822 1273 最长 ORF 54..3876 ModernaMrna1273 4004 3822 1273 最长 ORF 57..3879 ReconstructedBNT162b2 4175 3822 1273 最长 ORF 54..3876 S 3780 3780 1259 frame 0 编码区（无内部终止密码子） wuhCor1_ncbiGeneBGP_S 3822 3822 1273 frame 0 编码区（无内部终止密码子） 【SARS-CoV-2_NC_045512.2_CDS_all.fasta】去重后共 11 个基因 基因 长度(nt) 蛋白(aa) ORF1ab 21291 7096 S 3822 1273 ORF3a 828 275 E 228 75 M 669 222 ORF6 186 61 ORF7a 366 121 ORF7b 132 43 ORF8 366 121 N 1260 419 ORF10 117 38 合计 29265
3.2 指标计算函数与人源参考密码子表¶
定义 RSCU、CAI、GC / GC3 / GC3s、ENC 的计算函数，并载入 Kazusa 人源密码子表（TaxID 9606，组内相对频率）作为 CAI 与「最优密码子」的参考。
In [10]:
# ================== 3.2 指标计算函数与人源参考密码子表 ==================
from Bio.Seq import Seq
# 标准遗传密码表（NCBI table 1）
CODON2AA = {c: str(Seq(c).translate(table=1)) for c in CODONS}
AA2CODONS = {}
for c, aa in CODON2AA.items():
AA2CODONS.setdefault(aa, []).append(c)
SENSE_AA = [aa for aa in AA2CODONS if aa != "*"] # 20 种氨基酸
def count_codons(seq):
"""按读框统计 64 个密码子的频次"""
counts = dict.fromkeys(CODONS, 0)
for i in range(0, len(seq) - len(seq) % 3, 3):
codon = seq[i:i + 3]
if codon in counts:
counts[codon] += 1
return counts
def merge_counts(seqs):
"""把多条序列的密码子频次相加"""
total = dict.fromkeys(CODONS, 0)
for seq in seqs:
for codon, n in count_codons(seq).items():
total[codon] += n
return total
def rscu_from_counts(counts):
"""RSCU：某密码子观测值 ÷ 该氨基酸同义密码子的平均观测值（终止密码子不计）"""
rscu = {}
for aa in SENSE_AA:
codons = AA2CODONS[aa]
s = sum(counts[c] for c in codons)
if s == 0:
continue
for c in codons:
rscu[c] = counts[c] / s * len(codons)
return rscu
def cai_weights(reference, floor=0.5):
"""参考表的 CAI 权重 w = RSCU / RSCU_max = 组内频率 / 最高组内频率"""
weights = {}
for aa, codons in reference.items():
if aa == "*":
continue
fmax = max(codons.values())
for c, f in codons.items():
weights[c] = max(f / fmax, floor) # 防止 ln(0)
return weights
def cai(seq, weights):
"""Sharp & Li (1987) CAI：对 w 取几何平均，不计 Met / Trp / 终止密码子"""
logs = []
for i in range(0, len(seq) - len(seq) % 3, 3):
codon = seq[i:i + 3]
if codon in weights and len(AA2CODONS[CODON2AA[codon]]) > 1:
logs.append(np.log(weights[codon]))
return float(np.exp(np.mean(logs))) if logs else np.nan
def gc_metrics(seq):
"""GC 总含量、三个密码子位置 GC%、GC3s（同义位点第三位 GC%）"""
codons = [seq[i:i + 3] for i in range(0, len(seq) - len(seq) % 3, 3)]
gc = lambda s: 100 * sum(ch in "GC" for ch in s) / len(s) if s else np.nan
gc3s = [c[2] for c in codons
if len(AA2CODONS[CODON2AA[c]]) > 1 and CODON2AA[c] != "*"]
return {
"GC%": gc(seq[:len(codons) * 3]),
"GC1%": gc("".join(c[0] for c in codons)),
"GC2%": gc("".join(c[1] for c in codons)),
"GC3%": gc("".join(c[2] for c in codons)),
"GC3s%": gc("".join(gc3s)),
}
def enc(counts):
"""Wright (1990) 有效密码子数 ENC = 2 + 9/F̄2 + 1/F̄3 + 5/F̄4 + 3/F̄6"""
F = {}
for aa in SENSE_AA:
codons = AA2CODONS[aa]
k = len(codons) # 简并度
n = sum(counts[c] for c in codons)
if k == 1 or n <= 1: # Met / Trp 及样本量不足的家族跳过
continue
p2 = sum((counts[c] / n) ** 2 for c in codons)
F.setdefault(k, []).append((n * p2 - 1) / (n - 1))
if len(F) < 4:
return np.nan
Fbar = {k: float(np.mean(v)) for k, v in F.items()}
return 2 + sum(w / Fbar[k] for k, w in {2: 9, 3: 1, 4: 5, 6: 3}.items() if k in Fbar)
# ---------- 人源参考密码子使用谱（Kazusa, Homo sapiens, TaxID 9606） ----------
human_table = pct.get_codons_table("h_sapiens_9606") # {氨基酸: {密码子: 组内相对频率}}
human_weights = cai_weights(human_table)
human_rscu = {c: f * len(AA2CODONS[aa]) # RSCU = 组内频率 × 简并度
for aa, codons in human_table.items() if aa != "*"
for c, f in codons.items()}
human_best = {aa: max(codons, key=codons.get) # 每个氨基酸的人源最优密码子
for aa, codons in human_table.items() if aa != "*"}
print("人源参考表示例（组内相对频率）:")
for aa in ["L", "A", "V", "S", "R"]:
print(f" {aa}: " + " ".join(f"{c}={human_table[aa][c]:.3f}" for c in sorted(human_table[aa])))
人源参考表示例（组内相对频率）: L: CTA=0.070 CTC=0.200 CTG=0.400 CTT=0.130 TTA=0.080 TTG=0.130 A: GCA=0.230 GCC=0.400 GCG=0.110 GCT=0.270 V: GTA=0.120 GTC=0.240 GTG=0.460 GTT=0.180 S: AGC=0.240 AGT=0.150 TCA=0.150 TCC=0.220 TCG=0.050 TCT=0.190 R: AGA=0.210 AGG=0.210 CGA=0.110 CGC=0.180 CGG=0.200 CGT=0.080
3.3 计算各分析单元的密码子使用指标¶
对 SARS-CoV-2 全部 CDS、S 基因、各条疫苗 CDS 分别计算 CAI（以人源表为参考）、人源最优密码子比例、ENC 及各位点 GC 含量，结果保存为 codon_usage_summary.csv。
In [11]:
# ================== 3.3 计算各分析单元的密码子使用指标 ==================
VACCINE_NAMES = ["WHO_BNT162b2", "ModernaMrna1273", "ReconstructedBNT162b2"]
NATIVE_S_NAME = "wuhCor1_ncbiGeneBGP_S"
# 分析单元：名称 -> 参与计算的编码序列
analysis_units = {
"SARS-CoV-2（全部 CDS）": list(viral_cds.values()),
"SARS-CoV-2（S 基因）": [viral_cds["S"]],
"疫苗 CDS（3 条合计）": [vaccine_cds[n] for n in VACCINE_NAMES],
"BNT162b2（WHO）": [vaccine_cds["WHO_BNT162b2"]],
"mRNA-1273（Moderna）": [vaccine_cds["ModernaMrna1273"]],
"ReconstructedBNT162b2": [vaccine_cds["ReconstructedBNT162b2"]],
"S 基因（原生态，未优化）": [vaccine_cds[NATIVE_S_NAME]],
}
def optimal_codon_fraction(seq, best=human_best):
"""人源最优密码子（组内 RSCU 最高者）使用比例（%），Met/Trp/终止不计"""
hit = total = 0
for i in range(0, len(seq) - len(seq) % 3, 3):
codon = seq[i:i + 3]
aa = CODON2AA[codon]
if aa == "*" or len(AA2CODONS[aa]) == 1:
continue
total += 1
hit += codon == best[aa]
return 100 * hit / total if total else np.nan
rows, rscu_tables, counts_tables = [], {}, {}
for label, seqs in analysis_units.items():
merged = "".join(seqs) # 同一单元内多条 CDS 首尾相连（长度都是 3 的倍数）
counts = merge_counts(seqs)
counts_tables[label] = counts
rscu_tables[label] = rscu_from_counts(counts)
rows.append({
"分析单元": label,
"CDS 数": len(seqs),
"CDS 长度(nt)": len(merged),
"密码子数": len(merged) // 3,
"CAI(人源)": cai(merged, human_weights),
"人源最优密码子%": optimal_codon_fraction(merged),
"ENC": enc(counts),
**gc_metrics(merged),
})
summary = pd.DataFrame(rows).set_index("分析单元").round(3)
summary.to_csv("codon_usage_summary.csv", encoding="utf-8-sig")
print("已保存：codon_usage_summary.csv\n")
summary
已保存：codon_usage_summary.csv
Out[11]:
| CDS 数 | CDS 长度(nt) | 密码子数 | CAI(人源) | 人源最优密码子% | ENC | GC% | GC1% | GC2% | GC3% | GC3s% | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 分析单元 | |||||||||||
| SARS-CoV-2（全部 CDS） | 11 | 29265 | 9755 | 0.725 | 20.995 | 45.384 | 37.926 | 46.899 | 38.596 | 28.283 | 25.886 |
| SARS-CoV-2（S 基因） | 1 | 3822 | 1274 | 0.730 | 19.727 | 44.159 | 37.310 | 45.290 | 39.953 | 26.688 | 25.180 |
| 疫苗 CDS（3 条合计） | 3 | 11466 | 3822 | 0.961 | 82.358 | 26.860 | 58.756 | 50.262 | 40.188 | 85.819 | 85.592 |
| BNT162b2（WHO） | 1 | 3822 | 1274 | 0.952 | 77.787 | 28.999 | 56.986 | 49.765 | 40.188 | 81.005 | 80.674 |
| mRNA-1273（Moderna） | 1 | 3822 | 1274 | 0.981 | 91.500 | 22.037 | 62.297 | 51.256 | 40.188 | 95.447 | 95.429 |
| ReconstructedBNT162b2 | 1 | 3822 | 1274 | 0.952 | 77.787 | 28.999 | 56.986 | 49.765 | 40.188 | 81.005 | 80.674 |
| S 基因（原生态，未优化） | 1 | 3822 | 1274 | 0.730 | 19.727 | 44.159 | 37.310 | 45.290 | 39.953 | 26.688 | 25.180 |
3.4 三类密码子使用谱（RSCU）对比¶
把 Human、SARS-CoV-2、疫苗 CDS 的 RSCU 汇总到同一张表，并计算两两相关系数，用于衡量密码子偏好相似度（保存 rscu_comparison.csv、rscu_correlation.csv）。
In [12]:
# ================== 3.4 三类密码子使用谱（RSCU）对比 ==================
SENSE_CODONS = [c for c in CODONS if c not in STOP_CODONS] # 61 个有义密码子
rscu_compare = pd.DataFrame({
"Human（Kazusa）": pd.Series(human_rscu),
"SARS-CoV-2": pd.Series(rscu_tables["SARS-CoV-2（全部 CDS）"]),
"疫苗 CDS": pd.Series(rscu_tables["疫苗 CDS（3 条合计）"]),
"S 基因（原生）": pd.Series(rscu_tables["SARS-CoV-2（S 基因）"]),
}).reindex(SENSE_CODONS)
rscu_compare.insert(0, "氨基酸", [CODON2AA[c] for c in rscu_compare.index])
rscu_compare.insert(1, "人源最优", ["★" if human_best[CODON2AA[c]] == c else "" for c in rscu_compare.index])
rscu_compare.index.name = "密码子"
rscu_compare = rscu_compare.round(3)
# RSCU 之间的相关性：越接近 1 说明密码子偏好越相似
corr = rscu_compare.drop(columns=["氨基酸", "人源最优"]).corr().round(3)
rscu_compare.to_csv("rscu_comparison.csv", encoding="utf-8-sig")
corr.to_csv("rscu_correlation.csv", encoding="utf-8-sig")
counts_tables["SARS-CoV-2（全部 CDS）"] # 供后续绘图/导出使用
print("已保存：rscu_comparison.csv、rscu_correlation.csv\n")
print("RSCU 相关系数（61 个有义密码子）:")
print(corr.to_string(), "\n")
print("RSCU 对比表（★ 标记人源最优密码子）:")
rscu_compare
已保存：rscu_comparison.csv、rscu_correlation.csv
RSCU 相关系数（61 个有义密码子）:
Human（Kazusa） SARS-CoV-2 疫苗 CDS S 基因（原生）
Human（Kazusa） 1.000 -0.160 0.817 -0.109
SARS-CoV-2 -0.160 1.000 -0.353 0.965
疫苗 CDS 0.817 -0.353 1.000 -0.341
S 基因（原生） -0.109 0.965 -0.341 1.000
RSCU 对比表（★ 标记人源最优密码子）:
Out[12]:
| 氨基酸 | 人源最优 | Human（Kazusa） | SARS-CoV-2 | 疫苗 CDS | S 基因（原生） | |
|---|---|---|---|---|---|---|
| 密码子 | ||||||
| TTT | F | 0.92 | 1.404 | 0.286 | 1.532 | |
| TTC | F | ★ | 1.08 | 0.596 | 1.714 | 0.468 |
| TTA | L | 0.48 | 1.632 | 0.019 | 1.556 | |
| TTG | L | 0.78 | 1.065 | 0.000 | 1.111 | |
| TCT | S | 1.14 | 1.961 | 0.525 | 2.242 | |
| ... | ... | ... | ... | ... | ... | ... |
| GAG | E | ★ | 1.16 | 0.557 | 1.611 | 0.583 |
| GGT | G | 0.64 | 2.340 | 0.049 | 2.293 | |
| GGC | G | ★ | 1.36 | 0.715 | 3.398 | 0.732 |
| GGA | G | 1.00 | 0.826 | 0.439 | 0.829 | |
| GGG | G | 1.00 | 0.118 | 0.114 | 0.146 |
61 rows × 6 columns
3.5 可视化（一）：RSCU 密码子偏好对比¶
(a) 以人源 RSCU 为横轴的散点图，观察病毒与疫苗序列偏离对角线的方向与程度；(b) 常见氨基酸（L / V / A / P / G / R）的密码子偏好柱状图。
In [13]:
# ================== 3.5 可视化（一）：RSCU 密码子偏好对比 ==================
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 人源 RSCU 为横轴，看病毒 / 疫苗偏离对角线的情况
ax = axes[0]
for col, color, marker in [("SARS-CoV-2", "#c0392b", "o"), ("疫苗 CDS", "#2471a3", "s")]:
ax.scatter(rscu_compare["Human（Kazusa）"], rscu_compare[col],
s=36, alpha=0.75, c=color, marker=marker,
edgecolor="white", linewidth=0.5, label=col)
ax.plot([0, 3.6], [0, 3.6], ls="--", lw=1, c="grey", label="y = x")
ax.set_xlabel("Human RSCU（Kazusa 参考）")
ax.set_ylabel("RSCU")
ax.set_title("(a) 61 个有义密码子的 RSCU 分布\n人源-病毒 r = %.2f，人源-疫苗 r = %.2f"
% (corr.loc["Human（Kazusa）", "SARS-CoV-2"],
corr.loc["Human（Kazusa）", "疫苗 CDS"]))
ax.legend(frameon=True, loc="upper left")
# (b) 挑几个氨基酸看具体密码子的使用差异
show_aa = ["L", "V", "A", "P", "G", "R"]
sub = rscu_compare[rscu_compare["氨基酸"].isin(show_aa)]
x = np.arange(len(sub))
width = 0.27
ax = axes[1]
for i, (col, color) in enumerate([("Human（Kazusa）", "#7f8c8d"),
("SARS-CoV-2", "#c0392b"),
("疫苗 CDS", "#2471a3")]):
ax.bar(x + (i - 1) * width, sub[col], width, label=col, color=color, alpha=0.9)
ax.set_xticks(x)
ax.set_xticklabels(sub.index, rotation=0, fontsize=8)
ax.set_xlabel("密码子")
ax.set_ylabel("RSCU")
ax.set_title("(b) 常见氨基酸（L/V/A/P/G/R）的密码子偏好")
ax.legend(frameon=True, fontsize=8)
ax.set_xlim(-0.8, len(sub) - 0.2)
plt.tight_layout()
plt.savefig("RSCU_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
3.6 可视化（二）：CAI / ENC–GC3s / GC 分布¶
(a) 各分析单元的 CAI 条形图；(b) ENC–GC3s 图（含 Wright 期望曲线），评价密码子偏倚强度；(c) 密码子三个位点的 GC 含量，突出 GC3s 的差异。
In [15]:
# ================== 3.6 可视化（二）：CAI / ENC-GC3s / GC 分布 ==================
def unit_color(label):
"""灰色 = SARS-CoV-2 病毒，红色 = 未优化原生 S 基因，蓝色 = 优化后疫苗 CDS"""
if label.startswith("SARS-CoV-2"):
return "#7f8c8d"
if "原生态" in label:
return "#c0392b"
return "#2471a3"
fig, axes = plt.subplots(1, 3, figsize=(17, 5))
# (a) CAI
ax = axes[0]
vals = summary["CAI(人源)"]
ax.barh(range(len(vals)), vals, color=[unit_color(l) for l in vals.index], alpha=0.9)
ax.set_yticks(range(len(vals)))
ax.set_yticklabels(vals.index, fontsize=8)
ax.invert_yaxis()
ax.set_xlim(0, 1.15)
ax.set_xlabel("CAI（以人源密码子表为参考）")
ax.set_title("(a) 密码子适应指数 CAI")
for i, v in enumerate(vals):
ax.text(v + 0.01, i, f"{v:.3f}", va="center", fontsize=8)
# (b) ENC–GC3s 图（Wright 期望曲线：完全无偏倚的基因会落在虚线上）
ax = axes[1]
s = np.linspace(0.001, 0.999, 300)
ax.plot(s * 100, 2 + s + 29 / (s ** 2 + (1 - s) ** 2), ls="--", lw=1, c="grey", label="Wright 期望曲线")
for (gc3s, e), grp in summary.groupby(["GC3s%", "ENC"]): # 坐标相同者合并标注
ax.scatter(gc3s, e, s=70, color=unit_color(grp.index[0]), edgecolor="white", zorder=3)
ax.annotate(" / ".join(grp.index), (gc3s, e), fontsize=6.5,
xytext=(5, 3), textcoords="offset points")
ax.set_xlabel("GC3s (%)")
ax.set_ylabel("ENC")
ax.set_title("(b) ENC–GC3s 图（越靠下偏倚越强）")
ax.legend(fontsize=8)
# (c) 各密码子位置 GC 含量
ax = axes[2]
x = np.arange(len(summary))
for i, m in enumerate(["GC1%", "GC2%", "GC3s%"]):
ax.bar(x + (i - 1) * 0.27, summary[m], 0.27, label=m)
ax.set_xticks(x)
ax.set_xticklabels(summary.index, rotation=40, ha="right", fontsize=7)
ax.set_ylabel("GC (%)")
ax.set_title("(c) 密码子各位点 GC 含量")
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig("codon_usage_metrics.png", dpi=150, bbox_inches="tight")
plt.show()
4. 如何通过S蛋白氨基酸序列设计mRNA疫苗，请搜索合适的软件，生成设计的结果¶
4.0 设计思路与工具选择¶
从蛋白序列设计一条 mRNA 疫苗构件，通常包含以下步骤：
- 反向翻译 + 密码子优化：把氨基酸序列反译为 DNA/mRNA 编码序列（CDS），并在保证翻译产物不变的前提下， 把密码子替换为宿主（人源）偏好的同义密码子，同时规避不利于合成、体内稳定性或表达的序列特征 （酶切位点、强二级结构、长重复序列、局部 GC 含量异常等）。
- 拼装非编码元件：加上 5'UTR（含 Kozak 序列，促进核糖体识别起始密码子）、3'UTR（提高 mRNA 稳定性）、 poly(A) 尾（抵抗降解、促进翻译），5' 端在体外转录后还会加帽（Cap，不体现在序列中）。
- 评价设计结果：用第 3 节已经实现的 CAI / RSCU / GC / GC3 / ENC 等指标，检验优化后的序列是否 确实向人源密码子偏好靠拢，并与已上市疫苗（BNT162b2、mRNA-1273）的实测序列做对比。
工具选择：密码子优化使用开源库 DNAChisel
（Edinburgh Genome Foundry 出品，与第 3 节使用的 python-codon-tables 同一生态）。它把序列设计建模为
「在满足一组约束（constraints）的前提下，最大化目标函数（objectives）」的优化问题，
专门用于 DNA/mRNA 序列的合成生物学设计，是目前该场景下最常用的开源方案之一。
UTR 部分使用人 β-珠蛋白（HBB，NCBI RefSeq NM_000518.5）的真实 5'/3'UTR ——
这是学术界报道中常用于增强外源 mRNA 稳定性的经典元件，直接从 NCBI 下载获得（真实序列，非虚构）。
4.1 密码子优化：用 DNAChisel 生成优化后的 Spike CDS¶
In [16]:
# ================== 4.1 使用 DNAChisel 对 Spike 蛋白进行密码子优化 ==================
# pip install dnachisel # 若环境中未安装，先执行本行（去掉注释）
from Bio import SeqIO
from Bio.Seq import Seq
from dnachisel import (
DnaOptimizationProblem, reverse_translate,
EnforceGCContent, AvoidHairpins, AvoidPattern, UniquifyAllKmers,
CodonOptimize, EnforceTranslation,
)
# 设计起点：第 1 节下载的 NCBI 注释 Spike 蛋白 YP_009724390.1
spike_record = SeqIO.read("SARS-CoV-2_Spike_protein.fasta", "fasta")
spike_protein = str(spike_record.seq).rstrip("*")
print(f"设计起点：{spike_record.id}，{len(spike_protein)} aa")
# 1) 反向翻译得到一条初始 CDS（只保证翻译正确，密码子尚未优化，常包含稀有密码子/重复序列）
initial_cds = reverse_translate(spike_protein, table="Standard")
# 2) 用 DNAChisel 在约束下做密码子优化
problem = DnaOptimizationProblem(
sequence=initial_cds,
constraints=[
EnforceTranslation(), # 翻译结果必须与原蛋白完全一致
AvoidPattern("BsaI_site"), # 避免常用金门克隆酶切位点
AvoidPattern("BsmBI_site"),
AvoidPattern("AATAAA"), # 避免 CDS 内部出现隐藏的 polyA 信号
AvoidHairpins(stem_size=20, hairpin_window=200), # 避免强二级结构（发卡）
UniquifyAllKmers(k=15), # 避免长片段重复，降低同源重组/测序风险
EnforceGCContent(mini=0.4, maxi=0.65, window=100), # 局部 GC 含量控制在易合成范围内
],
objectives=[CodonOptimize(species="h_sapiens", method="use_best_codon")], # 向人源密码子偏好优化
logger=None,
)
problem.resolve_constraints()
problem.optimize()
optimized_cds = problem.sequence
assert str(Seq(optimized_cds).translate()) == spike_protein, "优化后翻译结果与原蛋白不一致！"
print(f"\n初始 CDS 长度：{len(initial_cds)} nt")
print(f"优化后 CDS 长度：{len(optimized_cds)} nt（长度不变，仅同义密码子替换）")
print(f"密码子替换比例：{sum(a != b for a, b in zip(initial_cds, optimized_cds)) / len(initial_cds):.1%}")
print("\n约束校验：")
print(problem.constraints_text_summary())
设计起点：YP_009724390.1，1273 aa
初始 CDS 长度：3819 nt
优化后 CDS 长度：3819 nt（长度不变，仅同义密码子替换）
密码子替换比例：40.2%
约束校验：
===> SUCCESS - all constraints evaluations pass
✔PASS ┍ EnforceTranslation[0-3819]
│ Enforced by nucleotides restrictions
✔PASS ┍ AvoidPattern[0-3819](pattern:BsaI(GGTCTC))
│ Passed. Pattern not found !
✔PASS ┍ AvoidPattern[0-3819](pattern:BsmBI(CGTCTC))
│ Passed. Pattern not found !
✔PASS ┍ AvoidPattern[0-3819](pattern:AATAAA)
│ Passed. Pattern not found !
✔PASS ┍ AvoidHairpins[0-3819](stem_size:20, hairpin_window:200)
│ Score: 0. Locations: []
✔PASS ┍ UniquifyAllKmers[0-3819(+)](k:15)
│ Passed: no nonunique 15-mer found.
✔PASS ┍ EnforceGCContent[0-3819](mini:0.40, maxi:0.65, window:100)
│ Passed !
4.2 获取真实 UTR 元件并拼装完整 mRNA 疫苗序列¶
In [17]:
# ================== 4.2 拼装完整 mRNA 序列 ==================
# 5'UTR / 3'UTR 取自人 β-珠蛋白 (HBB) mRNA（NCBI RefSeq NM_000518.5），
# 这是学术界报道中常用于增强外源 mRNA 稳定性的真实元件（非虚构占位序列）。
import requests
HBB_FASTA = "HBB_NM_000518.5.fasta"
if not Path(HBB_FASTA).exists():
url = ("https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
"?db=nuccore&id=NM_000518.5&rettype=fasta&retmode=text")
resp = requests.get(url, timeout=60) # 直连 NCBI，不使用代理
resp.raise_for_status()
Path(HBB_FASTA).write_text(resp.text, encoding="utf-8")
print(f"已下载：{HBB_FASTA}")
hbb = str(SeqIO.read(HBB_FASTA, "fasta").seq)
# HBB mRNA 注释：CDS = 51..494（1-based），故 5'UTR = 1..50，3'UTR = 495..628
HBB_UTR5 = hbb[0:50]
HBB_UTR3 = hbb[494:]
POLYA = "A" * 100 # poly(A) 尾，简化为 100 nt
print(f"5'UTR（HBB, {len(HBB_UTR5)} nt）：{HBB_UTR5}")
print(f"3'UTR（HBB, {len(HBB_UTR3)} nt）：{HBB_UTR3}")
# 拼装：5'UTR(含天然 Kozak 上下文) + 优化后的 Spike CDS + 3'UTR + poly(A)
mrna_designed = HBB_UTR5 + optimized_cds + HBB_UTR3 + POLYA
designed_record = SeqIO.SeqRecord(
Seq(mrna_designed),
id="Designed_SARS-CoV-2_Spike_mRNA",
description=(f"5'UTR(HBB)+codon-optimized Spike CDS+3'UTR(HBB)+polyA100; "
f"total {len(mrna_designed)} nt"),
)
SeqIO.write(designed_record, "SARS-CoV-2_Spike_mRNA_designed.fasta", "fasta")
print(f"\n设计完成，总长度：{len(mrna_designed)} nt")
print(f" 5'UTR {len(HBB_UTR5)} + CDS {len(optimized_cds)} + 3'UTR {len(HBB_UTR3)} + polyA {len(POLYA)}")
print("已保存：SARS-CoV-2_Spike_mRNA_designed.fasta")
5'UTR（HBB, 50 nt）：ACATTTGCTTCTGACACAACTGTGTTCACTAGCAACCTCAAACAGACACC 3'UTR（HBB, 134 nt）：GCTCGCTTTCTTGCTGTCCAATTTCTATTAAAGGTTCCTTTGTTCCCTAAGTCCAACTACTAAACTGGGGGATATTATGAAGGGCCTTGAGCATCTGGATTCTGCCTAATAAAAAACATTTATTTTCATTGCAA 设计完成，总长度：4103 nt 5'UTR 50 + CDS 3819 + 3'UTR 134 + polyA 100 已保存：SARS-CoV-2_Spike_mRNA_designed.fasta
4.3 评价设计结果：CAI / GC / GC3 / ENC，并与已上市疫苗对比¶
复用第 3 节定义的 cai、gc_metrics、enc、optimal_codon_fraction 等函数，
比较「优化前」「优化后」与已获取的 BNT162b2 / mRNA-1273 / 原生态 S 基因的密码子使用特征。
In [18]:
# ================== 4.3 设计结果评价 ==================
design_units = {
"本设计-优化前(reverse_translate)": initial_cds,
"本设计-优化后(DNAChisel)": optimized_cds,
"S 基因（原生态，未优化）": vaccine_cds[NATIVE_S_NAME],
"BNT162b2（WHO，已上市）": vaccine_cds["WHO_BNT162b2"],
"mRNA-1273（Moderna，已上市）": vaccine_cds["ModernaMrna1273"],
}
design_rows = []
for label, seq in design_units.items():
counts = count_codons(seq)
design_rows.append({
"分析单元": label,
"CDS 长度(nt)": len(seq),
"CAI(人源)": cai(seq, human_weights),
"人源最优密码子%": optimal_codon_fraction(seq),
"ENC": enc(counts),
**gc_metrics(seq),
})
design_summary = pd.DataFrame(design_rows).set_index("分析单元").round(3)
design_summary.to_csv("vaccine_design_comparison.csv", encoding="utf-8-sig")
print("已保存：vaccine_design_comparison.csv\n")
design_summary
已保存：vaccine_design_comparison.csv
Out[18]:
| CDS 长度(nt) | CAI(人源) | 人源最优密码子% | ENC | GC% | GC1% | GC2% | GC3% | GC3s% | |
|---|---|---|---|---|---|---|---|---|---|
| 分析单元 | |||||||||
| 本设计-优化前(reverse_translate) | 3819 | 0.682 | 0.000 | 20.000 | 28.332 | 42.969 | 39.984 | 2.042 | 0.000 |
| 本设计-优化后(DNAChisel) | 3819 | 0.981 | 91.259 | 21.887 | 60.958 | 47.997 | 39.984 | 94.894 | 94.787 |
| S 基因（原生态，未优化） | 3822 | 0.730 | 19.727 | 44.159 | 37.310 | 45.290 | 39.953 | 26.688 | 25.180 |
| BNT162b2（WHO，已上市） | 3822 | 0.952 | 77.787 | 28.999 | 56.986 | 49.765 | 40.188 | 81.005 | 80.674 |
| mRNA-1273（Moderna，已上市） | 3822 | 0.981 | 91.500 | 22.037 | 62.297 | 51.256 | 40.188 | 95.447 | 95.429 |
4.4 可视化：优化前后的密码子偏好变化¶
In [19]:
# ================== 4.4 可视化：优化前后对比 ==================
design_rscu = {
label: rscu_from_counts(count_codons(seq))
for label, seq in design_units.items()
}
design_rscu_df = pd.DataFrame(design_rscu).reindex(SENSE_CODONS)
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) CAI 优化前 vs 优化后 vs 已上市疫苗
ax = axes[0]
vals = design_summary["CAI(人源)"]
colors = ["#c0392b", "#27ae60", "#7f8c8d", "#2471a3", "#2471a3"]
ax.barh(range(len(vals)), vals, color=colors, alpha=0.9)
ax.set_yticks(range(len(vals)))
ax.set_yticklabels(vals.index, fontsize=8)
ax.invert_yaxis()
ax.set_xlim(0, 1.15)
ax.set_xlabel("CAI（以人源密码子表为参考）")
ax.set_title("(a) 密码子优化前后的 CAI 变化")
for i, v in enumerate(vals):
ax.text(v + 0.01, i, f"{v:.3f}", va="center", fontsize=8)
# (b) 优化前 vs 优化后的 RSCU 散点（以人源 RSCU 为横轴）
ax = axes[1]
ax.scatter(human_rscu_series := pd.Series(human_rscu).reindex(SENSE_CODONS),
design_rscu_df["本设计-优化前(reverse_translate)"],
s=32, alpha=0.7, c="#c0392b", label="优化前", edgecolor="white", linewidth=0.4)
ax.scatter(human_rscu_series,
design_rscu_df["本设计-优化后(DNAChisel)"],
s=32, alpha=0.7, c="#27ae60", label="优化后", edgecolor="white", linewidth=0.4)
ax.plot([0, 3.6], [0, 3.6], ls="--", lw=1, c="grey", label="y = x（与人源完全一致）")
ax.set_xlabel("Human RSCU（Kazusa 参考）")
ax.set_ylabel("RSCU")
ax.set_title("(b) 密码子优化前后向人源偏好靠拢的效果")
ax.legend(frameon=True, fontsize=8, loc="upper left")
plt.tight_layout()
plt.savefig("vaccine_design_optimization.png", dpi=150, bbox_inches="tight")
plt.show()
如何解读 (b) 图：CodonOptimize(method="use_best_codon") 对每种氨基酸都固定选用人源使用频率最高的那一个密码子，
因此优化后的 RSCU 只会取 0（未使用的同义密码子）或该氨基酸简并度对应的最大值（如 6 重简并的密码子家族会冲到 6），
而不是复现人源天然的 RSCU 分布形状——这与 mRNA-1273 / BNT162b2 等真实疫苗序列的密码子使用模式（第 3 节，见
rscu_compare 表）并不完全相同，但两者的 CAI 同样接近 1，说明「贴近人源偏好」这一优化目标本身已经达成。
若想更贴近天然 RSCU 分布，可将 method 改为 "match_codon_usage"（按人源频率随机采样而非总取最优）。