使用wtdbg2利用三代数据进行基因组de novo组装

wtdbg2软件下载与安装

$ wget https://github.com/ruanjue/wtdbg2/archive/v2.2.tar.gz -O ~/software_packages/wtdbg-v2.2.tar.gz
$ tar zxf /share/disk02_8T/home/chenlianfu/software_packages/wtdbg-v2.2.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/wtdbg2-2.2/
$ make -j 4
$ echo 'PATH=$PATH:/opt/biosoft/wtdbg2-2.2' >> ~/.bashrc
$ source ~/.bashrc

wtdbg2软件使用

进行基因组装
$ wtdbg2 -t 16 -i reads.fa.gz -fo prefix -L 5000

得到一致性序列
$ wtpoa-cns -t 16 -i prefix.ctg.lay -fo prefix.ctg.lay.fa

利用三代reads的比对结果对基因组序列进行打磨修正
$ minimap2 -t 16 -x map-pb -a prefix.ctg.lay.fa reads.fa.gz | samtools view -Sb - >prefix.ctg.lay.map.bam
$ samtools sort prefix.ctg.lay.map.bam prefix.ctg.lay.map.srt
$ samtools view prefix.ctg.lay.map.srt.bam | ./wtpoa-cns -t 16 -d prefix.ctg.lay.fa -i - -fo prefix.ctg.lay.2nd.fa

Linux系统下使用MegaCli软件对磁盘阵列进行操作

若管理的服务器比较多,服务器的硬盘经常会坏掉。可以使用megacli命令来查看磁盘阵列中的磁盘是否损坏。

1. CentOS6系统下安装MegaCli软件

在一些rpm网站中搜索megacli。然后下载megacli的rpm包进行安装:

# wget http://rpmfind.net/linux/Mandriva/devel/cooker/x86_64/media/non-free/release/megacli-8.02.21-1-mdv2012.0.x86_64.rpm
# rpm -ivh megacli-8.02.21-1-mdv2012.0.x86_64.rpm

megacli的使用需要依赖libsysfs库文件
# yum install libsysfs
# ln -s /usr/lib64/libsysfs.so.2.0.1 /usr/lib64/libsysfs.so.2.0.2

2. 使用megacli命令查看磁盘

# megacli -PDList -aALL 
此命令可以查看所有物理磁盘的信息。
若Media Error Count和Other Error Count的值大于0,则表示磁盘可能坏掉了。

# megacli -LDInfo -LALL -aAll 
此命令查看所有逻辑磁盘的信息。

使用高版本SMRT Analysis软件对Pacbio reads进行基因组De novo组装

1. SMRT Analysis简介

SMRT Analysis软件是由Pacbio Science公司开发的软件,现在是SMRT Link软件的一部分,用于对Pacbio Long Reads数据进行分析。其2.3版本仅能用于对Pacbio RSII测序仪获得的h5格式测序数据进行,该版本持续了较长一段时间。后来Pacbio测序仪增添Sequel平台后,SMRT Analysis软件更新到3.0,增加支持对Pacbio Sequel测序仪获得的bam格式数据进行分析。现在(2018.09.02)最新版本的SMRT Analysis到了5.1版本。

2. 下载并安装SMRT Link软件

将SMRT Link安装到CentOS 6.9系统中。安装软件前需要做一些准备:

首先,确定主机名是localhost并指向127.0.0.1的IP地址。刚安装完毕的系统一般都满足此要求。
保证/etc/hosts文件中含有以下一行信息:
127.0.0.1    localhost

其次,修改/etc/security/limits.conf文件,尾部增添如下4行:
chenlianfu      soft    nproc   10240
chenlianfu      hard    nproc   102400
chenlianfu      soft    nofile  10240   
chenlianfu      hard    nofile  102400
注意chenlianfu是安装SMRT Analysis软件所使用的用户的用户名。
该步骤用于增加chenlianfu用户的权限,这是与SMRT软件较高的资源消耗和较多的并行化相关的。
此外,注意此步骤过后,从新登陆

