使用minimap+miniasm对nanopore进行基因组组装

简介: 我们用来练手的文章发表在 Nature Communication ,"High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell", 非常不要脸的说,这篇文章是我师爷实验室发的。

我们用来练手的文章发表在 Nature Communication ,"High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell", 非常不要脸的说,这篇文章是我师爷实验室发的。

简单讲讲故事内容,就是他们实验室买了一台nanopore仪器,就是下面这台, 目前仪器价格国内是8K左右,当然测序的价格就另说了。如同买台PS4主机,还要买游戏,买个单反,你还得买镜头。仪器只是败家的开始!

img_53828e1e76a31220a3f5d9239cb8eec6.jpe
nanopore

他们认为三代测序目前有两大问题,测的还不够长以及不够准。nanopore解决了其中一个问题,不够长。Arabidopsis thaliana 当年用一代测序,虽然可以认为是组装的金标准了,但是还是有很多区域是BAC连BAC文库搞不定的,所以就用这台仪器把 Arabidopsis thaliana 测了一波。显然就测一个nanopore,还是已知序列的物种是不可能发文章的,于是他们又用Pacbio sequel测了一波。最后用bionano 光学图谱验证了一次(请大家自行计算要多少钱)。

光测序不行,还得组装对吧。传统的组装方法是想办法利用高深度和随机错误进行纠错,然后用纠错后的长序列进行组装,最后用二代进行纠错。对于一台不错的服务器(20W起步吧)大约花个十天半个月就行。作者或许认为买一台20多w的外设配合不到1w的测序仪可能是太蠢了,于是他用了比较Li Heng大神开发的工具,Minimap+miniasm进行组装,然后用racon+pillon进行纠错,用了一台Macbook Pro 15.6寸花了4天就搞定了,并且和常规工具比较,还算过得去哦。

下面就是正式的分析:

根据文章提供的项目编号"PRJEB21270", 在European Nucleotide Archive上找到下载地址。

img_3a585281e5241e7bf5123716dac19d3d.jpe
ENA搜索

进入这个页面之后,就可以去下载作者用到的所有数据,我们下载Sequel和MinIon和Illuminia的数据就好了,数据量加起来差不多30G。

## Sequal
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam.bai
## MinION
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116595/fastq/ont.fq.gz
# Illuminia MiSeq
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_1.fq.gz
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_2.fq.gz

拿到数据之后,我们就可以用作者提供的分析流程进行重复了。地址为https://github.com/fbemm/onefc-oneasm/wiki/Assembly-Generation

这就是大神的自信,把代码都给你,反正你也看不懂。当然我在重复的时候用的都是最新的软件,所以会有所不同

第一步:拿着80%~90%正确率的原始数据相互比对, 找序列之间的Overlap。这一步,我花了30分钟

time ~/opt/biosoft/minimap2/minimap2 -t 10 -x ava-ont ont.fq ont.fq > gzip -1 ont.paf.gz &

第二步:找到Overlap,就能够进行组装了。这一步我花了2分钟

time ~/opt/biosoft/miniasm/miniasm -f ont.fq ont.paf > ONTmin.gfa &
awk '/^S/{print ">"$2"\n"$3}' ONTmin.gfa | seqkit seq > ONTmin_IT0.fasta &

第三步: 原始的组装结果充满了错误,所以需要进行纠错。纠错分为两种,一种是用三代自身数据,一种是用二代数据进行纠错。当然这两步都是需要的

首先使用三代数据进行纠错,古语有云“事不过三”一般迭代个三次就差不多。这三步,差不多用了1个小时。

# Iteration 1
~/opt/biosoft/minimap2/minimap2 ONTmin_IT0.fasta ont.fq > ONTmin_IT0.paf &
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT0.paf ONTmin_IT0.fasta > ONTmin_IT1.fasta &
# Iteration 2
~/opt/biosoft/minimap2/minimap2 ONTmin_IT1.fasta ont.fq > ONTmin_IT1.paf
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT1.paf ONTmin_IT1.fasta> ONTmin_IT2.fasta
# Iteration 3
~/opt/biosoft/minimap2/minimap2 ONTmin_IT2.fasta ont.fq > ONTmin_IT2.paf
time ~/opt/biosoft/racon/build/bin/racon -t 10 ont.fq ONTmin_IT2.paf ONTmin_IT2.fasta > ONTmin_IT3.fasta

