本文檔采用 ARTIC Network 流程進(jìn)行實(shí)驗(yàn)错邦，可用于illumina平臺和nanopore平臺對SARS-Cov-2病毒進(jìn)行基因組測序佛点。
@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進(jìn)行測序也可以成功，在目標(biāo)reads豐度上SSIV似乎更好一些劈榨。
1.2 實(shí)驗(yàn)流程
- 采用RNA試劑盒手工提取樣本總RNA访递，也可以用機(jī)器自動提取。
RNA通過Qiagen RNeasy Mini Kit 手工提取試劑盒或者機(jī)器提取的都可以同辣，PCR實(shí)驗(yàn)結(jié)果CT無明顯差異拷姿。有文章比較過Qiagen不同的RNA提取試劑盒的宏基因組測序reads產(chǎn)出效果，可以供大家參考旱函。
- 取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 |
注意事項(xiàng)：樣本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蝌箍，相當(dāng)與摩爾數(shù) 25e-6μmol，因此根據(jù)表中量暴心，需要加2μl十绑。但根據(jù)試劑盒說明書中也是加1μl量，因此這里參考說明書酷勺。
推薦加入樣本前，在PCR儀上預(yù)熱扳躬。
- 反應(yīng)體系在PCR儀器中 65°C 5min脆诉。然后立即置于冰上至少1min。
- 按下表配置反應(yīng)體系贷币，添加到冰上的PCR管中击胜。置于PCR儀器中，42°C 50min役纹，70°C 10min偶摔，結(jié)束后5°C 保溫。這一步要設(shè)置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管內(nèi)液體 | 13 μl | 13 μl * 8 |
| Total | 20 μl | 20 μl * 8 |
以上配液操作在PCR工作站或生物安全柜內(nèi)進(jìn)行辰斋。使用前紫外處理30分鐘并通風(fēng)。
2. PCR擴(kuò)增
2.1 試劑與耗材
ARTIC Network流程里使用的是NEB Q5酶瘸味，試用Qiagen的多重PCR產(chǎn)品宫仗，效果也可以。
2.2 實(shí)驗(yàn)流程
- 將98組凍干引物 8000rpm 離心 10min旁仿。將每個(gè)引物配置成100uM濃度藕夫。
- 準(zhǔn)備2個(gè)EP管，標(biāo)記為P1和P2枯冈。將98組引物分為奇數(shù)和偶數(shù)組毅贮，每個(gè)組各98個(gè)100uM的引物，每管吸取5uL分別到Pool1和2中尘奏，最終每管含100uM濃度的引物共980uL滩褥。最終混合引物中每條引物濃度約為1uM。
- 震蕩混勻罪既，離心甩干后铸题，各取20uL P1/2 引物铡恕，加入180uL水。使得Pool1和Pool2引物終濃度為10uM丢间，每條引物終濃度約為0.1uM探熔。
以上配液操作在PCR工作站或生物安全柜內(nèi)進(jìn)行。使用前紫外處理30分鐘并通風(fēng)烘挫。
- 每個(gè)樣本取2個(gè)PCR管诀艰，使用NEB Q5 Hot Start DNA Polymerase試劑或者Qiagen Multiplex PCR試劑進(jìn)行PCR擴(kuò)增，配液參見表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)化其垄，每個(gè)引物的終濃度為0.015uM。
- NEB Q5 的PCR程序?yàn)椋?/li>
| 序號 | 溫度 | 時(shí)間 | 前往 |
|---|---|---|---|
| 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程序?yàn)椋?/li>
| 序號 | 溫度 | 時(shí)間 | 前往 |
|---|---|---|---|
| 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個(gè)循環(huán)卤橄。CT值介于21～35的根據(jù)分布自行設(shè)置绿满。
應(yīng)用本方法到其他病原檢測時(shí)，需要設(shè)計(jì)多重PCR的話窟扑，可以使用這里的引物設(shè)計(jì)工具獲得引物組喇颁。
3. Amplicon純化回收
3.1 試劑與耗材
3.2 實(shí)驗(yàn)流程
- 將每個(gè)樣本的Pool1和Pool2混合，因?yàn)闇y序文庫與上樣量有關(guān)嚎货，本實(shí)驗(yàn)流程設(shè)計(jì)的是針對8個(gè)以上樣本的上樣量橘霎。建議檢測樣本帶質(zhì)控DNA做為陰性對照，以便檢查是否有污染序列殖属。
- 將AMPure XP Beads 提前30min取出置于常溫姐叁，使用前在震蕩儀上混合均勻。
- 取50uL AMPure XP Beads 到 50uL 混合pooling Amplicon中洗显。震蕩混勻（增加震蕩時(shí)間能極大提高回收率）外潜，離心甩干。
1x 磁珠吸附250～1000bp效率較好墙懂。具體參見
- 室溫放置5min橡卤。
- 將EP管放置于磁力架上2min，直至液體澄清损搬。
- 吸棄液體碧库，用200uL新鮮配置的70～80%乙醇洗2次。
- 用10uL移液器吸棄殘液巧勤，開蓋1min嵌灰，讓乙醇揮發(fā)。
- 從磁力架上取下颅悉，用15～30uL Elution Buffer或者水重懸磁珠沽瞭，用槍頭或手指輕輕混勻，靜置2min剩瓶。
- 置于磁力架上驹溃，液體澄清后吸取30uL液體到一個(gè)新的EP管中城丧。取1uL進(jìn)行Qubit定量。
樣本濃度與PCR產(chǎn)物回收相關(guān)性：
- 10E5大約200ng/uL
- 10E4大約120ng/uL
- 10E3大約80ng/uL
- 10E2大約25ng/uL
- 10E1大約10ng/uL
4. 文庫構(gòu)建
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 實(shí)驗(yàn)流程
注意此流程是根據(jù)6~24個(gè)樣本混樣的量規(guī)劃的豌鹤。如果只做單個(gè)樣本亡哄，可以不加barcodes，但需要調(diào)整input DNA量布疙，建議加入DNA達(dá)40ng蚊惯，保證上機(jī)時(shí)能約有20ng DNA，達(dá)到最佳的芯片測序分子數(shù)灵临。
- 將定量的DNA稀釋成1ng/uL截型，該濃度是針對700bp長度的Amplicon，如果片段長400bp儒溉，建議濃度為2ng/uL宦焦，這樣分子濃度約為100～200fmol。
- 取0.2mLPCR管顿涣，根據(jù)下表進(jìn)行配液赶诊，進(jìn)行末端修復(fù)。
| 試劑 | 量 |
|---|---|
| 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后的所有樣本的混合到一個(gè)EP管中痛悯。取1uL進(jìn)行定量
- 取LoBind Tubes EP管余黎，按照下表進(jìn)行配液，室溫放置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管置于磁力架上。吸取澄清液體到一個(gè)新的EP管中绊寻，為測序文庫花墩。
用于上機(jī)測序的文庫DNA量為20ng
- 將30uL FLT 加入到1管FB中，震蕩混勻澄步。
- 打開測序芯片的priming port冰蘑，用1mL移液器垂直頂住priming port，慢慢旋轉(zhuǎn)移液器量程村缸，黃色納米孔保護(hù)液祠肥，直到?jīng)]有氣泡。
- 吸取800uL FLT/FB混合液梯皿，從Priming port中加入仇箱，可以不用加完，避免氣泡加入芯片中东羹。等待5min剂桥。
- 打開SpotOn Port，從Priming Port中加入200uL FLT/FB混合液属提∪ǘ海可也看到SpotOn Port 中有液滴因?yàn)镻riming Port加混合液而鼓起。
- 取一個(gè)EP管冤议，按照下表配液：
| 試劑 | 量 |
|---|---|
| SQB | 37.5 μl |
| LB | 25.5 μl |
| Final library | 12 μl |
| 總量 | 75 μl |
LB使用前再充分混勻
- 滴加75uL文庫到SpotOn Port中斟薇，液體完全吸收后關(guān)閉Spoton Port和Priming Port，將芯片裝入測序儀恕酸，打開MinKnow測序軟件進(jìn)行測序堪滨。
對于沒用用完nanopores的芯片，可以啟用清洗程序以后繼續(xù)使用：
- 取 1管 Wash Solution B置于室溫蕊温，震蕩混勻后放在冰上待用袱箱。
- 在 1個(gè) EP管中，加入20uL Wash Solution A和380uL Wash Solution B义矛。吹吸混勻（不要震蕩）犯眠，置于冰上。
- 暫停測序症革，但不要取出芯片筐咧。
- 確定Priming Port和SpotON Port關(guān)閉。用1000uL移液器，從Waste Port1中吸出所有液體量蕊，確定芯片納米孔處沒有殘液铺罢。
- 打開Priming Port，先確認(rèn)消除氣泡残炮。然后加入400Wash Solution洗液韭赘，注意避免加入氣泡。
- 關(guān)閉Priming Port势就，等30min后泉瞻，用1000uL移液器從Waste Port1中吸棄所有液體。
- 如果還要測其他樣本苞冯，可參考上樣流程操作袖牙。如果暫時(shí)不用，則按照下面流程進(jìn)行保存芯片舅锄。
- 將一管Storage Buffer(S)取出置于室溫鞭达，溶解后上下顛倒混勻。
- 打開用完的芯片Priming Port皇忿，1000uL移液器設(shè)置成200uL量程畴蹭，然后槍頭Priming Port頂住孔，旋轉(zhuǎn)移液器量程20uL鳍烁，以確定排除氣泡叨襟。
- 吸取500uL Storage Buffer(S)，從Priming Port中緩慢加入幔荒。關(guān)閉Priming Port芹啥。
- 關(guān)閉Priming Port和SpotON Port，用1000uL移液器從Waste Port1中吸取所有液體铺峭。
- 芯片置于冰箱冷藏。
下一次使用的芯片汽纠，最好設(shè)置成不同的電壓以獲得更好的結(jié)果卫键，參見
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擴(kuò)增子測序推薦使用PCR-Free的方式進(jìn)行虱朵。這里用NEB Ultra II DNA Library Prep Kit for Illumina(E7645)文庫構(gòu)建試劑盒進(jìn)行莉炉。如果使用Nextera XT之類的酶片段化試劑盒，后期數(shù)據(jù)分析時(shí)要考慮引物擴(kuò)增位點(diǎn)對于consensus序列的影響碴犬。
4.2.2 實(shí)驗(yàn)流程
- 將 pool1 和 pool2 的 PCR 產(chǎn)物混合后絮宁，使用Qubit進(jìn)行定量。NEB文庫構(gòu)建試劑盒適合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 ，運(yùn)行以下程序：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加入前要顛倒混勻亮靴。
| 試劑 | 量 |
|---|---|
| 前一步反應(yīng)液 | 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產(chǎn)品：單端（E7350），雙端（E7335,E7500,E7710,E7730,E7600,E7535,E6609）
- 用200uL槍的80uL量程吹戲10次混勻敌厘。離心甩干台猴。
- 置于金屬浴或PCR儀上（不開熱蓋）20°C 15min。
- 向反應(yīng)體系中加入3uL USER Enzyme（這個(gè)試劑包含在接頭試劑盒中）俱两。在PCR儀器上37°C 15min（熱蓋大于等于47°C）饱狂。
- 將液體轉(zhuǎn)移到深孔板內(nèi)，當(dāng)input DNA大于50ng時(shí)宪彩，進(jìn)行下面的DNA純化回收休讳。對于小于50ng DNA的樣本，只需要用磁株純化1次尿孔，用87uL量俊柔。
- 加入18uL AMPure Beads到反應(yīng)液中，吹吸10次混勻活合，注意槍頭中的液體要全部打出雏婶。室溫孵育5min。
- 將深孔板置于磁力架上白指，靜置5min待液體澄清留晚。吸出液體到新的孔中。
- 將深孔板從磁力架上取下告嘲，加入10uL AMPure Beads错维，吹吸10次混勻后，室溫孵育5min橄唬。
- 將深孔板置于磁力架上赋焕，靜置5min待液體澄清。將液體吸棄仰楚。
- 加入80%新鮮配置的乙醇200uL隆判，室溫放置30s后吸棄犬庇。重復(fù)洗1次。
- 開蓋風(fēng)干磁珠5min蜜氨。如果風(fēng)干時(shí)間過長械筛，會導(dǎo)致DNA回收率下降。
- 從磁力架上取下深孔板飒炎，17uL 10mM Tris-HCl或0.1xTE洗脫埋哟。
- 震蕩儀上1800rpm震蕩10min，室溫放置2min郎汪。
- 將深孔板置于磁力架上赤赊，5min后待液體澄清，取1uL進(jìn)行Qubit定量煞赢。
- 吸取15uL液體到新的PCR板抛计。
5. 數(shù)據(jù)分析
5.1 Nanopore 測序數(shù)據(jù)
5.1.1 軟件安裝
使用 conda 構(gòu)建獨(dú)立的運(yùn)行環(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來運(yùn)行g(shù)uppy
$ docker pull genomicpariscentre/guppy
# 如果服務(wù)器有GPU支持照筑，可以選擇guppy-gpu
$ docker pull genomicpariscentre/guppy-gpu
# 交互模式運(yùn)行container
$ docker run --name my_exp -it -v local_dir_of_fast5:/media/ genomicpariscentre/guppy /bin/bash
# 對fast5數(shù)據(jù)進(jìn)行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ù)期初由于其錯(cuò)誤率高而“聞名”，特別是對二聚體的區(qū)分支鸡。電信號數(shù)據(jù)隨著處理軟件和算法的不斷改進(jìn)冬念，準(zhǔn)確率得到不斷的提升。HDF5數(shù)據(jù)進(jìn)行basecalling的軟件目前很多牧挣，采用的算法也各異急前。官方的如albacore，guppy等瀑构，第三方工具也有很多裆针。
MinKnow軟件 將fast5數(shù)據(jù)復(fù)制到服務(wù)器，使用guppy進(jìn)行分析
guppy是一個(gè)是用
我們的實(shí)驗(yàn)室用一臺Workstation掛載MinION測序儀检碗，進(jìn)行測序工作÷肓冢基本配置是i7-7700k+16G+1TSSD/1T HD折剃。一般可以做到70%以上的MinKnow實(shí)時(shí)basecalling。但是如果要使用guppy進(jìn)行更準(zhǔn)確的basecalling的話像屋，就需要將數(shù)據(jù)同步到服務(wù)器上進(jìn)行怕犁。
# 在服務(wù)器上 rsync 同步數(shù)據(jù)到本地
(ncov)$ rsync -avz username@minknow_ip:/var/lib/minknown/path/to/fast5 .
# 或者在工作站上同步數(shù)據(jù)到服務(wù)器端
$ rsync -avz /var/lib/minknown/path/to/fast5 username@server_ip:~/data/fast5
這里設(shè)置服務(wù)器端 fast5 數(shù)據(jù)目錄存放于：$HOME/data/fast5；采用的是R9.4的測序芯片，各個(gè)不同測序芯片basecalling設(shè)置參數(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)時(shí)戈轿，跑默認(rèn)的basecalling基本上可以實(shí)時(shí)達(dá)到80%的狀態(tài)。如果需要用guppy做更準(zhǔn)確的basecalling阵子，我們會用inotify-tools之類的工具加上 rsync 實(shí)時(shí)同步到服務(wù)器思杯，在服務(wù)器上進(jìn)行。
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ā)生結(jié)構(gòu)變異的可能性較小挠进，可以直接使用比對參考基因組的方法獲得毒株基因組序列色乾。如果是研究新發(fā)的遠(yuǎn)源病毒或長時(shí)間演化的病毒序列，可能與參考基因組差異較大時(shí)领突，應(yīng)結(jié)合組裝與比對的方法獲得基因組序列暖璧。
$ conda create -n mapping
$ conda activate mapping
(mapping)$ conda install bwa samtools igv
5.2.2 分析流程
- 使用Miseq自帶的SAV軟件進(jìn)行basecalling，生成fastq數(shù)據(jù)復(fù)制到服務(wù)器端
- 可以采用bwa+samtools進(jìn)行參考序列比對獲得一致性序列君旦，或者采用bwa等工具比對獲得新冠病毒序列后在進(jìn)行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為例進(jìn)行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格式轉(zhuǎn)化成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
# 去除低質(zhì)量序列
$ gunzip -c s1.clean.fq.gz | NanoFilt -q 10 -l 400 --maxlength 700 | gzip > s1.clean.hq.fq.gz
# 獲得比對結(jié)果
$ 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工具進(jìn)行比對