然后,修改/etc/sysconfig/iptables文件,增加5432(PostGresql数据库需要)、9090和8243(软件网页端使用tomcat提供服务时需要9090和8423)端口。在该文件中正确的位置增加如下两行:
-A INPUT -m state --state NEW -m tcp -p tcp --dport 5432 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 9090 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 8243 -j ACCEPT
然后重启防火墙,让修改生效:
# /etc/init.d/iptables restart

最后,启动PostGresql数据库。
# /etc/init.d/postgresql initdb
# /etc/init.d/postgresql start
# chkconfig postgresql on

从Pacbio Science官网下载SMRT Link软件

将SMRT Link安装到 ~/biosoft 目录下
$ mkdir ~/biosoft
$ cd ~/biosoft
$ wget https://downloads.pacbcloud.com/public/software/installers/smrtlink_5.1.0.26412.zip
$ unzip smrtlink_5.1.0.26412.zip
    要求输入密码,可以从下载该软件的网页中找到。
$ ./smrtlink_5.1.0.26412.run
    根据提示一步一步进行配置,基本都是直接使用默认值,按Enter键即可。
    设置软件临时文件夹的时候,要选择一个较大分区所对应的目录。我将其设置为/home/chenlianfu/biosoft/tmp/smrtlink
$ /home/chenlianfu/biosoft/smrtlink/admin/bin/services-start
    启动SMRT Link软件。

3. 使用HGAP4进行基因组De novo组装

1. SMRT Link软件必须要在Chrome浏览器中打开。在Chrome浏览器中输入安装有SMRT Link软件机器的IP地址,并接:9090,例如:192.168.30.1:9090,就可以连接到软件的用户登陆界面,输入用户名admin,密码admin则会进入软件。

2. 点击Data Management,再点击VIEW OR IMPORT SEQUENCE DATA,点击IMPORT并从下拉菜单中点击Sequel Sequence Data,从弹出的浏览设置中选择路径/home/chenlianfu/biosoft/smrtlink/install/smrtlink-release_5.1.0.26412/bundles/smrtinub/install/smrtinub-release_5.1.0.25847/private/pacbio/canneddata/lambdaTINY/m150404_101626_42267_c100807920800000001823174110291514_s1_p0.subreadset.xml,点击IMPORT,运行一段时间后,导入数据成功。

3. 返回到软件主界面,点击SMRT Analysis,点击CREATE NEW ANALYSIS,在Analysis Application下拉菜单中选择Assembly(HGAP 4),Analysis Name栏随便填写字符串Lambda_HGAP4,在Genome Length栏将基因组大小修改为58000,在Data Sets中勾选lambda/007_tiny,最后点击START,开始程序运行。

4. 值得注意的是公司给的测序数据文件夹中常常不包含sts.xml、adapters.fasta和scraps bam等文件,而在总的xml文件中却包含了这些信息,需要将其对应的行(第7行到第15行)删除掉后再运行HGAP4,否则程序运行会失败。

5. 程序运行结果在/home/chenlianfu/biosoft/smrtlink/userdata/jobs_root/目录下的一个数字编号(按运行顺序从000001开始编号)的文件夹中。该文件夹中的的主要结果文件(也可以在网页中下载):
./tasks/pbcoretools.tasks.contigset2fasta-0/file.fasta 最终的基因组组装结果。
./tasks/pbcoretools.tasks.gather_fasta-1/file.fasta 按某个长度

使用TeamViewer突破内网封锁远程连接服务器

1. TeamViewer作用

经常由于学校对端口的封锁,在校外难以使用OpenSSH连接学校内网的服务器。这个时候,可以在学校内网的服务器上开启图形化界面运行TeamViewer软件,同时在外网的笔记本电脑上同时运行TeamViewer软件。在两个软件上输入相同的用户名和密码,则可以通过该软件连接学校内网的服务器桌面,从而对内网服务器进行操作。

2. 内网服务器上安装TeamViewer软件

我的服务器系统是CentOS 6.9 64位系统,在TeamViewer网站上下载Linux 64位版本的TeamViewer软件。目前最新是V13版本的TeamViewer,貌似支持CentOS7。对于CentOS 6则下载V12版本的TeamViewer

