本文檔采用 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)行比對