PacBio sequence error correction amd assemble via pacBioToCA

Illumina二代测序有个致命缺陷,说到底还是基于PCR扩增的,所以存在偏向性和对于高GC含量区无法扩增等系统误差,测序错误是不可避免的,其次就是测序长度短;但其价格便宜,通量非常高,准确性达99%,综合性价比也受到青睐。短序列的reads在做基因组装的时候,遇到大的重复片段就会很吃力。

10X Genomics

2015年备受瞩目的测序黑马:10X Genomics,是常规Illumina二代测序的升级版,由于开发出了一套巧妙的Barcoding建库方案,使得Illumina这种短读长二代测序能够得到跨度在30-100Kb的linked reads信息,与二代测序数据相结合,在Scaffold的组装上能够得到媲美三代测序的组装结果;


基本原理: 首先将每一条长片段的DNA分配至不同的油滴微粒中,通过专利的GEM建库技术,长片段DNA被切碎成适合测序的大小,并且来源于相同油滴(同一条长片段DNA)的DNA片段,会带上相同的一段DNA序列标记(Barcode),之后在Illumina系统上测序完成后,可以理论上再将来源相同的DNA序列独立拼接,得到原先的长片段DNA序列。
对于不同GC含量区其效果如何呢?2015年10月Nat Review Genetics文章Genetic variation and the de novo assembly of human genomes中总结的PacBio、10X Genomics以及Illumina技术在不同GC含量DNA区域的覆盖度分布:

10X Genomics技术相对于Illumina来说,有改进,但依旧是个拱形,而PacBio则是无偏倚的均一分布,10X的技术,其Coverage一样是受GC含量影响较大的,那么如果真要应用10X技术,那么必须注意目标DNA的GC含量分布最好能控制在30~70%。
但10Xgenome毕竟是升级版,其也存在一些特有的优势:

  • (1) 微量样本:仅需1ng基因组DNA即可进行长片段建库;
  • (2) 精确分区:由于拥有众多的barcode和Partions,可对DNA进行精确分区;
  • (3) 长片段信息:该技术可与Illumina测序仪进行无缝对接,利用短Reads可获得长达100Kb的片段;
  • (4) 基因组组装质量提升:利用长片段信息结合Illumina组装数据组装的ScaffoldN50长度比单纯用Illumina方法提高十几倍。
  • PacBio



    第三代测序中的PacBio单分子实时(Single Molecule Real-Time, SMRT)DNA测序可以实现超过99.999%(QV50)的高度精确测序,且不受DNA序列中GC和AT含量的影响,平均读长可达20kb(最长>60kb)。

    PacBio三代测序最大的死穴是:通量不足和单次(1X)测序错误率高(85%);但三代的错误是完全随机发生的,属于随机误差,可以靠覆盖度来自我纠错,如果通量不是限制因素,那么PacBio是目前最准确的测序方式:错误率可以无限接近罕见突变的发生率(即无法分辨是测序错误还是罕见突变)。2012年冷泉港实验室的Michael Schatz开发了一种纠错算法,用二代测序的短读长高精确数据对三代长读长数据进行纠错,这种称为”混合纠错拼接” (Hybrid error correction and de novo assembly of single-molecule sequencing reads)可以进一步提升PacBio测序精确度。

    PBcR: 混合纠错拼接


    PBcR: 混合纠错拼接 粉色长方形:单个PacBio RS reads;黑色竖线:测序错误;(a)由于测序错误碱基的存在使得两条reads就难确定是否在末端重叠;(b)高质量的短reads比对到存在错误的长reads;短reads中的黑色竖线表示 ‘mapping errors’ ,是长reads和短reads中测序错误的组合,此外双拷贝的重复序列的存在(灰色轮廓)导致在每一个拷贝中出现短reads的堆挤,为避免reads map到错误的重复区,仅保留最高比对值的短reads;(c)剩余的比对形成一致性序列(紫色长方形),长reads和短reads中共有的部分错误未能得到纠正;(d)overlap纠正后的长reads;(e) 最后的组装能够跨越重复区域。

    Illumina reads纠错覆盖度




    纠错的准确性和组装一致性在Illumina高质量reads达50X后开始收益递减,因此50X Illumina reads足够,纠错后PacBio长reads准确性将由85%提升至>99.9%,此时嵌合体和错误剪切reads分别为<2.5% 和 <1%。
    目前在P6C4试剂下,大约每SMRT Cell平均可以做到 600M~1G数据量。
    PacBio的长读长、无GC偏向性和无PCR扩增偏向性等独特优势有助于克服复杂的重复区域,从而跨越整个基因转录区,显著提升基因组和转录组的De Nove组装质量;

    Illumina二代+PacBio三代数据分析

    PBcR首先通过纠错来提升PacBio reads准确性,然后进行组装。PBcR的纠错和组装分为self-correction (using only PacBio RS data,自动运行fastqToCA) or correction with high-identity sequences(二代数据)。

    self-correction

    1
    PBcR -length 500 -partitions 200 -l lambda -s pacbio.spec -fastq pacbio.filtered_subreads.fastq genomeSize=50000 > run.out 2>&1

    高质量Illumina reads

    1
    2
    3
    4
    5
    6
    7
    #short read准备
    fastqToCA -libraryname illumina -technology illumina -reads illumina.fastq > illumina.frg
    #纠正
    pacBioToCA -length 500 -partitions 200 -l ec_pacbio -t 16 -s pacbio.spec \
    -fastq pacbio.filtered_subreads.fastq illumina.frg > run.out 2>&1
    #组装
    runCA -p asm -d asm -s asm.spec ec_pacbio.frg > asm.out 2>&1

    : 第一步short reads准备阶段请确认二代数据第四行质量编码值,一般是33,否则用-type参数指定,要不然会报错QV问题;
    纠正时PBcR需要安装AMOS和blasr依赖软件,输入文件short reads (illumina.frg)和long reads (pacbio.filtered_subreads.fastq);
    fastqToCA和PBcR两个中的libraryname需不同;
    fastqToCA生成的frg文件后面没有序列信息 ,是正确的;

    Spec files参数解释

    PBcR混合组装需要指定两个Spec配置文件: pacbio.spec(纠错)和asm.spec(组装)。这两个文件都包含特定的算法参数和计算机硬件参数,通常情况下算法参数可以忽略(此时将用软件默认值),但是计算机硬件参数需要根据实际情况调整。
    所有参数均为option = value形式,其中的value为布尔型(boolean),即true=1,false=0。
    具体关于specfile参数解释见PBcR:SpecFiles Options

    Spec files实例

    集群下参考pacbio.spec

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    #以下为grid计算
    stopAfter=overlapper

    # original asm settings
    utgErrorRate = 0.25
    utgErrorLimit = 4.5

    cnsErrorRate = 0.25
    cgwErrorRate = 0.25
    ovlErrorRate = 0.25

    merSize=14

    merylMemory = 128000
    merylThreads = 16

    ovlStoreMemory = 8192

    # grid info
    useGrid = 1
    scriptOnGrid = 1
    frgCorrOnGrid = 1
    ovlCorrOnGrid = 1

    sge = -V -S /bin/sh
    #sge = -V -A assembly
    sgeScript = -pe smp 16
    sgeConsensus = -pe smp 1
    sgeOverlap = -pe smp 4
    sgeFragmentCorrection = -pe smp 2
    sgeOverlapCorrection = -pe smp 1

    #ovlMemory=8GB --hashload 0.7
    ovlHashBits = 25
    ovlThreads = 4
    ovlHashBlockLength = 20000000
    ovlRefBlockSize = 50000000

    # for mer overlapper
    merCompression = 1
    merOverlapperSeedBatchSize = 500000
    merOverlapperExtendBatchSize = 250000

    frgCorrThreads = 2
    frgCorrBatchSize = 100000

    ovlCorrBatchSize = 100000

    ######################以下为非gird计算,useGrid = 0 ####################
    # non-Grid settings, if you set useGrid to 0 above these will be used
    merylMemory = 128000
    merylThreads = 4

    ovlStoreMemory = 8192

    ovlConcurrency = 6

    cnsConcurrency = 16

    merOverlapperThreads = 2
    merOverlapperSeedConcurrency = 6
    merOverlapperExtendConcurrency = 6

    frgCorrConcurrency = 8
    ovlCorrConcurrency = 16
    cnsConcurrency = 16

    集群下参考asm.spec

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    ######################以下为gird计算,useGrid = 1 ####################
    cnsErrorRate = 0.10
    ovlErrorRate = 0.10

    overlapper = ovl
    unitigger = bogart
    utgBubblePopping = 1

    merSize = 14

    merylMemory = 128000
    merylThreads = 16

    ovlStoreMemory = 8192

    # grid info
    useGrid = 1
    scriptOnGrid = 1
    frgCorrOnGrid = 1
    ovlCorrOnGrid = 1

    sge = -V -S /bin/sh
    sgeScript = -pe smp 16
    sgeConsensus = -pe smp 1
    sgeOverlap = -pe smp 4
    sgeFragmentCorrection = -pe smp 2
    sgeOverlapCorrection = -pe smp 1

    #ovlMemory=8GB --hashload 0.7
    ovlHashBits = 25
    ovlThreads = 6
    ovlHashBlockLength = 20000000
    ovlRefBlockSize = 5000000

    # for mer overlapper
    merCompression = 1
    merOverlapperSeedBatchSize = 500000
    merOverlapperExtendBatchSize = 250000

    frgCorrThreads = 2
    frgCorrBatchSize = 100000

    ovlCorrBatchSize = 100000

    ######################以下为非gird计算,useGrid = 0 ####################
    # non-Grid settings, if you set useGrid to 0 above these will be used
    merylMemory = 128000
    merylThreads = 12

    ovlStoreMemory = 8192

    ovlConcurrency = 8

    merOverlapperThreads = 6
    merOverlapperSeedConcurrency = 2
    merOverlapperExtendConcurrency = 2

    frgCorrConcurrency = 8

    ovlCorrConcurrency = 16
    cnsConcurrency = 16

    doToggle=0
    toggleNumInstances = 0
    toggleUnitigLength = 2000

    doOverlapBasedTrimming = 1
    doExtendClearRanges = 2

    输出结果

  • 最后组装结果文件夹9-terminator;
  • 主要输出文件是prefix.asm,以分层的数据结构提供组装结果的精确描述,包含生成的contig 和 scaffold 序列;
  • prefix.qc,关于组装结果的统计信息;
  • 参考资料

    当10X Genomics遇上PacBio——烫金开始剥落了
    第三代测序成本偏高是什么原因导致的?
    Computational Science Community Wiki: Sun Grid Engine: Job Arrays

    tiramisutes wechat
    欢迎关注