细菌基因组拼接工具比较(二)

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了细菌基因组拼接工具比较(二)相关的知识,希望对你有一定的参考价值。

参考技术A

     FastQC 是一款基于Java的软件,它可以快速地对测序数据进行质量评估,并生成网页版的报告。

    每个原始数据文件会生成两个文件,一个为 html网页文件 ,一个为 .zip*文 。打开 .html 可以查看序列的测序质量情况,分为三个等级:合格项为绿色√,警告是黄色的!,不合格为红色的×。(绿色越多越好)

     Cutadapt 是一款用python写的去接头软件,可以去除任意指定的接头和序列,同时也可以用于质量过滤。(实际上原始数据过滤会过滤两种序列:去接头引物、低质量序列)

    给定一组输入,KmerGenie首先计算多个K值的k-mer丰度直方图。然后,对每个K值,它预测数据集中不同基因组k-mer的数量,并返回这个数字最大化的k-mer长度。

    输出文件中主要看 *report.html 网页文件,会推算出最合适的kmer值与基因组大小。在这里,我的基因组得出的kmer=111。

     Velvet ,基于kmer的序列拼接工具,用于基因组的de novo组装,支持多种格式,局部拼接效果好,不善于连接成scaffold。

分为两步执行:

步骤一:利用velveth对数据构建一个hash表

步骤二:velvetg进行序列拼接

     SPAdes 是常用的序列拼接软件之一,支持illumina、PacBio、Nanopore、Sanger、Ion Torrent等测序数据的拼接,同样适合用于混合组装来改善拼接效果(尤其适用于Ion Torrent数据的拼接)

     SOAPdenovo 是华大基因开发的SOAP软件包的一部分,主要用于短序列reads拼接,尤其是illumina测序数据。从小的细菌基因组到大的动植物基因组,人基因组都适用。特点是使用起来比较简单,但是却拥有很好的拼接效果,尤其在基因组构建Scaffold方面,效果很好。

配置文件 config_file 的设置

     ABySS 用于从头拼接双端测序短序列或较大基因组序列。

综合contigs数量、总长度、N50 来看,velvet拼接contig效果最好

    除了Velvet软件之外,其他三款软件在拼接contigs的同时也可以拼接成scaffolds。同样综合scaffolds数量、总长度、N50 来看,从paired-end reads数据一次性拼接成scaffolds效果最好的是ABySS。但是在使用abyss的过程中没有找到调整线程数的参数,用时较长。

    Velvet拼接出的159个contigs已经很接近公司反馈的最终结果了,后面再用报告中提到的两款软件试试。

三代基因组拼接软件--Falcon篇

Falcon是PacBio开发的一款用于三代基因组devono拼接软件,该软件支持PBS,SLURM,SGE,TORQE, LSF作业调度系统。大体的工作原理是将长reads分隔成指定大小模块,这些模块彼此间再相互比对,进行自我纠错以及查找重叠信息,然后再用De Bruijn算法进一步拼接生成contigs。

 

原理简介:

1.Raw sub-reads overlapping for error correction(构建进行错误矫正的重叠信息)

 具体先由fasta2DB与DBsplit对原始reads进行建库分割,分割大小有参数可控,可见下文参数解析,再HPC.daligner生成多个job数的shell脚本放在run_jobs.sh里,接着再根据设定的job_type类型,即根据作业调度系统或者单一节点模式,PBS,SGE通过qsub,SLURM用sbatch将这些job投到相应的计算节点上去,再由每个rj开头的shell脚本对之前分割好的模块进行daligner比对生成相应的.las,再由LAchek,LAsort,LAmerge进行检查,排序,合并。

2.Pre-assembly and error correction (预组装与错误矫正)

 通过LA4Falcon将之前合并好的.las文件转成falcon识别的格式,再通过管道传递给fc_consensus进行预组装与错误矫正,最后生成fasta文件,所生成的fasta文件数应该与0_rawreads/preads/下的cns开头的文件夹数相同,最后的这些preads将被用于下面的组装。

3.overlapping detection of the error corrected reads (错误矫正后的重叠信息群的构建)

     与1相似也是通过fasta2DB与DBsplit对矫正过的reads进行建库分割,也是由HPC.daligner生成一个所有jobs数的集合脚本run_jobs.sh,再根据设定的job_type以及参数投放到节点上,再由daligner,LAcheck,LAsort,LAmerge进行比对,查找,排序,合并。

4.Overlap filtering(重叠信息过滤)

    过虑主要是为了减少计算以及后续拼接的复杂度,过虑的情况主要有:如果其中的一条reads跟一条完全一样就不能提供额外信息,故过滤之;序列两端的重叠信息冗余,不需要太多的重叠信息就可推断出重叠关系的;有些重复区域的会造成覆盖度较高,这些重复区域不能提供额外信息故可通过参数顾虑一些,但在覆盖度偏低的地方,又会造错误率过高,所以也可通过参数过滤掉低覆盖度(fc_ovlp_filter完成的)。