# wget https://dl.tvcdn.de/download/version_12x/teamviewer_12.0.93330.i686.rpm
# yum install teamviewer_13.2.13582.x86_64.rpm

然后再内网服务器进入图形化桌面,在桌面上的各种菜单中寻找TeamViewer图标启动软件,输入用户名和密码登陆软件软件或记录其ID和密码。

3. 外网笔记本电脑上TeamViewer软件安装和使用

在TeamViewer官网下载windows系统的软件安装并运行。登陆用户名和密码进入软件控制远程服务器,或输入ID和密码来控制远程服务器。

使用ossutil命令下载诺禾致源释放的数据

1. 诺禾致源释放数据的内容

在诺禾致源释放数据的邮件中,可以找到如下示例内容:

AccessKeyId: LTAI2SeW6yzC54B2
AccessKeySecret: XTpXkjeqATP29sft8YE30ueN6UafCX 
预设OSS路径: oss://novo-data-nj/customer-Cbd6p59N/ 
区域: 华东1(杭州)

2. 在CentOS系统上安装ossutil软件

$ wget  wget http://docs-aliyun.cn-hangzhou.oss.aliyun-inc.com/assets/attach/50452/cn_zh/1524643963683/ossutil64 -O ~/bin/ossutil
$ chmod 755 ~/bin/ossutil

3. 使用ossutil下载数据

先配置用户名和密码。可以在邮件中找到ID和Key。
$ ossutil config -e oss://novo-data-nj/customer-Cbd6p59N/ -i LTAI2SeW6yzC54B2 -k XTpXkjeqATP29sft8YE30ueN6UafCX
再下载数据到当前目录。可以在邮件中找到数据文件路径。
$ ossutil cp -r -f --jobs 3 --parallel 2 oss://novo-data-nj/customer-Cbd6p59N/ ./
    -r 表示复制目录;-f表示强制覆盖已有的文件结果;--jobs 3 表示从3个位置下载同一个文件;--parallel 2 表示同时下载2个文件。

启用windows10的linux子系统

1. 首先,在windows10中启用开发者模式。按windows+x键,打开系统管理菜单;点击应用和功能,再点击主页;点击更新和安全,再点击开发者选项,勾选开发人员模式。在联网情况下,系统会安装高级开发功能。

2. 然后,启用windows10的linux子系统。按windows+x键,打开系统管理菜单,点击页面右边的程序和功能;在弹出的页面中点击启用和关闭windows功能,在继续弹出的页面中勾选适用于Linux的Windows子系统,点击确定,并重启windows10系统。

3. 最后,按windows键,点击并打开应用商店;在搜索栏搜索linux,然后选择一个喜欢的linux发行版本,比如debian或ubuntu等,进行下载和安装。

4. 我选择了debian进行下载和安装,最后,提示输入用户名和密码。该用户自动是超级管理员用户,可以使用sudo su -命令切换到root用户。

5. 使用windows系统下的linux子系统,相比于使用linux虚拟机更加方便,更少消耗内存。

6. window10下的debian子系统是最小安装的系统,很多软件无法使用。需要先对apt源进行修改后再安装软件。对/etc/apt/sources.list文件内容进行修改,将内容修改成如下内容(163源):

deb http://mirrors.163.com/debian/ jessie main non-free contrib
deb http://mirrors.163.com/debian/ jessie-updates main non-free contrib
deb http://mirrors.163.com/debian/ jessie-backports main non-free contrib
deb-src http://mirrors.163.com/debian/ jessie main non-free contrib
deb-src http://mirrors.163.com/debian/ jessie-updates main non-free contrib
deb-src http://mirrors.163.com/debian/ jessie-backports main non-free contrib
deb http://mirrors.163.com/debian-security/ jessie/updates main non-free contrib
deb-src http://mirrors.163.com/debian-security/ jessie/updates main non-free contrib

再安装软件:

# apt-get update
# apt-get install openssh-client
# apt-get install openssh-server

# apt-get remove vim-common
# apt-get install vim

使用github进行项目开发与程序维护

1. github简介

我对GitHub不甚了解。现在很多生物信息学软件的Paper中都提供了软件的Github网址。

