材料与试剂

SARS-CoV-2基因组序列于2022年5月1日从GISAID数据库[13]获得。从人类以外的宿主收集和/或长度小于2,900个核苷酸且未知核苷酸超过0.05%的基因组被过滤掉了。超过 500 万个基因组被保留用于进一步分析。

设备

1.生物信息平台（CentOS，7.2.1511）
2.Windows 台式电脑（v11，2021 年 6 月 6 日）

软件和数据集

1.Perl（v5.34.0，2021 年 5 月 20 日）
2.Python（v3.9.0，2020 年 5 月 10 日）
3.Gcc（v8.4.1，2020 年 9 月 28 日）
4.SeqKit（v2.3.0，2022 年 12 月 8 日）
5.小地图2（v2.24，2021年12月26日）
6.ViralMSA（v3，2007 年 6 月 29 日）
7.CD-hit（v4.8.1，2018 年 1 月 9 日）
8.更改（v1.3.4，2016 年 10 月 30 日）
9.野兽（v1.10.4，2018 年 9 月 11 日）
10.Jdk（v17_windows-x64_bin，2021 年 6 月 16 日）
11.Jre（8u341-windows-x64，2021 年 5 月 6 日）
12.追踪器（v1.7.2，2018 年 5 月 1 日）
13.全球的共享所有流感数据倡议 (GISAID) ( https://gisaid.org/ )（访问日期，2022 年 1 月 5 日）
所有个性化脚本均已存放在 GitHub 中： https://github.com/XiaowenH/SARS-CoV-2-Classification （访问日期，2023 年 5 月 9 日）。图 1 是 GitHub 中文件的屏幕截图。

图1。存放在 GitHub 中的个性化脚本的屏幕截图
程序

A.设置工作环境

1.创建conda的Python环境。

$ wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ bash Miniconda3-latest-Linux-x86_64.sh -b -p ~/miniconda
$ 导出 PATH="/$HOME/miniconda/ bin:$ PATH"
$ conda create -n myenv python=3.9
$ conda 激活 myenv
#请注意“$HOME”对应的是你要在conda环境中安装的PATH（下同）。

2.安装以下软件：
a.在 Linux 平台上安装 Seqkit。

$ conda install -c bioconda seqkit

b.在Linux平台上安装minimap2。

$ git 克隆 https://github.com/lh3/minimap2
$ cd minimap2 && make

c.在 Linux 平台中安装 ViralMSA [14]。

$ wget“https://raw.githubusercontent.com/niemasd/ViralMSA/master/ViralMSA.py”
$ chmod a+x ViralMSA.py
$ sudo mv ViralMSA.py /usr/local/bin/ViralMSA.py

d.在 Linux 平台上安装 CD-HIT。

$ git 克隆 https://github.com/weizhongli/cdhit.git
$ cd cd-命中
$ 制作
$ cd cd-hit-auxtools
$ 制作

e.在 Linux 平台上安装 Alter。

$ git 克隆 https://github.com/sing-group/ALTER.git
$ cd 更改
$ /$HOME/apache-maven-3.8.4/bin/mvn _

f.在 Windows 平台中安装 BEAST [15]。

# 下载 BEAST.v1.10.4 在
http://beast.community/programs

g.上安装 BEAGLE 。

# 下载 BEAGLE 在
https://github.com/beagle-dev/beagle-lib

h.在 Windows 平台中安装 Tracer v1.7.2。

# 下载 BEAGLE 在
http://beast.community/tracer

B.下载病毒基因组和元数据

1.登录后从 GISAID 数据库 [16] 下载 SARS-CoV-2 基因组和 meta.tsv 文件 ( https://gisaid.org/ )。请注意，meta.tsv 中每列的内容可能会发生变化（例如，入藏号被放入 2022 年 5 月 1 日下载的元数据中的第 3 列，但在 2023 年 5 月 1 日下载的元数据中放入第 5 列）。
2.解压文件。

$ xz -dsequences_fasta_2022_04_29.tar.xz
$ tar -xvfsequences_fasta_2022_04_29.tar
$ xz -dmetadata_tsv_2022_04_29.tar.xz
$ tar -xvfmetadata_tsv_2022_04_29.tar

C.使用自定义脚本按登录号（例如 EPI_ISL_4405694）重命名基因组 Fasta 文件中的分类单元

$ perl rename_fasta_taxa_to_tsv_acc_column3.pl -t meta.tsv -f序列.fasta

D.获取 GISAID 参考序列

1.创建一个名为“accession.txt”的文件，粘贴到 EPI_ISL_402124 中，然后保存。

$ echo 'EPI_ISL_402124' >accession.txt
# EPI_ISL_402124 是GISAID 数据库和许多研究中使用的参考序列。

2.从总序列文件中获取参考序列。

$ cat序列.fasta | / $HOME/ seqkit grep -f accession.txt -t dna -j 10 -o reference_wiv04.fasta

E.检索具有完整、高覆盖率序列和准确日期的基因组，并从人类宿主中取样

基因组的登录号、宿主、完整性和覆盖度分别位于2022年4月29日元数据中的第3、8、18和19列。样本采集日期位于第4列。列号不同日期下载的元数据可能会有所不同。
1.过滤metadata.tsv以获取具有完整基因组、高覆盖度和来自人类宿主的种质。

$ awk -F '\t' '{ if( $18 == "True" && $19 == "True" && $8 == "Human") print}'metadata.tsv >global_2022_04_29_ human_complete_metadata.tsv

2.使用自定义脚本：dates_filter.pl，筛选具有准确样本采集日期的metadata.tsv。

$ perldates_filter.pl -t global_2022_04_29_ human_complete_metadata.tsv -o global_2022_04_29_ human_complete_dates_metadata.tsv

3.将选定的入藏号从元数据打印到文件中。

$ awk -F '\t' '{print $3}' global_2022_04_29_ human_complete_dates_metadata.tsv > global_2022_04_29_ human_complete_dates_accessions.txt

4.检索完整、高覆盖度和准确日期的基因组序列。

$ cat序列.fasta | / $ HOME/seqkit grep -f global_2022_04_29_ human_complete_dates_accessions.txt -t dna -j 10 -o global_2022_04_29_ human_complete_dates_genome.fasta

F.比对基因组序列

$ /$HOME/ViralMSA-master/ViralMSA.py -e [您的电子邮件] -s global_2022_04_29_ human_complete_dates_genome.fasta -o 输出 -r reference_wiv04.fasta --omit_ref

G.所有基因组分类

1.使用自定义脚本 fetch_核苷酸_from_alignments.pl 检索六个高度链接的位点（图 2）。

$ perl fetch_核苷酸_from_alignments.pl -f global_2022_04_29_ human_complete_dates_genome.fasta.aln
-r 241-241,3037-3037,8782-8782,14408-14408,23403-23403,28144-28144 -o global_2022_04_29_ human_complete_dates_genome.fasta_six_sites.txt

图 2. 基因组中六个连锁位点的输出格式。从 SARS-CoV-2 基因组中检索位点 241、3037、8782、14408、23403 和 28144 的六个核苷酸，与参考基因组 wiv04 EPI_ISL_402124 进行比对。

2.通过六个连锁突变位点将基因组分类为单倍型（图 3）。

$ cat global_2022_04_29_ human_complete_dates_genome.fasta_six_sites.txt | /$HOME/seqkit grep -s -i –p CCTCAC > global_2022_04_29_ human_complete_dates_DS.txt
$ cat global_2022_04_29_ human_complete_dates_genome.fasta_six_sites.txt | /$HOME/seqkit grep -s -i -p CCCCAT > global_2022_04_29_ human_complete_dates_DL.txt
$ cat global_2022_04_29_ human_complete_dates_genome.fasta_six_sites.txt | /$HOME/seqkit grep -s -i -p TTCTGT > global_2022_04_29_ human_complete_dates_GL.txt
$ /$HOME /seqkit stats global_2022_04_29_ human_complete_dates_DS.txt >haplotype.statistics.txt
$ /$HOME /seqkit stats global_2022_04_29_ human_complete_dates_DL.txt >>haplotype.statistics.txt
$ /$HOME / seqkit stats global_2022_04_29_ human_complete_dates_GL.txt >>haplotype.statistics.txt

图 3. 全球 SARS-CoV-2 基因组中单倍型的饼图。 DS (CCTCAC)、DL (CCCCAT) 和 GL (TTCTGT) 是按六个连锁位点分类的三种主要单倍型。

3.获取每个单倍型的登录号。

$ /$HOME/seqkit seq global_2022_04_29_ human_complete_dates_DS.txt -n > global_DS.accessions.txt
$ /$HOME/seqkit seq global_2022_04_29_ human_complete_dates_DL.txt -n > global_DL.accessions.txt
$ /$HOME/seqkit seq global_2022_04_29_ human_complete_dates_GL.txt -n > global_GL.accessions.txt

4.检索每个单倍型的比对基因组序列。

$ cat global_2022_04_29_ human_complete_dates_genome.fasta.aln | /$HOME/seqkit grep -f global_DS.accessions.txt -t dna -j 10 -o global_DS.aln
$ cat global_2022_04_29_ human_complete_dates_genome.fasta.aln | /$HOME/seqkit grep -f global_DL.accessions.txt -t dna -j 10 -o global_DL.aln
$ cat global_2022_04_29_ human_complete_dates_genome.fasta.aln | /$HOME/seqkit grep -f global_GL.accessions.txt -t dna -j 10 -o global_GL.aln

H.疫情早期收集的早期基因组分类（2020年4月上旬至下旬）

样本采集日期位于2022 年 4 月 29 日元数据的第 4 列中。
1.检索早期基因组的登录号。

$ awk -F '\t' '{ if( $4~/2019-12/) print}' global_2022_04_29_ human_complete_dates_metadata.tsv > global_early_meta.tsv
$ awk -F '\t' '{ if( $4~/2020-01/) print}' global_2022_04_29_ human_complete_dates_metadata.tsv >> global_early_meta.tsv
$ awk -F '\t' '{ if( $4~/2020-02/) print}' global_2022_04_29_ human_complete_dates_metadata.tsv >> global_early_meta.tsv
$ awk -F '\t' '{ if( $4~/2020-03/) print}' global_2022_04_29_ human_complete_dates_metadata.tsv >> global_early_meta.tsv
$ awk -F '\t' '{ if( $4~/2020-04/) print}' global_2022_04_29_ human_complete_dates_metadata.tsv >> global_early_meta.tsv
$ awk -F '\t' '{print $3}' global_early_meta.tsv >global_early_accessions.txt

2.检索早期基因组的比对序列。

$ cat global_2022_04_29_ human_complete_dates_genome.fasta.aln | /$HOME/seqkit grep -f global_early_accessions.txt -t dna -j 10 -o global_early_all_haplotype.fasta.aln

3.使用自定义脚本 fetch_核苷酸_from_alignments.pl 检索六个高度链接的位点。

$ perl fetch_核苷酸_from_alignments.pl -f global_early_all_haplotype.fasta.aln
-r 241-241,3037-3037,14408-14408,23403-23403,28144-28144 -o global_early.aln_six_sites.txt

4.对基因组进行分类。

$ cat global_early.aln_six_sites.txt | /$HOME/seqkit grep -s -i -p CCTCAC > global_early_DS.txt
$ cat global_early.aln_six_sites.txt | /$HOME/seqkit grep -s -i -p CCCCAT > global_early_DL.txt
$ cat global_early.aln_six_sites.txt | /$HOME/seqkit grep -s -i -p TTCTGT > global_early_GL.txt

5.获取每个单倍型的登录号。

$ /$HOME/seqkit seq global_early_DS.txt -n > global_early_DS.accessions.txt
$ /$HOME/seqkit seq global_early_DL.txt -n > global_early_DL.accessions.txt
$ /$HOME/seqkit seq global_early_GL.txt -n > global_early_GL.accessions.txt

6.检索每个单倍型的比对基因组序列。

$ cat global_early_all_haplotype.fasta.aln | /$HOME/seqkit grep -f global_ Early_DS.accessions.txt -t dna -j 10 -o global_ 早期_DS.aln
$ cat global_early_all_haplotype.fasta.aln | /$HOME/seqkit grep -f global_ Early_DL.accessions.txt -t dna -j 10 -o global_ Early_DL.aln
$ cat global_early_all_haplotype.fasta.aln | /$HOME/seqkit grep -f global_ Early_GL.accessions.txt -t dna -j 10 -o global_ 早期_GL.aln

I.以三种单倍型的早期基因组为例进行贝叶斯系统动力学分析

1.过滤掉未知数高于0.05%的基因组。
a.过滤掉未知核苷酸高于 0.05% 的 DS (CCTCAC) 单倍型基因组。

# 计算每个单倍型中每个序列的核苷酸组成
$ perl base_stats.pl global_early_DS.aln
# 选取未知数 < 0.05% 的种质
$ awk -F '\t' '{ if( $6 < 0.0005) print$1}' global_early_DS.aln.stats > global_early_DS.aln_0.0005N_ID.txt
# 获取未知数 <0.05% 的对齐基因组序列
$ cat global_early_DS.aln |/$HOME/seqkit grep -f global_early_DS.aln_0.0005N_ID.txt -t dna -j 10 -o global_early_DS_0.0005N.aln

b.过滤掉未知核苷酸含量高于 0.05% 的 DL (CCCCAT) 单倍型基因组。

# 计算每个单倍型中每个序列的核苷酸组成
$ perl base_stats.pl global_early_DL.aln
# 选取未知数 < 0.05% 的种质
$ awk -F '\t' '{ if( $6 < 0.0005) print$1}' global_early_DL.aln.stats > global_early_DL.aln_0.0005N_ID.txt
# 获取未知数 <0.05% 的对齐基因组序列
$ cat global_early_DL.aln |/$HOME/seqkit grep -f global_early_DL.aln_0.0005N_ID.txt -t dna -j 10 -o global_early_DL_0.0005N.aln

c.过滤掉未知核苷酸含量高于0.05%的GL（TTCTGT）单倍型的基因组。

# 计算每个单倍型中每个序列的核苷酸组成。
$ perl base_stats.pl global_early_GL.aln
# 选取未知数 < 0.05% 的种质
$ awk -F '\t' '{ if( $6 < 0.0005) print$1}' global_early_GL.aln.stats > global_early_GL.aln_0.0005N_ID.txt
# 获取未知数 <0.05% 的对齐基因组序列
$ cat global_early_GL.aln | /$HOME/seqkit grep -f global_early_GL.aln_0.0005N_ID.txt -t dna -j 10 -o global_early_GL_0.0005N.aln

d.以 0.9997 的阈值去除冗余基因组。

$ /$HOME/cd-hit-v4.8.1-2019-0228/cd-hit-est -i global_early_DS_0.0005N.aln -o global_early_DS_0.0005N_CDHit0.9997.fas -M 2500 -c 0.9997 -aL 1 -aS 1 -d 0
$ /$HOME/cd-hit-v4.8.1-2019-0228/cd-hit-est -i global_early_DL_0.0005N.aln -o global_early_DL_0.0005N_CDHit0.9997.fas -M 2500 -c 0.9997 -aL 1 -aS 1 -d 0
$ /$HOME/cd-hit-v4.8.1-2019-0228/cd-hit-est -i global_early_GL_0.0005N.aln -o global_early_GL_0.0005N_CDHit0.9997.fas -M 2500 -c 0.9997 -aL 1 -aS 1 -d 0

e.将格式更改为 NEX。

$ java -jar /$HOME/ALTER/alter-lib/target/ALTER-1.3.4-jar-with-dependency.jar -i global_early_DS_0.0005N_CDHit0.9997.fas - ia - o global_early_DS_0.0005N_CDHit0.9997.fas。 nex -of NEXUS -oo Windows -op MrBayes
$ java -jar /$HOME/ALTER/alter-lib/target/ALTER-1.3.4-jar-with-dependency.jar -i global_early_DL_0.0005N_CDHit0.9997.fas - ia - o global_early_DL_0.0005N_CDHit0.9997.fas。 nex -of NEXUS -oo Windows -op MrBayes
$ java -jar /$HOME/ALTER/alter-lib/target/ALTER-1.3.4-jar-with-dependency.jar -i global_early_GL_0.0005N_CDHit0.9997.fas - ia - o global_early_GL_0.0005N_CDHit0.9997.fas。 nex -of NEXUS -oo Windows -op MrBayes

f.从元数据中获取加入日期。

$ awk -F "\t" -v OFS="\t" '{print $ 3,$ 4}' global_early_meta.tsv > global_early_dates.txt

g.使用自定义脚本 date_convertor.pl 将日期转换为十进制日期。

$ perl date_convertor.pl global_early_dates.txt > global_early_decimal_dates

h.在Windows平台上创建BEAST XML文件。

#双击图标打开 BEAUti- v1.10.4
# 加载NEX格式的基因组比对文件
# 在文件global_early_decimal_dates中导入日期（样本采集日期）
# 设置分析选项，例如将“Clocks”设置为严格或宽松，将“Tree Prior”设置为“Bayesian Skyline”或“Bayesian SkyGrid”，将MCMC“Length of Chain”设置为12亿代。当所有参数的解释平方和 (ESS) 显着 (>200) 时，可以停止运行。参数全部设置完毕后，单击“生成 BEAST 文件”。

i.运行野兽。

# 双击图标打开BEAST v1.10.4
# 加载 BEAST XML 文件，单击“运行”。运行可能需要几天时间，直到所有 ESS 都显着。

j.通过 Tracer 可视化 tMRCA（图 4）。

# 双击图标打开 Tracer v1.7.2
# 导入BEAST创建的日志文件。 tMRCA 的后部分布可以通过单击“年龄（根）”和“边缘密度”来显示。 “估计”中提供了平均 tMRCA 和 95% HPD 间隔。

图 4.使用 Tracer 总结的早期 SARS-CoV-2 基因组的三种主要单倍型的最近共同祖先时间 ( tMRCA) 估计的屏幕截图。日期以十进制显示。

协议验证

该协议或其部分内容已在以下研究文章中使用和验证：
•胡等人。 [17]。基于 SARS-CoV-2 Spike-614 和 NS8-84 位点的基因组表征揭示了 COVID-19 大流行的两个主要可能爆发点。《公共科学图书馆 ONE》 18(6)：e0279221。 https://doi.org/10.1371/journal.pone.0279221 （图 4，面板 1；图 6）。
•关等人。 [12]。基因组分析显示，SARS-CoV-2 主要单倍型的分离和平行进化比其在中国的出现要早得多。 Science in One Health 2(2023)， https://doi.org/10.1016/j.soh.2023.100041 （图 1；图 3；图 4）。

致谢

我们衷心感谢所有数据贡献者，包括作者及其负责获取样本的原始实验室，以及负责生成基因序列和元数据并通过本研究所基于的 GISAID 计划进行共享的提交实验室。该研究得到国家重点研发计划、中央公益性科研机构基础研究基金项目JZ（1630052020022）和四川省科技厅项目LY（2019JDJQ0035）的资助。该协议在 PLoS ONE 18(6): e0279221 中有部分描述。土井： 10.1371/journal.pone.0279221 （Hu 等人[17]）和 Science in One Health (2023)，Doi: 10.1016/j.soh.2023.100041 ）（ Guan等人[12]）。

利益争夺

作者声明没有竞争利益。

道德考虑

该协议不涉及任何实验。

参考

1.童 Y.、刘 W.、刘 P.、刘 WJ、王 Q. 和高 GF (2021)。病毒的起源：发现需要时间、国际资源和合作。柳叶刀398（10309）：1401–1402。
2.Nahmias, A.、Weiss, J.、Yao, X.、Lee, F.、Kodsi, R.、 Schanfield , M.、Matthews, T.、Bolognesi, D.、Durack, D.、 Motulsky , A.,等人。（1986）。 1959 年，中非人类感染 HTLV III/LAV 样病毒的证据。柳叶刀327（8492）：1279-1280。
3.Worobey, M.、Gemmel, M.、 Teuwen , DE、 Haselkorn , T.、Kunstman, K.、Bunce, M.、Muyembe, JJ、Kabongo, JM、 Kalengayi , RM、Van Marck, E. 等人。（2008）。金沙萨 HIV-1 广泛多样性的直接证据自然455（7213）：661-664。
4.Zaki, AM、van Boheemen , S.、 Bestebroer , TM、Osterhaus, AD 和 Fouchier, RA (2012)。从沙特阿拉伯一名肺炎患者身上分离出新型冠状病毒。 N. 英格兰。 J. Med。 367（19）：1814-1820。
5.Azhar, EI、El- Kafrawy , SA、Farraj, SA、Hassan, AM、Al-Saeed, MS、Hashem, AM 和 Madani, TA (2014)。 MERS 冠状病毒从骆驼向人类传播的证据。 N. 英格兰。 J. Med。 370（26）：2499-2505。
6.Müller, MA、Corman, VM、 Jores , J、Meyer, B、Younan, M、Liljander, A、Bosch、BJ、 Lattwein , E、Hilali, M、Musa, BE 等人。（2014）。东非骆驼中的 MERS 冠状病毒中和抗体，1983 年至 1997 年。新发传染病20(12):2093–2095。
7.朱宁、张东、王文、李X、杨本、宋J、赵X、黄本、施文、卢荣、等人。（2020）。 2019 年来自中国肺炎患者的新型冠状病毒。 N. 英格兰。 J. Med。 382（8）：727-733。
8.吴F.、赵S.、于B.、陈YM、王W.、宋ZG、胡Y.、陶ZW、田JH、裴YY等。（2020）。一种与中国人类呼吸道疾病相关的新型冠状病毒。自然579（7798）：265–269。
9.李Q.、管X.、吴P.、王X.、周L.、童Y.、任R.、梁KS、刘EH、黄JY等。（2020）。中国武汉新型冠状病毒感染肺炎的早期传播动态。 N. 英格兰。 J. Med。 382（13）：1199-1207。
10.Pekar, J.、Worobey, M.、Moshiri, N.、Scheffler, K. 和 Wertheim, JO (2021)。湖北省 SARS-CoV-2 指示病例的时间安排。科学372(6540)：412–417。
11.Cheng, C. 和Zhang, Z. (2023)。 SARS-CoV-2在世界范围内的分化早于中国大陆。科学。中国人寿科学． 66（6）：1440-1443。
12.管S.、胡X.、易G.、姚L.和张J. (2023)。 SARS-CoV-2 单倍型的基因组分析：主要单倍型的分离和平行进化发生的时间远远早于它们在中国的出现。科学与健康2：100041。
13.Khare, S.、Gurry, C.、Freitas, L.、Schultz, MB、Bach, G.、Diallo, A.、Akite, N.、Ho, J.、Lee, RT、Yeo, W. 等。（2021）。 GISAID 在大流行应对中的作用。中国疾控中心周刊3(49)：1049–1051。
14.莫希里，N.（2020）。 ViralMSA：大规模可扩展的参考指导的病毒基因组多序列比对。生物信息学（英国牛津） 37（5）：714-716。
15.Suchard , MA、 Lemey , P.、Baele, G.、Ayres, DL、Drummond, AJ 和 Rambaut, A. (2018)。使用 BEAST 1.10 进行贝叶斯系统发生和系统动力学数据集成。病毒进化。 4（1）：vey016。
16.Shu, Y. 和 McCauley, J. (2017)。 GISAID：共享所有流感数据的全球倡议——从愿景到现实。欧洲监视22(13)：e30494。
17. 胡X.，穆Y.，邓R.，易G.，姚L.和张J.(2023)。基于 SARS-CoV-2 Spike-614 和 NS8-84 位点的基因组表征揭示了 COVID-19 大流行的两个主要可能爆发点。《公共科学图书馆一号》 18(6)：e0279221。