5.Constructing graph from overlap (通过重叠信息构图)

利用fc_ovlp_filter过滤后生成的pread.ovl构建一系列的图信息,边的信息存放在sg_edges_list,然后再进一步连接成unitigs存放在utg_data中,再根据unitigs构建contig图存放在ctg_path中。

6.Constructing contig from graph (通过图信息构建 contigs)

根据图信息以及序列构建contigs。

####################################################################  

三、输入参数解析:主要在fc_run.cfg文件中

[General]

(Falcon支持LOCAL,SLURM,PBS,SGE,TORQUE,LSF等工作类型,需要结合计算机资源以及拼接数据的大小来选择,这里介绍local,slurm与pbs的配置参数)

#job_type=local

#job_type=SLURM

job_type=PBS

input_fofn = input.fofn

#input_fofn = preads.fofn

 

input_type = raw(指定输入数据类型)

#input_type = preads

 

# The length cutoff used for seed reads used for initial mapping

length_cutoff = 10000 用于错误矫正的序列的最低长度

 

#The length cutoff used for seed reads used for pre-assembly

length_cutoff_pr = 10000用于预组装的序列的最低长度

 

#the job of queue

jobqueue = yourqueue

#############################################    

#THE CONFIGURE OF LOCAL

#############################################

#sge_option_da = -pe smp 8 -q %(jobqueue)s【stage-0】

#sge_option_la = -pe smp 2 -q %(jobqueue)s【stage-0】

#sge_option_cns = -pe smp 8 -q %(jobqueue)s【stage-0】

#sge_option_pda = -pe smp 6 -q %(jobqueue)s【stage-1】

#sge_option_pla = -pe smp 2 -q %(jobqueue)s【stage-1】

#sge_option_fc = -pe smp 8 -q %(jobqueue)s [stage-2]

##################################################################

#        THE CONFIGURE OF PBS

##################################################################

sge_option_da = -l nodes=1:ppn=12 -q %(jobqueue)s

sge_option_la = -l nodes=1:ppn=12 -q %(jobqueue)s

sge_option_pda = -l nodes=1:ppn=12 -q %(jobqueue)s

sge_option_pla = -l nodes=1:ppn=12 -q %(jobqueue)s

sge_option_fc = -l nodes=1:ppn=12 -q %(jobqueue)s

sge_option_cns = -l nodes=1:ppn=12 -q %(jobqueue)s

(这是PBS作业调度系统sge_option参数的参考配置信息,具体可以根据集群资源进行设置。)

####################################################################

pa_concurrent_jobs = 20

ovlp_concurrent_jobs = 20

cns_concurrent_jobs = 20

(提交的任务的并发数,根据计算资源来确定)

 

pa_DBsplit_option = -x500 -s200

ovlp_DBsplit_option = -x500 -s200

(-s分割数据块的大小即block的大小,越大需要的内存越大,每个模块间相互比对的时间也越长)

 

pa_HPCdaligner_option = -v -B10 -t16 -e.70 -l1000 -s1000 -T12

-B与jobs总数有关对内存影响不大,e取值范围【0.7,1】raw read设低点。-T线程数

 

ovlp_HPCdaligner_option = -v -B10 -t32 -h60 -e.96 -l500 -s1000 -T12

(参考上面的pa_HPCdaligner_option,e值此时为矫正过的可设高点)

 

falcon_sense_option = --output_muti --min_idt 0.70 --min_cov 4 --max_n_read 200 --n_core 12

(--min_idt 用于装配拼接的序列间比对的最低相似度,--min_cov 种子序列的最低覆盖度,--max_n_read用于错误矫正的最高的read数)

 

overlap_filter_setting = --max_diff 100 --max_cov 100 --min_cov 20 --bestn 12 --n_core 24

(--max_diff序列两端覆盖度的最大差异,--max_cov序列两端的最大覆盖度--minx_cov序列两端的最小覆盖度,--bestn 默认为12,--n_core线程数)

fc_ovlp_to_graph_option = --min_len (length_cutoff_pr,default4000) --min_idt(default 96)

####################################################################

四、运行:

fc_run.py fc_run.cfg

####################################################################

五,参考信息:

https://github.com/PacificBiosciences/FALCON/wiki/Manual

http://blog.sciencenet.cn/blog-3255195-1015620.html 

以上是关于细菌基因组拼接工具比较(二)的主要内容,如果未能解决你的问题,请参考以下文章

如何在genbank查找某个细菌的基因序列?

CRISPR/Cas9|InParanoid|orthoMCL|PanOCT|pan genome|meta genome|Core gene|CVTree3|

易基因:全基因组ChIP-seq分析揭示细菌转录因子PhoB的基因内结合位点|mBio

基因变异(codevs 3194)

细菌多位点序列分型(Multilocus sequence typing,MLST)的原理及分型方法

练习题 | 细菌繁殖