个人理解,GitHub能很方便用于开源软件的存储和代码维护。软件开发人员在Github网站注册一个账户后,相当于在Github服务器上有了账户和密码;开发某个项目的软件后,可以将软件源代码通过该账号和密码上传到Github服务器上;每次软件维护和代码上传,都会详细记录更新的内容;其它人可以将代码下载,并进行修改,形成软件的分支。

2. 在Github软件上注册用户

Github注册页面注册一个账户,Username: train-chenlianfu;Email address:910432211@qq.com;Password: train123456;再点击Create an account,完成注册。

3. 安装git软件

在CentOS系统上,虽然自带git软件,但是推荐使用高版本git,否则可能版本过低,导致不能兼容https而不能使用。

下载并安装最新版本git软件:

$ wget https://github.com/git/git/archive/v2.17.0.tar.gz -O git-v2.17.0.tar.gz
$ tar zxf git-v2.17.0.tar.gz 
$ cd git-2.17.0/
$ make configure
$ ./configure --prefix=/opt/sysoft/git-2.17.0
$ make -j 24
$ make install
$ echo 'PATH=/opt/sysoft/git-2.17.0/bin/:$PATH' >> ~/.bashrc
$ source ~/.bashrc
$ cd .. && rm -rf git-2.17.0

4. 在GitHub创建一个项目

新注册账号会弹出需要验证邮箱地址的信息,在邮箱收件箱中点击链接,进入GitHub主页;

点击GitHub主页右上角的加号“+”,弹出一个菜单,点击New repository;

在Repository name栏输入一个软件名称,例如:blast_tools;

再勾选Initialize this repository with a README,表示给本项目直接生成一个文件README.md;

在Add a license下拉菜单中选择GNU General Public License v3.0,表示本项目是开源项目,软件中自动增加一个文件LICENSE;

点击Create repository,创建项目。

5. 使用git管理代码的上传和更新

从GitHub上下载blast_tools软件:

$ git clone https://github.com/train-chenlianfu/blast_tools.git
$ cd blast_tools/
$ ls
LICENSE  README.md

在当前工作目录下生成文件夹blast_tools,该文件夹中包含两个文件:LICENSE和README.md。
其中README.md中仅包含一行内容:“# blast_tools”;该文件内容能直接显示在GitHub网站的本项目主页上;#号表示该行使用加粗加黑加大字体显示。

其实,在当前文件夹下,还存在一个隐藏文件夹.git,该文件夹

在blast_tools/文件夹中写软件代码,例如编写了perl程序blast.pl,并修改了README.md文件内容。

点击GitHub主页最右上角的用户图标,在下拉菜单中选择Settings,进入用户设置;在左边一栏中点击SSH and GPG keys;再点击New SSH key;在Title一栏中随意输入字符train_ssh;再在Key一栏中输入OpenSSH的公钥(点击进入本博客OpenSSH公钥生成方法的介绍文章);最后点击Add SSH key;添加公钥后,则可以将本地的软件代码通过git命令同步到GitHub网址上。

将本地修改或增加的文件上传到GitHub上:

$ git config user.email "910432211@qq.com"
$ git config user.name "train-chenlianfu"
  以上两行命令设置该项目对应的GitHub项目的用户名和邮箱地址,用于修改.git/config配置文件内容。

$ git add *
  对当前工作目录中的所有文件提出更改。当然git add命令后可以接指定的多个文件或文件夹。

$ git rm -r files
  删除文件或文件夹。若直接使用系统的rm命令,则仅删除本地的。使用git rm命令则可以删除githhub上的文件。

$ git commit -m "first submit"
  提交本次改动,本次改动的所有文件,在GitHub中全部被批注成"first submmit"。
  程序会自动检测有有多少个文件有改动,插入多少字符,删除多少字符等。这些改动信息会存于本地.git文件夹中。

$ git push origin master
  将改动信息和文件推送到origin(GitHub项目)的master分支(一个项目可能有多个分支,此处将本地改动推送到主分支)上。可能会弹出界面需要输入用户名和密码。
  origin是一个远程项目,在.git/config配置文件中有定义origin指向了GitHub的网址https://github.com/train-chenlianfu/blast_tools.git。 

