1. 测序前的准备 搜集物种相关信息,比如基因组大小,杂合度, 1.1 获取基因组大小 基因组大小的获取关系到对以后组装结果的大小的正确与否判断;基因组太大(10Gb),超出了目前denovo组装基因组软件的对机器内存的要求,从客观条件上讲是无法实现组装的。 一般物种的基因组大小可以从( http://www.genomesize.com/ )这个数据库查到。如果没有搜录,需要考虑通过实验(流式细胞仪)获得基因组大小。 1.1.1 流式细胞仪估计基因组大小的例子: Yoshida, S., J. K. Ishida, et al. (2010). A fulllength enriched cDNA library and expressed sequence tag analysis of the parasitic weed, Striga hermonthica ant Biol 10: 55. 1.1.2 基于福尔根染色估计基因组大小的描述: 这本书比较经典,重点推荐:Gregory, T. (2005). The evolution of the genome, Academic Press. 1.1.3 定量pcr估计基因组大小的例子: Wilhelm, J., A. Pingoud, et al. (2003). Real-time PCR based method for the estimation of genome sizes. Nucleic Acids Res 31(10): e56. Jeyaprakash, A. and M. A. Hoy (2009). The nuclear genome of the phytoseiid Metaseiulus occidentalis (Acari: Phytoseiidae) is among the smallest known in arthropods. Exp Appl Acarol 47(4): 263-273. 1.1.4 Kmer估计基因组大小的例子: Kim, E. B., X. Fang, et al. (2011). Genome sequencing reveals insights into physiology and longevity of the naked mole rat. Nature 479(7372): 223-227. 1.2 杂合度估计 杂合度对基因组组装的影响主要体现在不能合并姊妹染色体,杂合度高的区域,会把两条姊妹染色单体都组装出来,从而造成组装的基因组偏大于实际的基因组大小。 一般是通过SSR在测序亲本的子代中检查SSR的多态性。杂合度如果高于0.5%,则认为组装有一定难度。杂合度高于1%(?)则很难组装出来。 杂和度估计一般通过kmer分析来做,这里有一个例子: http://www.nature.com/nature/journal/vaop/ncurrent/full/nature11413.html 降低杂合度可以通过很多代近交来实现。 杂合度高,并不是说组装不出来,而是说,装出来的序列不适用于后续的生物学分析。比如拷贝数、基因完整结构。 1.3 是否有遗传图谱可用 随着测序对质量要求越来越高和相关技术的逐渐成熟,遗传图谱也快成了denovo基因组的必须组成。构建 遗传图构建相关概念可以参考这本书(The handbook of plant genome mapping: genetic and physical mapping ) 1.4 生物学问题的调研 这一步也是很重要的 2. 测序样品准备 确定第一步没问题,就意味着这个物种是可以尝试测序的。测序样品对一些物种也是很大问题的,某些物种取样本身就是一个挑战的问题。 基因组测序 用的样品最好是来自于同一个个体,这样可以降低个体间的杂和对组装的影响。大片段对此无要求。 3. 测序策略的选择 一般都是用不同梯度的插入片段来测序,小片段(200,500,800)和大片段(1k, 2kb 5kb 10kb 20kb 40kb)。如果是杂合度高和重复序列较多的物种,可能要采取fosmid-by-fosmid或者fosmid pooling的策略。 不言而喻,后者花费是相当高的。 4. 基因组组装 4.1 组装相关综述 : Li, Z., Y. Chen, et al. (2012). Comparison of the two major classes of assembly algorithms: overlap-layout-consensus and de-bruijn-graph. Brief Funct Genomics 11(1): 25-37. Treangen, T. J. and S. L. Salzberg (2012). Repetitive DNA and nextgeneration sequencing: computational challenges and solutions. Nat Rev Genet 13(1): 36-46. http://www.cbcb.umd.edu/research/assembly_primer.shtml Schatz, M. C., J. Witkowski, et al. (2012). Current challenges in de novo plant genome sequencing and assembly. Genome Biol 13(4): 243 Baker, M. (2012). De novo genome assembly: what every biologist should know . Nat Methods 9(4): 333-337. (重点推荐) Compeau, P. E., et al. (2011). How to apply de Bruijn graphs to genome assembly. Nat Biotechnol 29(11): 987-991. Birney, E. (2011). Assemblies: the good, the bad, the ugly. Nat Methods 8(1): 59-60. Schatz, M. C., et al. (2010). Assembly of large genomes using second-generation sequencing. Genome Res 20(9): 1165-1173. 4.2 测序数据预处理 Illumina数据预处理终结者 数据纠错 Brief Bioinform 测序数据错误校正方法测评 Yang, X., et al. (2013). A survey of error-correction methods for next-generation sequencing. Brief Bioinform 14(1): 56-66. Kelley, D. R., M. C. Schatz, et al. (2010). Quake: quality-aware detection and correction of sequencing errors. Genome Biol 11(11): R116. 4.3 组装软件比较 Vezzi, F., et al. (2012). Reevaluating assembly evaluations with feature response curves: GAGE and assemblathons. PLoS One 7(12): e52210. Salzberg, S. L., A. M. Phillippy, et al. (2012). GAGE: A critical evaluation of genome assemblies and assembly algorithms. Genome Res 22(3): 557-567. Zhang, W., et al. (2011). A practical comparison of de novo genome assembly software tools for nextgeneration sequencing technologies. PLoS One 6(3): e17915. Narzisi, G. and B. Mishra (2011). Comparing de novo genome assembly: the long and short of it. PLoS One 6(4): e19175. Lin, Y., et al. (2011). Comparative Studies of de novo Assembly Tools for Nextgeneration Sequencing Technologies. Bioinformatics. Hayden, E. C. (2011). Genome builders face the competition. Nature 471(7339): 425. Finotello, F., et al. (2011). Comparative analysis of algorithms for whole-genome assembly of pyrosequencing data. Brief Bioinform. Earl, D. A., et al. (2011). Assemblathon 1: A competitive assessment of de novo short read assembly methods. Genome Res. 4.4 组装质量评估 Schatz, M. C., et al. (2011). Hawkeye and AMOS: visualizing and assessing the quality of genome assemblies. Brief Bioinform. RibaGrognuz, O., et al. (2011). Visualization and quality assessment of de novo genome assemblies. Bioinformatics. 个人见解: 目前大基因组的denovo组装主流软件还是ALLPATH-LG, SOAPdenovo ALLPATH-LG的优点是:组装的连续性最好,准确性最好,但是消耗内存较大,不是太好使用 SOAPdenovo的优点是:速度快,消耗的内存可以接受,组装的连续性还可以,但是错误相对要多一些。 当然,上述评述并不是在所有情况下的,对不同物种,不同数据,他们的表现可能会不一样。 基于Overlap-layout的方法的组装软件首推CABOG,这是当年用来组装果蝇基因组的原型。另外,快要发布的MSR-CA貌似也不错,其整合了上述所有软件的优点,来势很猛啊。 5. 基因组注释 Yandell, M. and D. Ence (2012). A beginner's guide to eukaryotic genome annotation. Nat Rev Genet 13(5): 329-342. 6. 基因组可视化 Nielsen, C. B., M. Cantor, et al. (2010). Visualizing genomes: techniques and challenges. Nat Methods 7(3 Suppl): S5-S15. 7. 进化分析 Yang, Z. and B. Rannala (2012). Molecular phylogenetics: principles and practice. Nat Rev Genet 13(5): 303-314. 8. 经典案例 Colbourne, J. K., M. E. Pfrender, et al. (2011). The ecoresponsive genome of Daphnia pulex. Science 331(6017): 555-561. Kim, E. B., X. Fang, et al. (2011). Genome sequencing reveals insights into physiology and longevity of the naked mole rat. Nature 479(7372): 223-227. Grbic, M., T. Van Leeuwen, et al. (2011). The genome of Tetranychus urticae reveals herbivorous pest adaptations. Nature 479(7374): 487-492. 转自: http://seq.cn/forum.php?mod=viewthreadtid=4607reltid=17953pre_thread_id=0pre_pos=8ext =
利用 k-mer 组装的基因组的软件现在已经很多了,例如soapdenovo, velvet 等等。 PLoB 上已经有不少关于 velvet 的软件的介绍,今天再次谈谈 velvet 这个软件,主要是推荐一些 velvet 配套的其他软件。 1、 VAGUE 这是一个基于 velvet 的但是带有图形化界面的基因组组装软件。目前该软件支持linux和mac。关于 VAGUE 的介绍可以读读下面的英文信息。 VAGUE is a vague acronym for "Velvet Assembler Graphical Front End", which means it is a GUI for the Velvet de novo assembler. The command line version of Velvet can be complicated for beginners to use, but VAGUE makes it clear and simple. 上面是软件的使用截图,更多软件的信息请访问:http://www.vicbioinformatics.co/software.vague.shtml 2、 VelvetK 利用基于K-mer组装的软件的时候我们常常关心一个问题,K-mer值设置多少比较好呢?于是乎,我们常做的一件事情就是把很多个K-mer值都拿来试一下,看看哪个K-mer比较好。这不失为一种好办法,但是当数据量很大的时候确实有点浪费时间。 我记得之前在 velvet 的manual里面作者推荐了一个经验公式(如下文),该公式推荐的K-mer与reads覆盖度和reads长度相关。 All coverage values in Velvet are provided in k-mer coverage, i.e. how many times has a k-mer been seen among the reads. The relation between k-mer coverage C k and standard (nucleotide-wise) coverage C is C k = C(L−k +1)/L where k is your hash length, and L you read length 这个经验公式也是一个定性的公式了。今天这里我们给大家推荐的另外一个软件 VelvetK ,这个软件可以告诉你 k-mer 选择多少比较合适。 下面是作者对 VelvetK 的描述介绍: VelvetK can estimate the best k-mer size to use for your Velvet de novo assembly. It needs two inputs: the estimated genome size, and all your sequence read files. The genome size can be supplied as as a number (eg. 3.5M) or as a FASTA file of a closely related genome. The reads can be FASTA or FASTQ, and may optionally be compressed with GZIP or BZIP2. 了解更多信息可以从这里进入:http://www.vicbioinformatics.com/software. velvet k.shtml 下面是 velvet K的使用 参数 : 3、 VelvetOptimiser 上面介绍了在使用 velvet 过程中K-mer的 参数 如何设置,这只是一部分。下面给大家推荐另一款软件, VelvetOptimiser 。这款软件的作用是,帮助你优化 velvet 的 参数 ,针对你的数据,提供一个合理的 参数 ,这些 参数 包括,K, -exp_cov, -cov_cutoff。 下面是该软件的描述: VelvetOptimiser is a multithreaded Perl script for automatically optimising the three primary parameter options (K, -exp_cov, -cov_cutoff) for the Velvet de novo sequence assembler. 想知道更多关于 VelvetOptimiser 的资料和介绍,请从这里进入: http://www.vicbioinformatics.com/software. velvet optimiser.shtml 转自:http://www.plob.org/2012/11/21/4797.html
1. kmer值 kmer值哪个比较好很难说,对不同的数据,用不同的kmer值会有很不同的结果。最好的办法就是测试不同的kmer,然后看结果的N50,找到N50最 高的kmer。不过SOAPdenovo最新的版本已经支持最长127bp的kmer了,所以要从20多测试到127,显然不太可能。下面是文档中对 kmer的说明。 How to set K-mer size? The program accepts odd numbers between 13 and 127. Larger K-mers would have higher rate of uniqueness in the genome and would make the graph simpler, but it requires deep sequencing depth and longer read length to guarantee the overlap at any genomic location. 根据我的实际使用经验,如果你的read足够长,覆盖度足够高,kmer设的越高越好。 但是实际情况是,测序的覆盖度经常不够,或者用早期的GA平台测出来read长度只有35bp,或者为了节省成本,在mate-pair library(长片段insert的文库,一般2kb)测序时双端只有70bp,甚至40bp之类的,情况比较复杂。 一般来说,我尽量使用更高的kmer,如果我有100bp的pair-end,50bp的mate-pair,而且覆盖度挺高,我就用到kmer=45左 右,如果mate-pair只有40bp,kmer=35左右。如果mate-pair更短,只有35bp,kmer值就再降一点。 但是覆盖度不够时,我一般还是使用kmer=25来拼接。 SOAP推出了一个新的工具,prepare模块,似乎就是为了解决混合长度read的问题,你可以先用很高的kmer进行contig的拼接,只使用来 自180bp,300bp,500bp双端100bp的pair-end文库的reads。之后使用这些contig进行scaffolding,你可以 为这些contig重新构建一个短的kmer graph,然后整合来自mate-pair文库的短read进行scaffolding.只是遗憾的是,这个模块我尝试过,不行。这个模块的运行 SOAP应该要更详细说明才行。 7) Data Preparation Module generates necessary data for SOAPdenovo to run "map" and "scaff" steps from Contigs generated by SOAPdenovo or other assemblers with various length of kmer. options: -g Prefix of output. -K Kmer length. -c Input Contigs FASTA. (Filename cannot be prefix.contig) 2. config文件中的一个重要参数reverse_seq 这个参数在很多使用soapdenovo进行拼接的人当中都会设置错误,因为默认是0,下面是从其他人博客中看到的: “reverse_seq This option takes value 0 or 1. It tells the assembler if the read sequences need to be complementarily reversed. Illumima GA produces two types of paired-end libraries: a) forward-reverse, generated from fragmented DNA ends with typical insert size less than 500 bp; b) forward-forward, generated from circularizing libraries with typical insert size greater than 2 Kb. The parameter “reverse_seq” should be set to indicate this: 0, forward-reverse; 1, forward-forward. RF: first read of fragment pair is sequenced as anti-sense (reverse), and second read is in the sense strand (forward); FR: first read is in the sense strand (forward);second read of fragment pair is sequenced as anti-sense (reverse)” “read/1,read/2哪个是正向?哪个是反向的? ,这个是不能确定的,在华大这边建库小插入片段(2000bp)是打断后直接建库,测两端,没有反转;(=2000bp)的文库是打断后环化,环化后再打断测,这时称为reverse_seq,在soapdenovo里面reverse_seq=1 转录组的不用管这个了,都是小于2k的” 根据我们的实际经验,如果你有mate-pair的文库,那你得咨询建库的人,中间是否经过了环化这个步骤,如果有,则必须把revser_seq设置为 1。之前有个孩子不看参数,全部默认0,拼出来的N50只有30k,惨不忍睹,把这个参数纠正后,N50就有150k了。 其实这个参数就算设置正确,还是存在一个很严重的问题, 因 为 在illumina mate-pair library建库的过程中,由于实验方法上的技术缺陷,很多dna片断并没有被成功的环化,这些没有被环化的片断测序后实际上还是pair end,也就是说两个read的方向是FR,而非RF,这时你用了reverse_seq=1,这些read的方向就是错的。所以,如果有可能,特别是你 已经有reference的时候,尽量先对mate-pair文库的read进行筛选,把那些insert-size远小于理论值的,方向不正确的 read删掉,否则你的拼接将会引入大量的错误。 3.其他几个对N50有重要影响的参数 -M mergeLevel(default 1,min 0, max 3): the strength of merging similar sequences during contiging 这个参数似乎在上一篇博文中没有解释,并且SOAP的提示也很简单。默认情况下M=1,M的值可以设置,从0到3。这个参数的作用是在拼contig的过程中,对buble合并和分解的一个重要参数。详细可见下面一段文字: We used Dijkstra’s algorithm to detect bubbles, which is similar to the ‘‘Tour-bus’’ method in Velvet.We merged the detected bubbles into a single path if the sequences of the parallel paths were very similar; that is, only had a single base pair difference or had fewer than four base pairs difference with 90% identity. 简单说,就是一些差别很小的kmer,例如只有1个mismatch,在M值高的时候,将会被合并到一个node里面。 这个对于杂合度高的基因组意义比较大,因为两个allel在拼接时可能会被独立的拼出来,如果你把M值调高,这些allel就可以被合并在一起。 但是如果一个基因组repeat很多的话,M值设高了就会把很多差别很小的repeat序列合并在一起,那么相当多的序列将无法被用于构建edges,这 时你的N50就会比较差,降低这个值可以显著的提高N50,但是吧,由于测序错误的存在,或者repeat特别多时候,降低M值又会导致误拼的概率大大提 升,这个问题很难说得清,有些人在拼接时候尽量提高M值,把尽可能多的repeat与及可能的测序错误都合并到一个区域,这样可以保证基因组其他区域拼接 的准确率,但是代价是N50的显著降低。 一般情况下,我还是使用了默认的M=1。要是碰上repeat多,杂合度又高的基因组,我只能建议远离illumina,远离SOAPdenovo,呵呵。 -R 也是跟repeat分解有关的参数,用他可以把一些由短片段重复而被合并在一起node分解开,一般情况下也可以明显的提高N50,但是一样的,repeat很多的情况下,这个参数要慎用,误拼的概率会大大提升。 config文件里的asm_flag参数 1 (reads only used for contig assembly), 2 (only used for scaffold assembly) and 3 (used for both contig and scaffold assembly). 如果你的pair-end数据不够,那就让mate-pair文库的序列也用于contig拼接,asm_flag=3。也可以比较明显的提高N50。 4.提醒 以上这些都是我的个人经验,我相信大部分情况下还是适用的,但是每个基因组的情况的差别很大,所以一定还是要多测试,才能得到比较好的拼接结果。短 read拼接最大的问题就是处理repeat序列,N50提高了,某种程度上拼接准确率也要下降。因为repeat被放到哪个scaffold,似乎在 SOAP中没有特别优化,只是随机扔个位置,这样就会出现不该连的scaffold被连在一起,甚至在contig内部都会出现错误,于是人为造成了很多 假的structural variation。不过一般来说,insertion,deletetion之类的错误还可以忍受,但是还是有不少 inversion,translocation,这种就是比较严重的错误了
我们都知道,测序本身并不难,难就难在基因组的后续组装拼接,因为它涉及到大量需要考虑的问题,如重复、到位、覆盖率等等,于是如何有效的得到最后的序列或者有意义的Scaffold是做基因组面临的一个很大问题。不同的人去做会得到不同的结果,如N50、N90,scaffold数量等等。 下面简单介绍一下SOAPdenovo组装的一般过程: Schematic overview of the assembly algorithm. (A) Genomic DNA was fragmented randomly and sequenced using paired-end technology.Short clones with sizes between 150 and 500 bp were amplifiedand sequenced directly; while long range (2–10 kb) paired-end libraries were constructed by circularizing DNA, fragmentation, and then purifying fragments with sizes in the range of 400–600 bp for cluster formation. (B) The raw or precorrected reads were then loaded into computer memory and de Bruijn graph data structure was used to represent the overlap among the reads. (C) The graph was simplified by removing erroneous connections (in red color on the graph) and solving tiny repeats by readpath: (i) Clipping the short tips, (ii) removing low-coverage links,(iii) solving tiny repeats by read path, and (iv) merging the bubbles thatwere caused by repeats or heterozygotes of diploid chromosomes. (D) On the simplified graph, we broke the connections at repeat boundaries and output the unambiguous sequence fragments as contigs. (E)We realigned the reads onto the contigs and used the paired-end information to join the unique contigs into scaffolds. (F) Finally, we filled in the intrascaffold gaps,which were most likely comprised by repeats, using the paired-end extracted reads. 更多内容 http://www.seq.cn/forum.php?mod=viewthreadtid=3451
转自; http://www.jb51.net/article/18092.htm @echo off :: 把多行文本拼接成用;连接的一行 :: nul 不能省略,省略掉就无法运行下去 for /f "tokens=*" %%i in (源文件.txt) do set /p "var=%%i;" nul 目标文件.txt exit 另一种方法: set tmpstr= setlocal enabledelayedexpansion for /f "tokens=*" %%i in (1.txt) do set tmpstr=!tmpstr! %%i echo %tmpstr% pause exit 还有一种不启用变量延迟的方法,能兼容除英文双引号外的所有特殊字符: @echo off :: code by jm 2006-12-14 for /f "delims=" %%i in (1.txt) do call set "var=%%var%%%%i" echo "%var%" pause