1.PBJelly2 简介

PBJelly2 用于利用 Pacbio 数据进行基因组补洞和 scaffold 连接。

2.安装 PBJelly2

安装 HDF：
$ wget http://www.hdfgroup.org/ftp/HDF5/current/bin/linux-centos6-x86_64/hdf5-1.8.15-patch1-linux-centos6-x86_64-shared.tar.gz
$ tar zxf hdf5-1.8.15-patch1-linux-centos6-x86_64-shared.tar.gz
$ mv hdf5-1.8.15-patch1-linux-centos6-x86_64-shared /opt/biosoft/hdf5-1.8.15-patch1
$ export HDF5INCLUDEDIR=/opt/biosoft/hdf5-1.8.15-patch1/include/
$ export HDF5LIBDIR=/opt/biosoft/hdf5-1.8.15-patch1/lib/
$ export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/biosoft/hdf5-1.8.15-patch1/lib/
$ export C_INCLUDE_PATH=$C_INCLUDE_PATH:/opt/biosoft/hdf5-1.8.15-patch1/include/

安装 Blasr（需求 HDF v1.8.0或以上版本）：
$ wget https://github.com/PacificBiosciences/blasr/archive/master.zip -O blasr.zip
$ unzip blasr.zip
$ mv blasr-master/ /opt/biosoft/blasr
$ cd /opt/biosoft/blasr
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/blasr/alignment/bin/' >> ~/.bashrc

安装Python模块 Networkx v1.7
$ wget https://pypi.python.org/packages/source/n/networkx/networkx-1.7.tar.gz
$ tar zxf networkx-1.7.tar.gz
$ cd networkx-1.7
$ sudo /usr/local/bin/python setup.py install

安装Python模块 pyparsing
$ wget https://pypi.python.org/packages/source/p/pyparsing/pyparsing-2.0.3.tar.gz
$ tar zxf pyparsing-2.0.3.tar.gz
$ cd pyparsing-2.0.3
$ sudo /usr/local/bin/python setup.py install

安装Python模块 numpy
$ wget https://pypi.python.org/packages/source/n/numpy/numpy-1.9.2.tar.gz
$ tar zxf numpy-1.9.2.tar.gz 
$ cd numpy-1.9.2
$ sudo /usr/local/bin/python setup.py install

安装Python模块 h5py
$ wget https://pypi.python.org/packages/source/h/h5py/h5py-2.5.0.tar.gz
$ cd h5py-2.5.0
$ export LIBRARY_PATH=/opt/biosoft/hdf5-1.8.15-patch1/lib/:$LIBRARY_PATH
$ /usr/local/bin/python setup.py build
$ sudo /usr/local/bin/python setup.py install

安装Python模块 pysam
$ https://pypi.python.org/packages/source/p/pysam/pysam-0.8.3.tar.gz
$ tar zxf pysam-0.8.3.tar.gz 
$ cd pysam-0.8.3
$ sudo /usr/local/bin/python setup.py install

安装Python模块 intervaltree
$ wget https://pypi.python.org/packages/source/i/intervaltree/intervaltree-2.0.4.tar.gz
$ tar zxf intervaltree-2.0.4.tar.gz
$ cd intervaltree-2.0.4
$ sudo /usr/local/bin/python setup.py install

安装 PBJelly2
$ wget http://sourceforge.net/projects/pb-jelly/files/latest/download -O PBSuite.tar.gz
$ tar zxf PBSuite.tar.gz -C /opt/biosoft/
$ SWEETPATH=`ls /opt/biosoft/PBSuite* -d`
$ echo "perl -p -i -e 's#export SWEETPATH=.*#export SWEETPATH=$SWEETPATH#' $SWEETPATH/setup.sh" | sh
$ echo "source $SWEETPATH/setup.sh" >> ~/.bashrc

3. 使用 PBJelly2 进行补洞

首先创建配置文件 Protocol.xml，内容如下:

