SARS-CoV-2 Genome Sequencing Protocol

本文檔采用 ARTIC Network 流程進(jìn)行實(shí)驗(yàn)错邦,可用于illumina平臺和nanopore平臺對SARS-Cov-2病毒進(jìn)行基因組測序佛点。

@indexofire indexofire@gmail.com

1. cDNA獲取

1.1 試劑與耗材

SSIV的酶無RNaseH活性灰伟,獲得長片段(~12kb)的全長cDNA效果更好。我們使用SSIII蓉驹,有中等活性的RNaseH活性埂蕊,獲得的cDNA進(jìn)行測序也可以成功,在目標(biāo)reads豐度上SSIV似乎更好一些劈榨。

不同酶的關(guān)鍵逆轉(zhuǎn)錄性能比較

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 試劑與耗材

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剂桥。
芯片Priming port和SpotON port Cover位置介紹
  • 打開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)行比對

參考資料

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個(gè)濱河市金砍,隨后出現(xiàn)的幾起案子局蚀,更是在濱河造成了極大的恐慌,老刑警劉巖捞魁,帶你破解...
    沈念sama閱讀 223,002評論 6 519
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件至会,死亡現(xiàn)場離奇詭異,居然都是意外死亡谱俭,警方通過查閱死者的電腦和手機(jī)奉件,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,357評論 3 400
  • 文/潘曉璐 我一進(jìn)店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來昆著,“玉大人县貌,你說我怎么就攤上這事〈斩” “怎么了煤痕?”我有些...
    開封第一講書人閱讀 169,787評論 0 365
  • 文/不壞的土叔 我叫張陵,是天一觀的道長接谨。 經(jīng)常有香客問我摆碉,道長,這世上最難降的妖魔是什么脓豪? 我笑而不...
    開封第一講書人閱讀 60,237評論 1 300
  • 正文 為了忘掉前任巷帝,我火速辦了婚禮,結(jié)果婚禮上扫夜,老公的妹妹穿的比我還像新娘楞泼。我一直安慰自己驰徊,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 69,237評論 6 398
  • 文/花漫 我一把揭開白布堕阔。 她就那樣靜靜地躺著棍厂,像睡著了一般。 火紅的嫁衣襯著肌膚如雪超陆。 梳的紋絲不亂的頭發(fā)上牺弹,一...
    開封第一講書人閱讀 52,821評論 1 314
  • 那天,我揣著相機(jī)與錄音侥猬,去河邊找鬼例驹。 笑死,一個(gè)胖子當(dāng)著我的面吹牛退唠,可吹牛的內(nèi)容都是我干的鹃锈。 我是一名探鬼主播,決...
    沈念sama閱讀 41,236評論 3 424
  • 文/蒼蘭香墨 我猛地睜開眼瞧预,長吁一口氣:“原來是場噩夢啊……” “哼屎债!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起垢油,我...
    開封第一講書人閱讀 40,196評論 0 277
  • 序言:老撾萬榮一對情侶失蹤盆驹,失蹤者是張志新(化名)和其女友劉穎,沒想到半個(gè)月后滩愁,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體躯喇,經(jīng)...
    沈念sama閱讀 46,716評論 1 320
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,794評論 3 343
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片舍扰。...
    茶點(diǎn)故事閱讀 40,928評論 1 353
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡粱挡,死狀恐怖防泵,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情,我是刑警寧澤,帶...
    沈念sama閱讀 36,583評論 5 351
  • 正文 年R本政府宣布焦履,位于F島的核電站,受9級特大地震影響雏逾,放射性物質(zhì)發(fā)生泄漏嘉裤。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,264評論 3 336
  • 文/蒙蒙 一栖博、第九天 我趴在偏房一處隱蔽的房頂上張望屑宠。 院中可真熱鬧,春花似錦笛匙、人聲如沸侨把。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,755評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽秋柄。三九已至,卻和暖如春蠢正,著一層夾襖步出監(jiān)牢的瞬間骇笔,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,869評論 1 274
  • 我被黑心中介騙來泰國打工嚣崭, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留笨触,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 49,378評論 3 379
  • 正文 我出身青樓雹舀,卻偏偏與公主長得像芦劣,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個(gè)殘疾皇子说榆,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,937評論 2 361

推薦閱讀更多精彩內(nèi)容