先写到这儿,后续分支什么的还没怎么理解。以后需要的时候再了解。

将WordPress博客迁移到全新Aliyun服务器上

1. 配置系统和基础软件安装

配置Aliyun服务器时,选择64位的CentOS 6操作系统。再使用yum安装系统软件:

# yum install httpd mysql mysql-devel mysql-server php php-mysql iptables*

生成防火墙配置文件/etc/sysconfig/iptables,内容如下:

# Firewall configuration written by system-config-firewall
# Manual customization of this file is not recommended.
*filter
:INPUT ACCEPT [0:0]
:FORWARD ACCEPT [0:0]
:OUTPUT ACCEPT [0:0]
-A INPUT -m state --state ESTABLISHED,RELATED -j ACCEPT
-A INPUT -p icmp -j ACCEPT
-A INPUT -i lo -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 22 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 80 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 8080 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 3306 -j ACCEPT
-A INPUT -j REJECT --reject-with icmp-host-prohibited
-A FORWARD -j REJECT --reject-with icmp-host-prohibited
COMMIT

重启防火墙

# /etc/init.d/iptables restart

初始化mysql数据库并启动mysql数据库

# /bin/cp /usr/share/mysql/my-large.cnf /etc/my.cnf
# /etc/init.d/mysqld start
# /usr/bin/mysql_secure_installation

修改httpd配置文件/etc/httpd/conf/httpd.conf,修改几处内容:

AddType text/html .shtml
AddType text/plain .sh
AddType text/plain .pl
AddType text/plain .py
AddHandler cgi-script .cgi
NameVirtualHost *:80

尾部添加:

<VirtualHost *:80>
    DocumentRoot /home/chenlianfu/wordpress
    ServerName www.chenlianfu.com
</VirtualHost>
<Directory "/home/chenlianfu/wordpress">
    AllowOverride None
    Options MultiViews FollowSymLinks ExecCGI
    Order allow,deny
    Allow from all
</Directory>

修改apache对用户chenlianfu家目录的访问权限

# chmod 750 /home/chenlianfu/
# usermod -aG chenlianfu apache

启动httpd服务

/etc/init.d/httpd start

2. 下载并安装WordPress

$ cd /home/chenlianfu/
$ wget https://cn.wordpress.org/wordpress-4.8.1-zh_CN.tar.gz
$ tar zxf wordpress-4.8.1-zh_CN.tar.gz

3. 迁移数据库

将原来wordpress的mysql数据库复制到/var/lib/mysql/目录下,即完成了数据的迁移。要注意修改文件权限:

# chown -R mysql:mysql /var/lib/mysql/wordpress_chenlianfu/
# chown -R apache:chenlianfu /home/chenlianfu/wordpress
# chmod 775 /home/chenlianfu/wordpress
# cd /home/chenlianfu/wordpress
# find ./ -perm 755 -exec chmod 775 {} \;
# find ./ -perm 644 -exec chmod 664 {} \;

再将域名解析指向Aliyun新服务器的IP后,输入www.chenlianfu.com网址,简单进行WordPress配置后则结束了。

4. 每天自动备份wordpress数据到邮箱

使用SAMR对蛋白组数据表达量进行差异分析

1. SAMR简介

SAM(Significance Analysis of Microarrays)在基因芯片数据时代中被开发出来进行基因表达量差异分析。该算法也能用于进行RNA-Seq数据的基因表达量差异分析,但貌似较少人会用它进行RNA-Seq数据分析。

最近在一篇对蛋白组数据差异分析软件进行比较的文章中,SAM结果表现最优。本文对基于R软件的SAM算法软件SAMR的使用进行简单讲述。

2. SAMR软件的安装和启动

打开R软件,输入如下命令安装SAMR相关的包:

$ R
> install.packages(c("samr", "matrixStats", "GSA", "shiny", "openxlsx"))
> source("http://bioconductor.org/biocLite.R")
> biocLite("impute")

要注意的是,对openxlsx包的安装可能会失败(我使用的是R-3.2.0版本),则在上述命令中去除对openxlsx的安装,选择手动下载并安装老版本的openxlsx包:

$ wget https://cran.r-project.org/src/contrib/Archive/openxlsx/openxlsx_2.4.0.tar.gz
$ R CMD INSTALL openxlsx_2.4.0.tar.gz

启动SAMR软件:

$ R
> library(shiny)
> runGitHub("SAM", "MikeJSeo")

启动SAMR软件,则会自动打开CentOS 6.8系统自带的火狐浏览器,进入软件的网页界面。

3. SARM软件使用

3.1 输入文件准备

软件的输入文件必须是xlsx格式的EXCEL文件。进行蛋白组表达量数据进行分析,其文件内容要求如下:

1. 第一列是基因Name,第二列是基因ID,每个基因ID应该独一无二。
2. 第一行表示样品名,第一行第一列和第一行第二列是空着的,从第三列开始表示样品名,且样品名仅能是数字1或2,代表两个不同的样品。都使用数字1表示样品1的多个重复,都使用数字2表示样品2的多个重复。虽然SARM能分析多个样品的数据,但是其结果格式不好,一般是进行两两比较。
3. 若有数据缺失,则某行某列不填入数据,或填入非数字(推荐 NA)代替。
4. 若有多组样品进行比较,则需要准备多个excel文件,而不是一个excel文件,多个sheet。

3.2 在网页中进行SAMR操作

按如下步骤进行操作:

1. 点击 Browse,选中目标excel文件。
2. 在 Minimum fold change 一栏中填写一个最小的差异倍数值,比如 1.5 或 2 等。
3. 在 Data type 一栏选择默认的 Array 。
4. 在 Response Type 一栏选择 Two class unpaired。这个选择项有很多,与输入文件的数据格式相关,我们是对两个样品进行比较,故选择此选项。输入文件第一列不是1和2,是不同的数值(比如血压),此处则选择Quantitative。
5. 在 Analysis Type 一栏选择默认的 Standard 即可。
6. 在 Test Statistic 一栏选择默认的 T-statistic 即可。
7. 在 Median Center the arrarys 一栏选择 Yes。选择Yes表示软件会根据中位数来对数据进行标准化,即让每一列的中位数都变成零。这样每个样品的数据都具有相同的中位数,相当于进行了样品间的标准化。软件推荐在输入数据之前进行标准化(不如TMM算法进行标准化),若输入的数据是标准化后的数据,则此栏选择默认的 No。
8. 其它参数选择一般情况下选择默认的即可。
9. 点击左上角的 Run 按钮进行分析,分开得到结果,在右边各个标签栏进行查看结果。
10. 在 Delta Table 页面下查看 Median FDR 值,从上往下找到该值变为0.05时的delta值,再在左侧 Delta 一栏选择 manually Enter Delta, 然后再在 Delta value 一栏填入该值,再进行计算。然后点击Current Settings标签页,看其False Discovery Rate的值,手动调整 Delta 值,直到 FDR 值最接近并低于 0.05 为止。Delta值越大,FDR值越小。
11. 在 Paste the filepath to save the output 一栏填入需要输出的文件夹路径,该路径一定要存在,在其下一栏填入输出文件的前缀,比如 out 。 在点击 Save,然后会输出结果文件 out.xlsx,该excel文件有多个sheet分别多个标签页中的结果。

4. SARM 软件算法和原理

这个比较复杂,我也不怎么搞懂,包括结果中一些Delta, Score(d), q-value, FDR, localfdr等,有点糊涂。等以后有时间搞更明白了,再添加解释。

Delta value:
Score(d):T检验的的d值,由 Numerator(r) / Denominator(s+s0) 得到。

FALCON

1. 安装FALCON

$ cd /opt/biosoft
$ export GIT_SYM_CACHE_DIR=~/.git-sym-cache # to speed things up
$ git clone git://github.com/PacificBiosciences/FALCON-integrate.git
$ cd FALCON-integrate
$ git checkout master  # or whatever version you want
$ git submodule update --init # Note: You must do this yourself! No longer via `make init`.
$ make init
$ source env.sh
$ make config-edit-user
$ make -j all
$ make test  # to run a simple one

2. FALCON运行示例

先设置好环境变量
$ cd /opt/biosoft/FALCON-integrate/
$ source env.sh