<jellyProtocol>
    <reference>基因组fasta文件的路径</reference>  
    <outputDir>输出文件路径</outputDir>
    <blasr>-minMatch 8 -minPctIdentity 70 -bestn 1 -nCandidates 20 -maxScore -500 -nproc 24 -noSplitSubreads</blasr>
    <input baseDir="输入Pacbio数据文件所在的文件夹">
        <job>Pacbio数据文件名称</job>
    </input>
</jellyProtocol>

然后依次运行下6步：

$ Jelly.py setup Protocol.xml
$ Jelly.py mapping Protocol.xml
$ Jelly.py support Protocol.xml
$ Jelly.py extraction Protocol.xml
$ Jelly.py assembly Protocol.xml -x "--nproc=24"
$ Jelly.py output Protocol.xml

--nproc 参数设置运行线程数。

输出结果文件为 jelly.out.fasta 。

4. 使用 PBJelly2 进行 scaffold 连接

1. PacBioToCA 简介

PacBioToCA 是 Celera Assembler 的一个模块，用于Pacbio数据的校正。因此，使用 PacBioToCA 需要安装 Celera Assembler。
利用 PacBioToCA 进行校正时：如果 Pacbio 数据量 > 20x 时，可以进行自我校正，而不需要其它更为精确的数据；推荐 Pacbio 数据量 > 50x 时进行自我校正；通常 > 15x 的 Pacbio 数据即可得到很好的组装结果。

2. 安装 Celera Assembler

