本文檔采用 ARTIC Network 流程進行實驗，可用于illumina平臺和nanopore平臺對SARS-Cov-2病毒進行基因組測序。
@indexofire indexofire@gmail.com
1. cDNA獲取
1.1 試劑與耗材
- RNA提取試劑
- Invitrogen SuperScript? IV First-Strand Synthesis System
- 0.2mL PCR管
SSIV的酶無RNaseH活性，獲得長片段(～12kb)的全長cDNA效果更好。我們使用SSIII，有中等活性的RNaseH活性，獲得的cDNA進行測序也可以成功，在目標reads豐度上SSIV似乎更好一些。
1.2 實驗流程
- 采用RNA試劑盒手工提取樣本總RNA，也可以用機器自動提取。
RNA通過Qiagen RNeasy Mini Kit 手工提取試劑盒或者機器提取的都可以，PCR實驗結果CT無明顯差異。有文章比較過Qiagen不同的RNA提取試劑盒的宏基因組測序reads產出效果，可以供大家參考。
- 取EP管，依次裝入以下試劑，分裝到0.2mL PCR中，加入RNA模板，PCR管用移液器吹吸混勻，離心將管壁殘液甩下。
| 試劑 | 1x | 8x |
|---|---|---|
| 50μM random hexamers | 1 μl | 8μl |
| 10mM dNTPs mix (10mM each) | 1 μl | 8μl |
| Template RNA | 11 μl | 11μl * 8 |
| Total | 13 μl | 13μl * 8 |
注意事項：樣本RNA做過SARS-CoV-2 qPCR測試，可以根據(jù)CT值初步判斷RNA含量。如果RNA量在CT18～35之間，可以直接加入模板RNA。如果在12～15之間，100倍稀釋；如果在15～18之間，10倍稀釋。對于臨床病例來說，大部分樣本CT值在18～35之間。
使用random hexamers是為了最大程度獲得全長cDNA，SSIV自帶的random hexamers濃度是50ng/μl，6mers長度，如果加1μl，相當與摩爾數(shù) 25e-6μmol，因此根據(jù)表中量，需要加2μl。但根據(jù)試劑盒說明書中也是加1μl量，因此這里參考說明書。
推薦加入樣本前，在PCR儀上預熱。
- 反應體系在PCR儀器中 65°C 5min。然后立即置于冰上至少1min。
- 按下表配置反應體系，添加到冰上的PCR管中。置于PCR儀器中，42°C 50min，70°C 10min，結束后5°C 保溫。這一步要設置105°C熱蓋。
| 試劑 | 1x | 8x |
|---|---|---|
| 5x SSIV Buffer | 4 μl | 32 μl |
| 100mM DTT | 1 μl | 8 μl |
| RNaseOUT RNase Inhibitor | 1 μl | 8 μl |
| SSIV Reverse Transcriptase | 1 μl | 8 μl |
| 前一步PCR管內液體 | 13 μl | 13 μl * 8 |
| Total | 20 μl | 20 μl * 8 |
以上配液操作在PCR工作站或生物安全柜內進行。使用前紫外處理30分鐘并通風。
2. PCR擴增
2.1 試劑與耗材
ARTIC Network流程里使用的是NEB Q5酶，試用Qiagen的多重PCR產品，效果也可以。
2.2 實驗流程
- 將98組凍干引物 8000rpm 離心 10min。將每個引物配置成100uM濃度。
- 準備2個EP管，標記為P1和P2。將98組引物分為奇數(shù)和偶數(shù)組，每個組各98個100uM的引物，每管吸取5uL分別到Pool1和2中，最終每管含100uM濃度的引物共980uL。最終混合引物中每條引物濃度約為1uM。
- 震蕩混勻，離心甩干后，各取20uL P1/2 引物，加入180uL水。使得Pool1和Pool2引物終濃度為10uM，每條引物終濃度約為0.1uM。
以上配液操作在PCR工作站或生物安全柜內進行。使用前紫外處理30分鐘并通風。
- 每個樣本取2個PCR管，使用NEB Q5 Hot Start DNA Polymerase試劑或者Qiagen Multiplex PCR試劑進行PCR擴增，配液參見表1/2。
| NEB Q5 Hot Start DNA Polymerase | Pool1 | Pool2 |
|---|---|---|
| 5X Q5 Reaction Buffer | 5 μl | 5 μl |
| 10 mM dNTPs | 0.5 μl | 0.5 μl |
| Q5 Hot Start DNA Polymerase | 0.25 μl | 0.25 μl |
| Primer Pool 1 or 2 (10μM) | 3.6 μl | 3.6 μl |
| Nuclease-free water | 13.15 μl | 13.15 μl |
| cDNA | 2.5uL | 2.5uL |
| Total | 25 μl | 25 μl |
| Qiagen Multiplex PCR Kit | Pool1 | Pool2 |
|---|---|---|
| 2X Qiagen Multiplex PCR Master Mix | 12.5 μl | 12.5 μl |
| Primer Pool 1 or 2 (10μM) | 3.6 μl | 3.6 μl |
| 5X Q Solution | 2.5 μl | 2.5 μl |
| Nuclease-free water | 3.9 μl | 3.9 μl |
| cDNA | 2.5uL | 2.5uL |
| Total | 25 μl | 25 μl |
根據(jù)體系優(yōu)化，每個引物的終濃度為0.015uM。
- NEB Q5 的PCR程序為：
| 序號 | 溫度 | 時間 | 前往 |
|---|---|---|---|
| 1. | 98°C | 30s | |
| 2. | 98°C | 15s | |
| 3. | 65°C | 5min | 前往2, 循環(huán)25～35cycles |
| 4. | 4°C | Hold |
- Qiagen Multiplex Kit的PCR程序為：
| 序號 | 溫度 | 時間 | 前往 |
|---|---|---|---|
| 1. | 95°C | 15min | |
| 2. | 94°C | 30s | |
| 3. | 58°C | 90s | |
| 4. | 72°C | 45s | 前往2, 循環(huán)25～35cycles |
| 5. | 72°C | 10min | |
| 6. | 4°C | Hold |
Ct18-21的樣本用25循環(huán), Ct 35的樣本用35個循環(huán)。CT值介于21～35的根據(jù)分布自行設置。
應用本方法到其他病原檢測時，需要設計多重PCR的話，可以使用這里的引物設計工具獲得引物組。
3. Amplicon純化回收
3.1 試劑與耗材
3.2 實驗流程
- 將每個樣本的Pool1和Pool2混合，因為測序文庫與上樣量有關，本實驗流程設計的是針對8個以上樣本的上樣量。建議檢測樣本帶質控DNA做為陰性對照，以便檢查是否有污染序列。
- 將AMPure XP Beads 提前30min取出置于常溫，使用前在震蕩儀上混合均勻。
- 取50uL AMPure XP Beads 到 50uL 混合pooling Amplicon中。震蕩混勻（增加震蕩時間能極大提高回收率），離心甩干。
1x 磁珠吸附250～1000bp效率較好。具體參見
- 室溫放置5min。
- 將EP管放置于磁力架上2min，直至液體澄清。
- 吸棄液體，用200uL新鮮配置的70～80%乙醇洗2次。
- 用10uL移液器吸棄殘液，開蓋1min，讓乙醇揮發(fā)。
- 從磁力架上取下，用15～30uL Elution Buffer或者水重懸磁珠，用槍頭或手指輕輕混勻，靜置2min。
- 置于磁力架上，液體澄清后吸取30uL液體到一個新的EP管中。取1uL進行Qubit定量。
樣本濃度與PCR產物回收相關性：
- 10E5大約200ng/uL
- 10E4大約120ng/uL
- 10E3大約80ng/uL
- 10E2大約25ng/uL
- 10E1大約10ng/uL
4. 文庫構建
4.1 Nanopore 測序文庫
4.1.1 試劑與耗材
- Nanopore EXP-NBD104/114
- Nanopore SQK-LSK109
- Nanopore 測序芯片制備試劑盒 EXP-FLP002
- Eppendorf DNA LoBind Tubes
- NEBNext? Ultra? II DNA Library Prep Kit for Illumina E7645
4.1.2 實驗流程
注意此流程是根據(jù)6~24個樣本混樣的量規(guī)劃的。如果只做單個樣本，可以不加barcodes，但需要調整input DNA量，建議加入DNA達40ng，保證上機時能約有20ng DNA，達到最佳的芯片測序分子數(shù)。
- 將定量的DNA稀釋成1ng/uL，該濃度是針對700bp長度的Amplicon，如果片段長400bp，建議濃度為2ng/uL，這樣分子濃度約為100～200fmol。
- 取0.2mLPCR管，根據(jù)下表進行配液，進行末端修復。
| 試劑 | 量 |
|---|---|
| DNA amplicons | 5 μl |
| Nuclease-free water | 7.5 μl |
| Ultra II End Prep Reaction Buffer | 1.75 μl |
| Ultra II End Prep Enzyme Mix | 0.75 μl |
| 總量 | 15 μl |
- 室溫放置10min。置于PCR儀上，65°C 5min，立即置于冰上至少1min
- 在PCR管中加入以下試劑：
| 試劑 | 量 |
|---|---|
| 連接NBXX的混合液 | 15 μl |
| NBXX barcode | 2.5 μl |
| Ultra II Ligation Master Mix | 17.5 μl |
| Ligation Enhancer | 0.5 μl |
| 總量 | 35.5 μl |
NBXX barcodes 為 EXP-NBD104(01~12)和EXP-NBD114(13～24) 試劑中的 barcodes
- 室溫放置15min。置于PCR儀上，70°C 10min，立即置于冰上至少1min
70°C 10min 目的是抑制DNA Ligase活性，避免接頭形成二聚體。
- 將連接barcodes后的所有樣本的混合到一個EP管中。取1uL進行定量
- 取LoBind Tubes EP管，按照下表進行配液，室溫放置15min。
| 試劑 | 量 |
|---|---|
| Barcoded amplicon pools | 30 μl |
| NEBNext Quick Ligation Reaction Buffer (5X) | 10 μl |
| AMII adapter mix | 5 μl |
| Quick T4 DNA Ligase | 5 μl |
| 總量 | 50 μl |
input DNA的量根據(jù)barcods數(shù)量決定，數(shù)量在40ng～160ng(8～24 barcods)。
- 加入50uL AMPure XP Beads，輕輕混勻。室溫放置5min
- 置于磁力架上2min，液體澄清后吸棄液體。
- 磁力架上取下EP管，加入200uL SFB，重懸磁珠。
- EP管放到磁力架上靜置2min，液體澄清后吸棄。用SFB再洗1次。
- 加入15uL EB重懸磁珠，室溫放置2min。
- 將EP管置于磁力架上。吸取澄清液體到一個新的EP管中，為測序文庫。
用于上機測序的文庫DNA量為20ng
- 將30uL FLT 加入到1管FB中，震蕩混勻。
- 打開測序芯片的priming port，用1mL移液器垂直頂住priming port，慢慢旋轉移液器量程，黃色納米孔保護液，直到沒有氣泡。
- 吸取800uL FLT/FB混合液，從Priming port中加入，可以不用加完，避免氣泡加入芯片中。等待5min。
- 打開SpotOn Port，從Priming Port中加入200uL FLT/FB混合液。可也看到SpotOn Port 中有液滴因為Priming Port加混合液而鼓起。
- 取一個EP管，按照下表配液：
| 試劑 | 量 |
|---|---|
| SQB | 37.5 μl |
| LB | 25.5 μl |
| Final library | 12 μl |
| 總量 | 75 μl |
LB使用前再充分混勻
- 滴加75uL文庫到SpotOn Port中，液體完全吸收后關閉Spoton Port和Priming Port，將芯片裝入測序儀，打開MinKnow測序軟件進行測序。
對于沒用用完nanopores的芯片，可以啟用清洗程序以后繼續(xù)使用：
- 取 1管 Wash Solution B置于室溫，震蕩混勻后放在冰上待用。
- 在 1個 EP管中，加入20uL Wash Solution A和380uL Wash Solution B。吹吸混勻（不要震蕩），置于冰上。
- 暫停測序，但不要取出芯片。
- 確定Priming Port和SpotON Port關閉。用1000uL移液器，從Waste Port1中吸出所有液體，確定芯片納米孔處沒有殘液。
- 打開Priming Port，先確認消除氣泡。然后加入400Wash Solution洗液，注意避免加入氣泡。
- 關閉Priming Port，等30min后，用1000uL移液器從Waste Port1中吸棄所有液體。
- 如果還要測其他樣本，可參考上樣流程操作。如果暫時不用，則按照下面流程進行保存芯片。
- 將一管Storage Buffer(S)取出置于室溫，溶解后上下顛倒混勻。
- 打開用完的芯片Priming Port，1000uL移液器設置成200uL量程，然后槍頭Priming Port頂住孔，旋轉移液器量程20uL，以確定排除氣泡。
- 吸取500uL Storage Buffer(S)，從Priming Port中緩慢加入。關閉Priming Port。
- 關閉Priming Port和SpotON Port，用1000uL移液器從Waste Port1中吸取所有液體。
- 芯片置于冰箱冷藏。
下一次使用的芯片，最好設置成不同的電壓以獲得更好的結果，參見
4.2 illumina Miseq 測序文庫
4.2.1 試劑與耗材
- Nextera XT LIbrary Preparation Kit
- NEB Ultra II DNA Library Prep Kit for Illumina
- AMPure XP Beads
- Miseq Reagent v2 PE150
對于濃度較高的樣本，illumina擴增子測序推薦使用PCR-Free的方式進行。這里用NEB Ultra II DNA Library Prep Kit for Illumina(E7645)文庫構建試劑盒進行。如果使用Nextera XT之類的酶片段化試劑盒，后期數(shù)據(jù)分析時要考慮引物擴增位點對于consensus序列的影響。
4.2.2 實驗流程
- 將 pool1 和 pool2 的 PCR 產物混合后，使用Qubit進行定量。NEB文庫構建試劑盒適合500pg~1ug的起始input DNA量。
- 取新的 PCR 管，按照下表配液。用200uL槍的50uL量程吹吸10次混勻。離心甩干。
| 試劑 | 量 |
|---|---|
| NEBNext Ultra II End Prep Enzyme Mix | 3 μl |
| NEBNext Ultra II End Prep Reaction Buffer | 7 μl |
| DNA | 50 μl |
| 總量 | 60 μl |
- 放置于PCR儀器上，熱蓋溫度大于等于75°C ，運行以下程序：20°C 30min，65°C 30min，4°C hold。
- 按照下表稀釋 adaptor。建議初始DNA量大于100ng，如200ng左右。
| Input DNA | Adaptor | Adaptor 工作濃度 |
|---|---|---|
| 1ug-101ng | 不稀釋 | 15uM |
| 100ng-5ng | 1:10稀釋 | 1.5 uM |
| <5ng | 1:25稀釋 | 0.6uM |
- 根據(jù)下表配置Adaptor連接體系。NEBNext Ultra II Ligation Master Mix加入前要顛倒混勻。
| 試劑 | 量 |
|---|---|
| 前一步反應液 | 60 μl |
| NEBNext Ultra II Ligation Master Mix | 30 μl |
| NEBNext LIgation Enhancer | 1 μl |
| NEBNext Adaptor for illumina | 2.5 μl |
| 總量 | 93.5 μl |
可以采用的Adaptor產品：單端（E7350），雙端（E7335,E7500,E7710,E7730,E7600,E7535,E6609）
- 用200uL槍的80uL量程吹戲10次混勻。離心甩干。
- 置于金屬浴或PCR儀上（不開熱蓋）20°C 15min。
- 向反應體系中加入3uL USER Enzyme（這個試劑包含在接頭試劑盒中）。在PCR儀器上37°C 15min（熱蓋大于等于47°C）。
- 將液體轉移到深孔板內，當input DNA大于50ng時，進行下面的DNA純化回收。對于小于50ng DNA的樣本，只需要用磁株純化1次，用87uL量。
- 加入18uL AMPure Beads到反應液中，吹吸10次混勻，注意槍頭中的液體要全部打出。室溫孵育5min。
- 將深孔板置于磁力架上，靜置5min待液體澄清。吸出液體到新的孔中。
- 將深孔板從磁力架上取下，加入10uL AMPure Beads，吹吸10次混勻后，室溫孵育5min。
- 將深孔板置于磁力架上，靜置5min待液體澄清。將液體吸棄。
- 加入80%新鮮配置的乙醇200uL，室溫放置30s后吸棄。重復洗1次。
- 開蓋風干磁珠5min。如果風干時間過長，會導致DNA回收率下降。
- 從磁力架上取下深孔板，17uL 10mM Tris-HCl或0.1xTE洗脫。
- 震蕩儀上1800rpm震蕩10min，室溫放置2min。
- 將深孔板置于磁力架上，5min后待液體澄清，取1uL進行Qubit定量。
- 吸取15uL液體到新的PCR板。
5. 數(shù)據(jù)分析
5.1 Nanopore 測序數(shù)據(jù)
5.1.1 軟件安裝
使用 conda 構建獨立的運行環(huán)境。
$ git clone --recursive https://github.com/artic-network/artic-ncov2019.git
$ conda env create -n ncov -f artic-ncov2019/environment.yml
$ conda activate ncov
(ncov)$ conda list
如果使用docker images來運行guppy
$ docker pull genomicpariscentre/guppy
# 如果服務器有GPU支持，可以選擇guppy-gpu
$ docker pull genomicpariscentre/guppy-gpu
# 交互模式運行container
$ docker run --name my_exp -it -v local_dir_of_fast5:/media/ genomicpariscentre/guppy /bin/bash
# 對fast5數(shù)據(jù)進行basecalling
$ guppy_basecaller -i $HOME/data/fast5 -s output_folder \
> -c dna_r9.4.1_450bps_hac.cfg --barcode_kits EXP-NBD104 \
> --num_callers 40 --trim_barcodes
5.1.2 分析流程
與illumina原始數(shù)據(jù)是圖像文件不同，Nanopore的原始數(shù)據(jù)以HDF5格式保存。HDF是"Hierarchical Data Format"首字母縮寫，HDF5是一種純文本，類似json的數(shù)據(jù)格式記錄電信號。查看HDF5原始數(shù)據(jù)可以用 hdf5_tools 工具。數(shù)據(jù)文件后綴一般用.fast5表示。
Nanopore 測序數(shù)據(jù)期初由于其錯誤率高而“聞名”，特別是對二聚體的區(qū)分。電信號數(shù)據(jù)隨著處理軟件和算法的不斷改進，準確率得到不斷的提升。HDF5數(shù)據(jù)進行basecalling的軟件目前很多，采用的算法也各異。官方的如albacore，guppy等，第三方工具也有很多。
MinKnow軟件 將fast5數(shù)據(jù)復制到服務器，使用guppy進行分析
guppy是一個是用
我們的實驗室用一臺Workstation掛載MinION測序儀，進行測序工作?；九渲檬莍7-7700k+16G+1TSSD/1T HD。一般可以做到70%以上的MinKnow實時basecalling。但是如果要使用guppy進行更準確的basecalling的話，就需要將數(shù)據(jù)同步到服務器上進行。
# 在服務器上 rsync 同步數(shù)據(jù)到本地
(ncov)$ rsync -avz username@minknow_ip:/var/lib/minknown/path/to/fast5 .
# 或者在工作站上同步數(shù)據(jù)到服務器端
$ rsync -avz /var/lib/minknown/path/to/fast5 username@server_ip:~/data/fast5
這里設置服務器端 fast5 數(shù)據(jù)目錄存放于：$HOME/data/fast5；采用的是R9.4的測序芯片，各個不同測序芯片basecalling設置參數(shù)可以參見這里
(ncov)$ guppy_basecaller -c dna_r9.4.1_450bps_fast.cfg \
> -i $HOME/data/fast5 -s run_name -x auto -r
我們的MinION測序儀連接的是一臺Ubuntu 16.04的Workstation，測序儀最大測序狀態(tài)時，跑默認的basecalling基本上可以實時達到80%的狀態(tài)。如果需要用guppy做更準確的basecalling，我們會用inotify-tools之類的工具加上 rsync 實時同步到服務器，在服務器上進行。
artic 的 demultiplex 命令是
$ porechop --verbosity 2 --untrimmed -i "re_all.fastq" -b ./tmp5vx2iwn0 --native_barcodes --discard_middle --require_two_barcodes --barcode_threshold 80 --threads 8 --check_reads 10000 --barcode_diff 5 > re_all.fastq.demultiplexreport.txt
5.2 Illumina 測序數(shù)據(jù)
5.2.1 軟件安裝
暴發(fā)疫情中病毒毒株發(fā)生結構變異的可能性較小，可以直接使用比對參考基因組的方法獲得毒株基因組序列。如果是研究新發(fā)的遠源病毒或長時間演化的病毒序列，可能與參考基因組差異較大時，應結合組裝與比對的方法獲得基因組序列。
$ conda create -n mapping
$ conda activate mapping
(mapping)$ conda install bwa samtools igv
5.2.2 分析流程
- 使用Miseq自帶的SAV軟件進行basecalling，生成fastq數(shù)據(jù)復制到服務器端
- 可以采用bwa+samtools進行參考序列比對獲得一致性序列，或者采用bwa等工具比對獲得新冠病毒序列后在進行de novo組裝。
# 最基本的比對分析流程
(mapping)$ bwa index MN909847.3.fasta
(mapping)$ bwa mem -k 45 -R "@RG\tID:HZ-1_1\tSM:HZ-1" MN909847.3.fasta \
> HZ-1_S1_R1.fastq.gz HZ-1_S1_R2.fastq.gz | samtools view -bS > HZ-1.bam
(mapping)$ samtools sort -O bam -o HZ-1.sorted.bam -@4 HZ-1.bam
(mapping)$ samtools index HZ-1.sorted.bam
(mapping)$ igv HZ-1.sorted.bam
# 以bwa為例進行reads過濾后spades組裝
$ bwa index MN908947.3.fasta
$ bwa mem -k 47 -R "@RG\tID:HZ-1\tSM:HZ-1" MN908947.3.fasta \
> HZ-1_S1_L001_R1_001.fastq.gz HZ-1_S1_L001_R2_001.fastq.gz -t 40 | \
> samtools view -bS - | samtools sort -@40 -O bam -o HZ-1.sorted.bam
$ samtools index HZ-1.sorted.bam
$ samtools faidx MN908947.3.fasta
$ bcftools mpileup -f MN908047.3.fasta HZ-1.sorted.bam | \
> bcftools call --ploidy 1 -mv -Oz > HZ-1.vcf.gz
$ tabix HZ-1.vcf.gz
$ cat MN908947.3.fasta | bcftools consensus HZ-1.vcf.gz > Hz-1_consensus.fasta
# de novo 組裝
$ bwa index MN908947.3.fasta
$ bwa mem -k 45 -R "@RG\tID:HZ-1\tSM:HZ-1" MN908947.3.fasta \
> HZ-1_S1_L001_R1_001.fastq.gz HZ-1_S1_L001_R2_001.fastq.gz -t 4 | \
> samtools view -bS > HZ-1.bam
# 提取雙端都比對到參考序列的reads
$ samtools view -bF 12 HZ-1.bam > HZ-1.mapped.bam
# bam格式轉化成fastq格式
$ bamToFastq -i HZ-1.mapped.bam -fq HZ-1.mapped_R1.fastq -fq2 HZ-1.mapped_R2.fastq
# de novo assembly
$ shovill --trim --outdir test --R1 HZ-1.mapped_R1.fastq -fq2 HZ-1.mapped_R2.fastq --ram 16 --cpus 40
$ bwa mem -t 40 -x ont2d MN908947.3.fasta test/contigs.fa | samtools view -bS - | samtools sort -o test.sorted.bam -
$ samtools index test.sorted.bam
自定義分析流程
針對MinKnow生成的數(shù)據(jù)
# 生成樣本的 fastq
$ for i in *.fastq; do zcat $i >> s1.fq; done
$ gzip s1.fq
$ ls *.fastq | parallel --max-args=1 cat {1} | gzip > all.fq.gz
# 去除可能的接頭序列
$ porechop -t 4 -i s1.fq.gz --format fastq.gz -o s1.clean.fq.gz
# 去除低質量序列
$ gunzip -c s1.clean.fq.gz | NanoFilt -q 10 -l 400 --maxlength 700 | gzip > s1.clean.hq.fq.gz
# 獲得比對結果
$ bwa index MN908947.3.fasta
$ bwa mem -t 4 -x ont2d MN908947.3.fasta s1.clean.hq.fq.gz | \
> samtools view -bS - | \
> samtools sort -o s1.clean.hq.sorted.bam -
$ samtools index s1.clean.hq.sorted.bam
$ samtools faidx MN908947.3.fasta
# vcf calling
$ bcftools mpileup --threads 4 --fasta-ref MN908047.3.fasta s1.clean.hq.sorted.bam | \
> bcftools call --ploidy 1 -mv -Oz > s1.vcf.gz
$ tabix s1.vcf.gz
# 獲得 consensus 序列
$ bcftools consensus -f MN908947.3.fasta -H 1 s1.vcf.gz -o s1_consensus.fasta
# 或者使用minimap工具進行比對