之后使用二代数据进行纠错。二代数据虽然短,但是测序质量高,所以一般都要用它进行纠错。推荐用30X PCR free的illuminia 测序数据。

Step 1: 数据预处理,过滤低质量短读,去接头。工具很多,常用的是trimmomatic,cutadapter. 我安利一个国内海普洛斯搞的一个工具fastp。

# data clean
fastp -q 30 -5 -l 100 -i il_1.fq.gz -I il_2.fq.gz -o i1_clean_1.fq -O i1_clean_2.fq 

这里标准为:平均质量高于Q30,对5‘端进行低质量碱基删除,保留大于100bp的短读

Step2: 比对,这一步基本都只用了bwa了

# align
bwa index ONTmin_IT3.fasta
bwa mem -t 8 ONTmin_IT3.fasta il_clean_1.fastq il_clean_2.fastq | samtools sort -@ 8 > ONTmin_IT3.bam

step3: 使用比对后的BAM文件进行纠错

# short read consensus call
java -Xmx16G -jar pilon-1.22.jar --genome ONTmin_IT3.fasta --frags ONTmin_IT3.bam --fix snps --output ONTmin_IT4

二代纠错的时间明显比之前的久,需要一天时间。

大家拿出自己的笔记本实际感受下呗

参考文献

  • nanopore组装拟南芥: High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell
  • 不纠错组装: Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences
  • 三代组装软件评测: Comprehensive evaluation of non-hybrid genome assembly tools for third-generation PacBio long-read sequence data
目录
相关文章
|
8天前
|
数据可视化
数量生态学冗余分析(RDA)分析植物多样性物种数据结果可视化|数据分享
数量生态学冗余分析(RDA)分析植物多样性物种数据结果可视化|数据分享
33 0
|
15天前
基因组组装:Hifiasm 使用教程
基因组组装:Hifiasm 使用教程
26 1
|
7月前
|
存储 JSON Java
GATK4重测序数据怎么分析?
GATK4重测序数据怎么分析?
|
10月前
|
算法 索引 Python
宏基因组之基因组装
宏基因组组装,即把短的reads拼装成连续的序列contig,再根据PE等关系将contig拼装成scaffold。与单个基因组组装不同,宏基因组组装最终得到的是环境样品中全部微生物的混合scaffold。理想情况下一条scaffold对应一个物种的全基因组。但由于序列太短或者覆盖度不够,很难拼出一条完整的基因组。针对高通量测序数据,出现了多种拼接算法和软件。
424 0
|
10月前
|
数据采集 算法
测序质控和基因组组装原理
测序质控和基因组组装原理
|
10月前
|
数据库
利用massdatabase包提取物种KEGG通路与基因/化合物对应信息
最近手头处理一批代谢组数据, 想基于几十个关键差异代谢物代谢物进行下KEGG富集,能想到有两种方式解决,一种常用方式就是基于MetaboAnalyst在线富集,另一种就是解析出该物种的通路与代谢物的对应关系文件,然后用Y叔叔的Clusterprofiler包富集。经一番搜索,massdatabase包可帮我们轻松获得这个文件。 作者:凯凯何_Boy 链接:https://www.jianshu.com/p/654784925903 来源:简书 著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
167 1
|
10月前
|
存储 算法 数据可视化
利用TCseq包进行基因表达趋势分析
TCseq包提供了一个统一的套件去处理不同时序类型的数据分析,可以应用于转录组或者像ATAC-seq,Chip-seq的表观基因组时序型数据分析。该包主要的集中于不同时间点的差异分析,时间趋势分析及可视化作图。
250 0
|
10月前
|
数据采集 设计模式 存储
全基因组重测序流程【超细致!!】
全基因组重测序流程【超细致!!】
|
机器学习/深度学习 算法 数据可视化
|
弹性计算 网络安全 文件存储
BCS三代基因组装
使用阿里云批量计算服务优化的WDL-Canu方案高效、经济地进行三代基因组装。
BCS三代基因组装