下载基因组大小只有200kb的greg200k数据,以及基因组大小有4.6Mb的E. coli数据
$ cd FALCON-examples/.git-sym/
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/greg200k-sv2.tar | tar xvf -
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar | tar xvf -

对200kb基因组进行组装
$ cd /opt/biosoft/FALCON-integrate/FALCON-examples/run/greg200k-sv2
$ fc_run.py fc_run.cfg

对4.6Mb基因组进行组装
$ cd /opt/biosoft/FALCON-integrate/FALCON-examples/run/ecoli
$ perl -p -i -e 's/^#job_type/job_type/' fc_run.cfg
$ fc_run.py fc_run.cfg

生成的最终结果文件为 2-asm-falcon/p_ctg.fa 。

3. FALCON工作原理

1. 对原始subreads进行Overlapping分析
2. 进行预组装,对Subreads进行校正
3. 对校正后的Subreads进行Overlapping分析
4. 过滤不可靠的Overlaps结果
5. 根据可靠的Overlap构建Grahp
6. 由Grahp得到contigs序列

4. FALCON的配置文件

FALCON软件的输入文件是Pacbio测序得到的Fasta文件。将这些Fasta文件的路径写入到文件 input.fofn 中用于软件输入。

FALCON软件运行的难点在于其配置文件fc_run.cfg的编写。一个fc_run.cfg示例如下:

[General]
## 设置任务提交方式
# jobtype的默认值应该是SGE,表示使用SGE集群进行计算;
# local表示使用本地主机运行FALCON,兼容性最好,毕竟集群的部署非常麻烦;
job_type = local

## 数据输入
# 设置输入文件为input.fofn,该文件中包含有PacBio数据的fasta文件。此外也可以输入dexta格式的文件。dexta格式是Pacbio测序结果h5的一种压缩结果,能极大压缩测序数据文件的大小。可以使用FALCON软件附带的命令dexta将fasta文件转换为dexta格式。
input_fofn = input.fofn
# 设置输入数据为原始测序数据,或者为修正后的数据。
input_type = raw
#input_type = preads

## 对PacBio数据进行校正:对Pacbio raw subreads进行overlapping、consensus和pre-assembly分析,对reads进行校正。
# 选择长度高于指定阈值的reads进行分析
length_cutoff = 12000
# 选择长度高于指定阈值的reads进行预组装
length_cutoff_pr = 12000
# 调用fasta2DB命令将reads数据分割成多份Blocks。
# -s50 参数表示每份数据含有50Mbp的数据量,该参数默认值为200,默认参数情况下,使用daligner进行比对需要消耗约16G内存。
# -x500 参数表示read长度低于500bp的会被忽略掉。
# -a 参数表示来同一个自零级波导孔的第二个read不会被忽略掉。
pa_DBsplit_option = -x500 -s50
# 调用daligner对所有Blocks进行Overlapping分析。
# -v 参数表示输出daligner程序运行信息
# -B4 参数表示每个daligner命令对4个Blocks进行Overlapping分析。该参数的值越大,则每个daligner命令的计算量越大,但是总的任务数越少。该参数等同于以前的-dal
# -k -w -h 参数设置匹配的kmer相关参数,其默认值分别为 14,6,35 。
# -T4 表示每个daligner使用4个线程进行计算,该参数默认值是4,该参数可以设置成2,4,8,16,32...。
# -M32 表示每个daligner命令使用32G内存,加入该参数起到限制内存使用的作用,对大基因组比较有用。
# -t16 参数表示过滤掉覆盖度高于16的kmer,这些kmer会导致内存使用过多。默认设置下,daligner可以根据-M参数的值自动计算本参数的值。
# -l1000 参数表示忽略长度低于1000bp的reads。
# -s1000 参数表示记录比对结果时以每1000bp为一个记录点,相比于默认值100,能少很多记录点。
# 使用daligner的默认参数能很好地处理raw pacbio数据。
# 而对corrected pacbio数据,推荐使用-k20 -h60 -e.85参数。
pa_HPCdaligner_option = -v -B4 -k14 -T8 -t16 -e.70 -l1000 -s1000
# FALCON使用fc_consensus.py调用C语言写的程序来根据daligner进行Overlapping分析的结果进行consensus分析,从而对subreads进行校正。
# --min_cov 参数设置当read序列上某位点覆盖度低于指定阈值时,对read进行打断或截短,默认值是6。
# --min_cov_aln 参数设置当read序列的平均覆盖度低于指定阈值时,直接略过该read,默认值是10。
# --min_n_read 和 --max_n_read 参数设定比对结果中包含的reads数在此范围内才能得到consensus结果,其默认值分别是10和500。对于基因组重复程度较高的情况,要设置较低的--max_n_read值来减少对重复区域进行consensus分析的计算消耗。
# --min_idt 参数设置最小identity的比对结果能用于reads校正。
# --n_core 参数设置允许的线程数,默认值是24。
falcon_sense_option = --output_multi --min_idt 0.70 --min_cov 4 --max_n_read 200 --n_core 8
# 设置运行daligner任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程,因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源,容易导致宕机。
da_concurrent_jobs = 10
la_concurrent_jobs = 2
# 设置运行fc_consensus.py任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程,因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源,容易导致宕机。
cns_concurrent_jobs = 20
# 设置在SGE集群运行的并发数
sge_option_da = -pe smp 8 -q jobqueue
sge_option_la = -pe smp 2 -q jobqueue
sge_option_fc = -pe smp 80 -q jobqueue
sge_option_cns = -pe smp 16 -q jobqueue