推荐使用源码安装，直接下载二进制版本，可能提示错误：/lib64/libc.so.6: version `GLIBC_2.14′ not found

$ wget http://sourceforge.net/projects/wgs-assembler/files/wgs-assembler/wgs-8.3/wgs-8.3rc2.tar.bz2
$ tar jxf wgs-8.3rc2.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/wgs-8.3rc2/kmer
以下使用 make，不要使用 -j 参数进行多线程编译，会出错。
$ make && make install && cd ..
$ cd src && make && cd ..
$ echo 'PATH=$PATH:/opt/biosoft/wgs-8.3rc2/Linux-amd64/bin/' >> ~/.bashrc

如果进行 Pacbio 数据的自我校正，则需要 JAVA 1.8 版本
$ wget http://javadl.sun.com/webapps/download/AutoDL?BundleId=106240 -O jre1.8.0_45.tar.gz
$ tar zxf jre1.8.0_45.tar.gz -C /opt/biosoft
$ echo 'PATH=/opt/biosoft/jre1.8.0_45/bin/:$PATH' >> ~/.bashrc
$ source ~/.bashrc

3. PacBioToCA 的使用

3.1 将错误率低的测序数据转换为 FRG 数据

Celera Assembler 软件输入的测序数据的文件格式是 FRG 格式。PacBioToCA 输入的错误率低的数据需要 FRG 格式。
1. 使用 fastqToCA 将 fastq 数据转换为 FRG 格式
一般情况下是利用错误率较低的 Illumina 数据来对 Pacbio 数据进行校正。 因此需要将 Illumina 的 Fastq 文件转换成 FRG 格式。 FRG 文件是文本文件，其中包含有 fastq 文件的路径和一些测序信息，比如插入片段长度、测序方向和碱基质量格式等。
fastqToCA 的常用示例和常用参数：

$ fastqToCA -insertsize 500 50 -libraryname pe500 -mates reads.1.fastq,reads.2.fastq > pe500.frg

-insertsize
    设置插入片段长度，为2个值，前者为插入片段长度，后者为标准差。
-libraryname
    文库的ID
-reads
    单端测序的 fastq 文件路径。
-mates
    双端测序的 fastq 文件路径。若双端数据交叉在一个文件中，则该参数后只有一个路径；若双端数据在2个fastq文件中，则2个路径用逗号分割。
-technology  default:illumina
    测序技术。可选的值有：none, sanger, 454, illumina(读长 < 160bp), illumina-long(任意读长), moleculo, pacbio-ccs, pacbio-corrected, pacbio-raw。
-type  default:sanger
    碱基质量格式。可选的值有：sanger, solexa, illumina。
-innie
    测序方向为 5'-3' <-> 3'-5'，适合小片段测序文库。
-outtie
    测序方向为 3'-5' <-> 5'-3'，适合大片段文库。

2. 使用 sffToCA 将 SFF 数据为 FRG 格式
Roche/454 数据格式为 SFF。
sffToCA 的常用示例和常用参数：

$ sffToCA  -libraryname LIB -output NAME IN.SFF ...

-insertsize
    设置插入片段长度，为2个值，前者为插入片段长度，后者为标准差。
-libraryname
    文库的ID。
-output
    输出文件的前缀。

3.2 使用 pacBioToCA/PBcR 进行 Pacbio 数据的修正

pacBioToCA 能对 Pacbio 数据进行校正，也能进行调用 runCA 进行基因组组装。
pacBioToCA 的输入文件为：

1. pacbio.spec   该文件包含一些参数信息，以 参数=值 表示。必须提供该文件，该文件内容为空，则全部使用默认参数。
2. pacbio.fastq  Pacbio测序的fastq文件，对该文件中的数据进行校正。
3. seq.frg        Illumina或454数据由 fatsqToCA 转换得到的 FRG 格式文件。如果不提供该文件，则进行自我校正。

pacBioToCA 的常用例子(示例数据)和参数：

$ pacBioToCA -libraryname out_prefix -s pacbio.spec -fastq pacbio.fastq seq1.frg seq2.frg

必须参数：
-libraryname STRING
    设定输出文件的前缀。
-s FILE
    从指定的参数文件获取额外的参数。若该文件内容为空，则使用默认的参数。
-fastq FILE
    输入需要校正的 Pacbio 数据。

可选参数：
-length INT    default:500
    pacBioToCA 对长度大于此值的测序 reads 进行修正。
-partitions INT    default:200
    将用于校正的数据分成指定份数，进行计算。内存不够，则需要增大该值。若内存足够，减小该值，可以减少计算时间。若该值太小(< 总线程数)，软件会中途运行runCorrection.sh报错并停止运行。
-threads INT    default：max
    使用的线程数。默认为所有可用的线程数。
-shortReads
    如果用于辅助校准的序列长度都 <= 100bp，则使用该参数。
-genomeSize INT    default:0
    设置基因组大小。设置改值后，在其它默认条件下，如果校正后的数据量 >25x，则会进行基因组组装。不设置该值，可能等同于设置该值为 5000000 （从软件输出信息中查知）。
-maxCoverage INT    default:40
    默认下，最多校正 40x 的数据量。
-maxGap INT
    默认下，如果有一段 pacibo-read 区域没有 short-read 覆盖，则会打断 pacbio-read。此值设定所能允许的没有 short-read 覆盖的最大长度。当然，若此段区域没有其它 pacbio-reads 重叠，则也会进行打断。

此外，在 pacbio.spec 文件中常设置的参数为：

若不进行基因，则需要设置：
assemble=0

软件默认设置最大内存使用量设置为系统最大可使用的内存量，若设定最大使用内存为 64G，则设置：
ovlMemory=64

使用 CA 8 版本对大基因组（>100Mbp）进行校正，设置：
maxGap=1500
ovlHashBlockLength=1000000000
ovlRefBlockLength=1000000000
使用Illumina数据进行校正：
blasr=-noRefineAlign -advanceHalf -noSplitSubreads -minMatch 10 -minPctIdentity 70 -bestn 24 -nCandidates 24
自我校正：
blasr=-minReadLength 200 -maxScore -1000 -maxLCPLength 16 -bestn 24 -nCandidates 24
ovlThreads=16
下面一个参数设置并行数。其值推荐设置为 总线程数 / ovlThreads值 。设置为2，表示每个 blasr 程序运行消耗 16 线程，并行运行 2 个 blasr 程序。
ovlConcurrency=2

若进行基因组组装。默认参数适用于单倍体数据或近亲杂交数据。对于大基因组或双倍体基因组，则需要设置：
asmUtgErrorRate=0.10
asmCnsErrorRate=0.10
asmCgwErrorRate=0.10
asmOBT=1
asmObtErrorRate=0.08
asmObtErrorLimit=4.5
utgGraphErrorRate=0.05
utgMergeErrorRate=0.05
ovlHashBits=24
ovlHashLoad=0.80