## 对校正后的reads进行overlapping分析,其参数和上一个步骤的参数一致。
ovlp_DBsplit_option = -x500 -s50
ovlp_HPCdaligner_option = -v -B4 -k20 -h60 -T8 -t32 -e.96 -l500 -s1000
# 设置对校正后reads运行daligner任务的并发数。
pda_concurrent_jobs = 10
pla_concurrent_jobs = 2
sge_option_pda = -pe smp 8 -q jobqueue
sge_option_pla = -pe smp 2 -q jobqueue

## 过滤overlaps
# 若reads首尾两端的覆盖度比平均覆盖度大很多,则表明reads首尾是重复序列;若reads首尾两端的覆盖度比平均覆盖度相差较小很多,则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。该步骤的参数设置不对,容易导致overlaps全部被过滤掉,得不到基因组组装的结果。
# --bestn设置报告reads此数目的最优overlaps。
# --min_cov和--max_cov表示所允许的reads首尾的覆盖度范围。对于通过length_cutoff得到 >20x 校正后的数据进行的基因组组装,可以设置--min_cov值为5,设置--max_cov为平均覆盖度的3倍。若用于基因组组装的数据量较少,可以设置该值为1或2。
# --max_diff设置所允许的首尾覆盖度值的最大差异。设置该参数的值为平均覆盖度的2倍。
# 对于特殊情况,可以设置更高的--max_cov和--max_diff值。
# 可以使用在1-preads_ovl目录下运行 fc_ovlp_stats.py --fofn merge-gather/las.fofn 导出overlap结果首尾的覆盖度结果,从而帮助设置以上参数
overlap_filtering_setting = --max_diff 50 --max_cov 75 --min_cov 5 --bestn 10

5. FALCON结果文件

0-rawreads/
    该目录存放对raw subreads进行overlpping分析与校正的结果;
    以job_前缀的文件夹的表示daligner任务的结果;
    以m_前缀的文件夹表示整合任务的结果,对raw subreads进行overlapping分析的比对结果是las文件,位于这些文件夹中;
    merge-gather中存放着las.fofn文件,该文件中存放着las文件的路径信息;
    preads目录下保存着校正后的reads信息,位于cns_前缀目录下的fasta文件。

1-preads_ovl/
    该目录存放对校正后reads进行overlapping的结果;
    以job_前缀的文件夹的表示daligner任务的结果;
    以m_前缀的文件夹表示整合任务的结果,对校正后reads进行overlapping分析的比对结果是las文件,位于这些文件夹中;
    merge-gather中存放着las.fofn文件,该文件中存放着las文件的路径信息;

2-asm-falcon/
    该目录是最终结果目录,包含draft contigs结果文件p_ctg.fa。