FileZilla在CentOS系统上的安装

由于编译或使用新版本需要高版本的 GCC 和 wxWidgets,因此,不推荐使用新版本的 FileZilla。Fileilla 官网仅提供了最新版本的下载链接。可以到sourceforge上下载旧版本。

$ wget http://sourceforge.net/projects/filezilla/files/FileZilla_Client/3.5.3/FileZilla_3.5.3_x86_64-linux-gnu.tar.bz2
$ tar jxf FileZilla_3.5.3_x86_64-linux-gnu.tar.bz2 -C /opt/
$ echo 'PATH=$PATH:/opt/FileZilla3/bin/' >> ~/.bashrc
$ source ~/.bashrc 
$ filezilla

AGP格式简单说明

AGP文件为NCBI数据上传要求的标准格式,用来描述小片段序列(比如contig)如何构成大片段序列(比如scaffold和chromosome)。详细的说明文档请见:http://www.ncbi.nlm.nih.gov/projects/genome/assembly/agp/AGP_Specification.shtml
AGP文件有9列,分别是:

1. 大片段的序列名(object)
2. 大片段起始(object_begin)
3. 大片段结束(object_end)
4. 该段序列在大片段上的编号(part_number)
    一般一个大片段由多个小片段和gap组成。此处则为这些小片段和gap在大片段上的编号。
5. 该段序列的类型(component_type)
    常用的是W、N和U。W表示WGS contig;N表示指定大小的gap;U表示不明确长度的gap,一般用100bp长度。
6. 小片段的ID或gap长度(component_id or gap_length)
    如果第5列不为N或U,则此列为小片段的ID。
    如果第5列是N或U,则此列为gap的长度。如果第5列为U,则此列值必须为100。
7. 小片段起始或gap类型(component_begin or gap_type)
    如果第5列是N或U,则此列表示gap的类型。常用的值是scaffold,表示是scaffold内2个contigs之间的gap。其它值有:contig,2个contig序列之间的unspanned gap,这样的gap由于没有证据表明有gap,应该要打断大片段序列;centromere,表示中心粒的gap;short_arm,a gap inserted at the start of an acrocentric chromosome;heterochromatin,a gap inserted for an especially large region of heterochromatic sequence;telomere,a gap inserted for the telomere;repeat,an unresolvable repeat。
8. 小片段结束或gap是否被连接(component_end or linkage)
    如果第5列是N或U,则此列一般的值为yes,表示有证据表明临近的2个小片段是相连的。
9. 小片段方向或gap的连接方法(orientation or linkage_evidence)
    如果第5列不为N或U,则此列为小片段的方向。其常见的值为 +、-或?。
    如果第5列是N或U,则此列表明临近的2个小片段能连接的证据类型。其用的值是paired-ends,表明成对的reads将小片段连接起来。其它值有:na,第8列值为no的时候使用;align_genus,比对到同属的参考基因组而连接;align_xgenus,比对到其它属的参考基因组而连接;align_trnscpt,比对到同样物种的转录子序列上;within_clone,gap两边的序列来自与同一个clone,但是gap没有paired-ends跨越,因此这种连接两边小片段无法确定方向和顺序;clone_contig,linkage is provided by a clone contig in the tiling path (TPF);map,根据连锁图,光学图等方法确定的连接;strobe,根据PacBio序列得到的连接;unspecified。如果有多中证据,则可以写上多种证据,之间用分号分割。

例子:
Scaffold from component (WGS)
Chromosome from scaffold (WGS)

Perl多线程的简单例子

perl多线程使用的一个简单例子:同时对多个ip进行ping的命令,检测各个ip地址的响应速度。

#!/usr/bin/perl
use strict;
# 声明调用的模块
use Thread;

my @ips = ("221.141.1.222", "89.46.101.122", "199.200.120.37");

my @threads;
# 创建线程并push到数组中。第一个参数是一个子程序的名字,后面的参数是输入到该子程序的数组内容。
foreach (@ips) {
    push @threads, threads->create(\&ping,$_);
}

# 子程序join到主线程,取得返回值。
foreach (@threads) {
    $_->join();
}

# 子程序,对ip地址ping5次,得到其平均值。
sub ping {
    $_ = shift;
    my $ping_log = `ping -c 5 $_`;
    my @ping_log = split /\n/, $ping_log;
    my $time = 0;
    foreach (@ping_log) {
        $time += $1 if /time=(.*) ms/;
    }
    my $avg_time = $time / 5;
    $time = "na" if $time == 0;
    print "$_\t$avg_time\n";
}

Circos的安装和简单使用

1. Circos简介

Circos用于画基因组圈图。其帮助文档:http://www.circos.ca/documentation/tutorials/

2. Circos的安装

Circs是perl写出来的程序,其正常使用需要依赖于一些Perl模块,特别是GD。GD的安装如下:

先安装libgd:
# unzip libgd-gd-libgd-00cd9583242e.zip 
# cd libgd-gd-libgd-00cd9583242e/
# ./bootstrap.sh 
# ./configure && make -j 4 && make install 
# gdlib-config
# cd ..
# rm -rf libgd-gd-libgd-00cd9583242e

再安装GD
# ln -s /usr/lib64/libgd.so.2.0.0 /usr/lib64/libgd.so.3
# cpan -i GD

下载circos并安装:

$ wget http://circos.ca/distribution/circos-tools-0.18.tgz
$ wget http://circos.ca/distribution/circos-0.66.tgz
$ wget http://circos.ca/distribution/circos-tutorials-0.66.tgz
$ wget http://circos.ca/distribution/circos-course-0.67.tgz

$ tar zxf circos-0.66.tgz
$ cd circos-0.66/bin
$ ./list.modules        该命令列出circos所需求的perl模块
$ ./test.modules        该命令检测这些所需求的perl模块是否正确安装上
$ sudo cpan -i Config::General Font::TTF List::MoreUtils Math::Bezier Math::Round Math::VecStat Params::Validate Readonly Regexp::Common Set::IntSpan Text::Format Clone
$ ./gddiag              检测GD是否能用于画图,最终能生成diag.png文件
$ ./circos --help       能给出简单的命令帮助文件,否则安装不成功
$ cd ../example/
$ ../bin/circos -conf etc/circos.conf  一个使用的例子

3. Circos的配置文件准备

Circos的使用主要通过输入一个配置文件。该配置文件的内容格式主要以各种区块表示,大区块中可以包含小区块。区块中以“变量 = 值”的方式来进行参数的设定。例如:

<links>

<link>
 file      = data/set1.txt
 color     = black
 ...
</link>

<link>
 file      = data/set2.txt
 color     = red
 ...
</link>

</links>

此外,有些配置信息一般不需要改动,比如颜色,字体等。我们一般将这类信息保存到一个独立的配置文件中。只需要在主配置文件中声明包含这些独立的配置文件名,即表示使用其配置信息。例如,最常用的放置到主配置文件尾部的数行:

设置生成的图片参数
<image>
<<include etc/image.conf>>
</image>
设置颜色,字体,填充模式的配置信息:
<<include etc/colors_fonts_patterns.conf>>
系统与debug参数:
<<include etc/housekeeping.conf>>

4. Circos的使用参数

-version
    查询circos版本
-modules
    检测perl模块
-conf <string>
    输入主配置文件
-outputdir <string>
    设置输出文件的路径
-outputfile <string>
    设置输出文件名,该参数的值以.png为后缀
-svg
    生成svg结果文件
-nosvg
    不生成svg结果文件

5. Circos配置文件详解

Circos的命令使用简单,但是配置文件极其复杂,以下从各个track进行详解

5.1 ideogram block 显示染色体

将染色体在圈图上展示出来,代表每个染色体的图形,称为ideogram。将以下配置信息放入一个单独的配置文件中,给其命名 ideogram.conf 。

<ideogram>

## 设定 ideograms 之间的空隙
<spacing>
# 设置圈图中染色体之间的空隙大小,以下设置为每个空隙大小为周长的 0.5%
default = 0.005r

# 也可以设置指定两条染色体之间的空隙
#<pairwise hsY;hs1>
# 以下设定为两条染色体之间的空隙约为圆的 20 度角。
#spacing = 20r
#</pairwise>

</spacing>

## 设定 ideograms 
# 设定 ideograms 的位置,以下设定 ideograms 在图离圆心的 90% 处
radius           = 0.90r
# 设定 ideograms 的厚度,可以使用 r(比例关系) 或 p(像素)作为单位
thickness        = 20p
# 设定 ideograms 是否填充颜色。填充的颜色取决于 karyotype 指定的文件的最后一列。
fill             = yes
# 设定 ideograms 轮廓的颜色及其厚度。如果没有该参数或设定其厚度为0,则表示没有轮廓。
stroke_color     = dgrey
stroke_thickness = 2p

## 设定 label 的显示
# 设定是否显示 label 。 label 对应着 karyotype 文件的第 4 列。如果其值为 yes,则必须要有 label_radius 参数来设定 label 的位置,否则会报错并不能生成结果。
show_label       = yes
# 设定 label 的字体
label_font       = default
# 设定 label 的位置
label_radius     = 1r+90p
# 设定 label 的字体大小
label_size       = 40
# 设定 label 的字体方向,yes 是易于浏览的方向。
label_parallel   = yes

</ideogram>

5.2 ticks block 以刻度形式显示染色体大小

将染色体的大小以刻度的形式在圈图上展示出来。将以下配置信息放入一个单独的配置文件中,给其命名 ticks.conf 。

# 是否显示 ticks
show_ticks         = yes
# 是否显示 ticks 的 lables
show_tick_labels    = yes

## 设定 ticks
<ticks>
## ticks 的设置
# 设定 ticks 的位置
radius           = 1r
# 设定 ticks 的颜色
color            = black
# 设定 ticks 的厚度
thickness        = 2p
# 设定 ticks' label 的值的计算。将该刻度对应位置的值 * multiplier 得到能展示到圈图上的 label 值。
multiplier       = 1e-6
# label 值的格式化方法。%d 表示结果为整数;%f 结果为浮点数; %.1f 结果为小数点后保留1位; %.2f 结果为小数点后保留2位。
format           = %d

## 以下设置了 2 个 ticks,前者是小刻度,后者是大刻度。
<tick>
# 设置每个刻度代表的长度。若其单位为 u,则必须要设置 chromosomes_units 参数。比如设置 chromosomes_units = 1000000,则如下 5u 表示每个刻度代表 5M 长度的基因组序列。
spacing        = 5u
# 设置 tick 的长度
size           = 10p
</tick>

<tick>
spacing        = 25u
size           = 15p
# 由于设置的是大刻度,以下用于设置展示 ticks' label。
show_label     = yes
# 设置 ticks' label 的字体大小
label_size     = 20p
# 设置 ticks' label 离 ticks 的距离
label_offset   = 10p
format         = %d
</tick>

</ticks>

5.3 links block 以曲线连接显示基因组内部区域之间的联系

基因组内部不同的序列区域之间有联系,将之使用线条进行连接,从而展示到圈图上。常见的是重复序列之间的连接。将以下配置信息放入一个单独的配置文件中,给其命名 links.conf 。

<links>

<link>
# 指定 link 文件的路径,其文件格式为:
# chr1  start   end     chr2    start   end
# hs1   465     30596   hs2     114046768       114076456
# 表明这两个染色体区域有联系,例如这个区域的序列长度>1kb且序列相似性>=90%。
file          = data/5/segdup.txt
# 设置 link 曲线的半径
radius        = 0.8r
# 设置贝塞尔曲线半径,该值设大后曲线扁平,使图像不太好看。
bezier_radius = 0r
# 设置 link 曲线的颜色
color         = black_a4
# 设置 link 曲线的厚度
thickness     = 2

<rules>
# 以下可以设置多个 rules,用来对 link 文件的每一行进行过滤或展示进行设定。每个 rule 都有一个 condition 参数;如果该 condition 为真,除非 flow=continue ,则不

# 如果 link 文件中该行数据是染色体内部的 link,则不对其进行展示
<rule>
condition     = var(intrachr)
show          = no
</rule>

# 设置 link 曲线的颜色与 ideogram 的颜色一致,否则为统一的颜色。
<rule>
# condition 为真,则执行该 block 的内容
condition     = 1
# 设置 link 曲线的颜色为第 2 条染色体的颜色。对应这 link 文件中第 4 列数据对应的染色体的名称
color         = eval(var(chr2))
# 虽然 condition 为真,但依然检测下一个 rule
flow          = continue
</rule>

# 如果 link 起始于 hs1,则其 link 曲线半径为 0.99r
<rule>
condition     = from(hs1)
radius1       = 0.99r
</rule>

# 如果 link 结束于 hs1,则其 link 曲线半径为 0.99r
<rule>
condition     = to(hs1)
radius2       = 0.99r
</rule>

</rules>

</link>

</links>

5.4 plots block 以直方图形式展示数据

将基因组序列的GC含量,表达量等以直方图的形式在圈图中展示出来。将以下配置信息放入一个单独的配置文件中,给其命名 plots_histogram.conf 。以下作了两个直方图,并对分别添上背景或网格线。

<plot>
# 设定为直方图
type = histogram
# 数据文件路径,为 4 列:
# chromosome	start	end	data
# hs1	0	1999999	180.0000
file = data/5/segdup.hs1234.hist.txt
# 设置直方图的位置,r1 要比 r0 大。直方图的方向默认为向外。
r1   = 0.88r
r0   = 0.81r
# 直方图的填充颜色
fill_color = vdgrey
# 默认下直方图轮廓厚度为 1px,若不需要轮廓,则设置其厚度为0,或在 etc/tracks/histogram.conf 中修改。
thickness = 0p
# 直方图是由 bins (条行框)所构成的。若 bins 在坐标上不相连,最好设置不要将其bins连接到一起。例如:
# hs1 10 20 0.5
# hs1 30 40 0.25
# 上述数据设置值为 yes 和 no 时,图形是不一样的。
extend_bin = no

# 以下添加 rule ,不在 hs1 上添加直方图。
<rules>
<<include exclude.hs1.rule>>
</rules>

# 设定直方图的背景颜色
<backgrounds>
show  = data

<background>
color = vvlgrey
</background>
<background>
color = vlgrey
y0    = 0.2r
y1    = 0.5r
</background>
<background>
color = lgrey
y0    = 0.5r
y1    = 0.8r
</background>
<background>
color = grey
y0    = 0.8r
</background>

</backgrounds>

</plot>

<plot>
type = histogram
# 此处直方图的数据文件第 4 列是多个由逗号分割的数值,需要制作叠加的直方图。
file = data/5/segdup.hs1234.stacked.txt
r1   = 0.99r
r0   = 0.92r
# 给 4 个值按顺序填充不同的颜色
fill_color  = hs1,hs2,hs3,hs4
thickness = 0p
orientation = in
extend_bin  = no

<rules>
<<include exclude.hs1.rule>>
</rules>

# 在直方图中添加坐标网格线
<axes>
show = data
thickness = 1
color     = lgrey

<axis>
spacing   = 0.1r
</axis>
<axis>
spacing   = 0.2r
color     = grey
</axis>
<axis>
position  = 0.5r
color     = red
</axis>
<axis>
position  = 0.85r
color     = green
thickness = 2
</axis>

</axes>

</plot>

5.5 plots block 以热图形式显示数据

基因组一个区域内有多组数据时,适合以热图形式显示数据。比如基因表达量。将以下配置信息放入一个单独的配置文件中,给其命名 plots_heatmap.conf 。

<plot>
# 绘制 heat map
type  = heatmap
# 设定数据文件路径。文件有 5 列
# chrID start   end     data    class
# hs1 0 1999999 113.0000 id=hs1
# hs1 0 1999999 40.0000 id=hs4
# hs1 0 1999999 20.0000 id=hs2
# hs1 0 1999999 7.0000 id=hs3
file  = data/5/segdup.hs1234.heatmap.txt
# 设定图形所处位置
r1    = 0.89r
r0    = 0.88r
# 设定热图的颜色。颜色为 hs3 ,以及相应带不同透明程度的 5 种颜色。
color = hs1_a5,hs1_a4,hs1_a3,hs1_a2,hs1_a1,hs1
# 设定 scale_log_base 参数。计算颜色的方法如下:
# f = (value - min) / ( max - min )    热图中每个方块代表着一个值,并给予相应的颜色标示。一系列的值 [min,max] 对应一系列的颜色 c[n], i=0..N
# n = N * f ** (1/scale_log_base)
# 由上面两个公式计算出代表颜色的 n 值。
# 若 scale_log_base = 1,则数值与颜色的变化是线性的;
# 若 scale_log_base > 1,则颜色向小方向靠近;
# 若 scale_log_base < 1,则颜色向大方向靠近。
scale_log_base = 5

<rules>
<<include exclude.hs1.rule>>

# 仅显示 id = hs1 的数据
<rule>
condition = var(id) ne "hs1"
show      = no
</rule>

</rules>

</plot>

<plot>
type  = heatmap
file  = data/5/segdup.hs1234.heatmap.txt
r1    = 0.90r
r0    = 0.89r
color = hs2_a5,hs2_a4,hs2_a3,hs2_a2,hs2_a1,hs2
scale_log_base = 5

<rules>
<<include exclude.hs1.rule>>

<rule>
condition = var(id) ne "hs2"
show      = no
</rule>

</rules>

</plot>

<plot>
type  = heatmap
file  = data/5/segdup.hs1234.heatmap.txt
r1    = 0.91r
r0    = 0.90r
color = hs3_a5,hs3_a4,hs3_a3,hs3_a2,hs3_a1,hs3
scale_log_base = 5

<rules>
<<include exclude.hs1.rule>>

<rule>
condition = var(id) ne "hs3"
show      = no
</rule>

</rules>

</plot>

<plot>
type  = heatmap
file  = data/5/segdup.hs1234.heatmap.txt
r1    = 0.92r
r0    = 0.91r
color = hs4_a5,hs4_a4,hs4_a3,hs4_a2,hs4_a1,hs4
scale_log_base = 5

<rules>
<<include exclude.hs1.rule>>

<rule>
condition = var(id) ne "hs4"
show      = no
</rule>

</rules>

</plot>

5.6 plots block 以文本形式显示数据

若需要在圈图上显示一些基因的名称,此时需要以文本形式显示数据。将以下配置信息放入一个单独的配置文件中,给其命名 plots_text.conf 。

<plot>
# 表示出文字
type  = text
# 数据文件路径
file  = data/6/genes.labels.txt
# 显示在图形中的位置
r1    = 0.8r
r0    = 0.6r
# 标签的字体
label_font = light
# 标签大小
label_size = 12p
# 文字边缘的大小,设置较小则不同单词就可能会连接到一起了。
# padding  - text margin in angular direction
# rpadding - text margin in radial direction
rpadding   = 5p
# 设置是否需要在 label 前加一条线,用来指出 lable 的位置。
show_links     = no
link_dims      = 0p,2p,5p,2p,2p
link_thickness = 2p
link_color     = black

<rules>
<<include exclude.hs1.rule>>

# 设置 rule ,对 label 中含有字母 a 或 b 的特异性显示
<rule>
condition  = var(value) =~ /a/i
label_font = bold
flow       = continue
</rule>
<rule>
condition  = var(value) =~ /b/i
color      = blue
</rule>
</rules>

</plot>

5.7 rules block 放置常用的规则配置

本例子中,很多track没有在1号染色体上展示,需要设置如下规则信息,将之写入到文件 exclude.hs1.rule 中

<rule>
condition = on(hs1)
show      = no
</rule>

5.8 主配置文件

在主配置文件 circos.conf 中,包含以上所需要的配置文件信息,则可以画出所需要的track。此外,可以设置一些全局的设置。

# 指定染色体组型的文件,该文件分为
karyotype = data/karyotype/karyotype.human.txt

# 设置长度单位,以下设置表示 1M 长度的序列代表为 1u。
chromosomes_units = 1000000

# 默认设置下是将 karyotype 文件中所有的染色体都展示出来。当然,也可能根据需要仅展示指定的 chromosomes, 使用如下的参数进行设置。
chromosomes_display_default = no
# 以下参数设置指定的 chromosomes 用于展示到圈图中。// 中是一个正则表达式,匹配的 chromosomes 用于展示到圈图中。其匹配的对象是 karyotype 文件中的第 3 列。也可以直接列出需要展示的 chromosomes, 例如:hs1;hs2;hs3;hs4 。
chromosomes                 = /hs[1-4]$/
# chromosomes                 = hs1;hs2;hs3;hs4

# 以下设置各个 ideograms 的大小。其总长度为 1 ,hs1 的长度为 0.5, hs2,hs3 和 hs4 这 3 个 chromosomes 的总长度为 0.5,并且这 3 个 chromosomes 的长度是分布均匀的。注意前者的单位是 r, 后者使用了正则表达式对应多个 chromosomes, 其单位于是为 rn 。
chromosomes_scale   = hs1=0.5r,/hs[234]/=0.5rn

# 使 hs2, hs3 和 hs4 在圈图上的展示方向是反向的。
chromosomes_reverse = /hs[234]/

# 设置各个 ideograms 的颜色
chromosomes_color   = hs1=red,hs2=orange,hs3=green,hs4=blue

# 默认下在 ideogram block 中统一设置了 ideogram 的位置,可以使用此参数调整指定 ideogram 的位置。
chromosomes_radius  = hs4:0.9r
# chromosomes_radius  = hs2:0.9r;hs3:0.8r;hs4:0.7r

# karyotype 文件最后一列指定了各个 chromosomes 的颜色,而使用 chromosomes_color 参数也能修改颜色。当然,使用如下方式进行颜色的修改,则更加直观。以下方式是对颜色重新进行定义。chr1,chr2,chr3 和 chr4 对应着 karyotype 文件最后一列的值,代表着颜色的类型。此处使用 color block 来对其进行重新定义。注意重新定义的时候需要加符号 * 
<colors>
chr1* = red
chr2* = orange
chr3* = green
chr4* = blue
</colors>

### 绘制 plot 图
<plots>

<<include plots_histogram.conf>>
<<include plots_heatmap.conf>>
<<include plots_text.conf>>

</plots>

<<include ideogram.conf>>
<<include ticks.conf>>
<<include links.conf>>

################################################################
# 插入必须的并不常修改的标准参数
<image>
<<include etc/image.conf>>
</image>
<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

5.9 使用 circos 命令画图

对配置文件设置完毕后,使用命令进行画图

$ ./bin/circos -conf circos.conf 

结果如下:
circos.png

使用 abyss 进行 scaffolding

1. ABySS 进行 scaffolding 的目的与优点

目的: 对其它基因组 denovo 的 assembly 结果,使用 abyss 再进行一次 scaffolding。
优点: 可以使用 RNA-seq 的转录子数据进行基因组的辅助组装。

2. ABySS 进行 scaffolding 的命令行

输入文件: assembly.fasta, 2000.1.fastq, 2000.2.fastq, 5000.1.fastq, 5000.2.fastq。
输入文件是基因组的组装结果,和 3 对 mate-paired Illumina 数据。

2.1 对 assembly.fasta 进行序列改名

去除序列之间的换行
fasta_no_blank.pl assembly.fasta > 11; mv 11 assembly.fasta
给序列按顺序重命名
perl -e '$num = 0; while (<>) {if (/^>/) { s/>(.*)/>$num/; print; $num ++; } else { print } }' assembly.fasta > ledodes-6.fa

2.2 将 mate-paired 数据比对到基因组序列上

根据比对结果,得到 mate-paired library 的 insertSize 信息(以 .hist 为后缀的文件)和 序列之间的连接、距离与顺序信息 (以 .dist.dot 为后缀的 graph 文件)。

abyss-map -j24 -l87 3000.1.fastq 3000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp1-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp1-6.dist.dot mp1-6.hist
abyss-map -j24 -l87 8000.1.fastq 8000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp2-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp2-6.dist.dot mp2-6.hist

以上命令行中参数:

-j24
    使用 24 个线程运行
-l87
    使用的 kmer值 为 87
-s200
    sedd contigs的最小长度为 200bp
-n10
    所允许连接两条序列的最小的pairs的数目

2.3 进行 scaffolding

abyss-scaffold -k87 -s200 -n5 -g ledodes-6.path.dot ledodes-6.fa mp1-6.dist.dot mp2-6.dist.dot > ledodes-6.path
PathConsensus -k87 -p0.9 -s ledodes-7.fa -g ledodes-7.adj -o ledodes-7.path ledodes-6.fa ledodes-6.fa ledodes-6.path
cat ledodes-6.fa ledodes-7.fa \
  | MergeContigs -k87 -o ledodes-8.fa - ledodes-7.adj ledodes-7.path
ln -sf ledodes-8.fa ledodes-scaffolds.fa
PathOverlap --overlap --dot -k87 ledodes-7.adj ledodes-7.path > ledodes-8.dot

2.4 使用转录子序列进行 rescaffolding

bwa index ledodes-8.fa
bwa mem -a -t2 -S -P -k87 ledodes-8.fa transcripts.fasta \
  |gzip > long1-8.sam.gz
abyss-longseqdist -k87 long1-8.sam.gz \
  |grep -v "l=" >long1-8.dist.dot
abyss-scaffold -k87 -s200 -n1 -g ledodes-8.path.dot ledodes-8.dot long1-8.dist.dot > ledodes-8.path
PathConsensus -k87 -p0.9 -s ledodes-9.fa -g ledodes-9.adj -o ledodes-9.path ledodes-8.fa ledodes-8.dot ledodes-8.path
cat ledodes-8.fa ledodes-9.fa \
  | MergeContigs -k87 -o ledodes-10.fa - ledodes-9.adj ledodes-9.path
ln -sf ledodes-10.fa ledodes-long-scaffs.fa

磁盘IO检测工具 iostat 的简单使用说明

1. iostat 的常用例子和常用参数

iostat 的常用例子:

$ iostat -c -d 1 10
$ iostat -c -x -m 1 10
$ iostat -c -x /dev/sda /dev/sdb 1 10
上个命令表示对 /dev/sda 和 /dev/sdb 两个磁盘进行 I/O 统计,每秒统计一次,共统计10次。

iostat 的常用参数:

-c
    显示 CPU 的使用情况
-d
    显示设备的使用情况
-x
    显示扩展统计数据。若设置了 -x 参数,则 -d 参数失效
-k
    使用 Kb 作为单位
-m
    使用 Mb 作为单位

2. iostat 使用注意事项

1. 需要进行多次统计,多次统计的结果是实时的结果。

2. 使用 -c 参数能得到 CPU 的使用统计。 %iowait 表示系统有 I/O request 的 CPU 空闲时间。

3. 使用 -d 参数,其统计结果中 tps(the number of transfers per second) 表示设备每秒的传输次数。

4. 使用 -x 参数, %util 表示磁盘使用率,即处理 I/O request 的时间的比率。该值接近 100%,表明 I/O 成了瓶颈。若此时 CPU 的使用率比预期低,则表示磁盘性能不够,导致 CPU 不能充分利用。

使用 GCE 进行基因组大小评估

1. GCE 简介

GCE(Genome Characteristics Estimation) 是华大基因用于基因组评估的软件,其参考文献为:Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects。下载地址:ftp://ftp.genomics.org.cn/pub/gce

GCE 软件包中主要包含 kmer_freq_hash 和 gce 两支程序。前者用于进行 kmer 的频数统计,后者在前者的结果上进行基因组大小的准确估算。

2. GCE 下载和安装

$ wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
$ tar zxf gce-1.0.0.tar.gz -C /opt/biosoft

3. kmer_freq_hash 的使用

kmer_freq_hash 的常用例子:

$ /opt/biosoft/gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
  -k 21 -l reads.list -a 10 -d 10 -t 24 -i 50000000 -o 0 -p species &> kmer_freq.log

kmer_freq_hash 的常用参数:

-k <17>
    设置 kmer 的大小。该值为 9~27,默认值为 17 。
-l string
    list文本文件,其中每行为一个fastq文件的路径。
-t int
    使用的线程数,默认为 1 。
-i int
    初始的 hash 表大小,默认为 1048576。该值最好设置为 (kmer 的种类数 / 0.75)/ 线程数。如果基因组大小为 100M,测序了 40M 个 reads,reads 的长度为 100bp,测序错误率为 1%,kmer的大小为 21,则kmer的种类数为100M+40M*100*1%*21=940M,若使用24线程,则该参数设置为 i=940M/0.75/24=52222222。
-p string
    设置输出文件的前缀。
-o int
    是否输出 k-mer 序列。1: yes, 0: no,默认为 1 。推荐选 0 以节约运行时间。
-q int
    设置fastq文件的phred格式,默认为 64。该值可以为 33 或 63。
-c double
    设置k-mer最小的精度,该值位于 0~0.99,或为 -1。 -1 表示不对 kmer进行过滤。设置较高的精度,可以用于过滤低质量 kmer。精度是由 phred 格式的碱基质量计算得来的。
-r int
    设置获取 k-mer 使用到的 reads 长度。默认使用 reads 的全长。
-a int
    忽略read前面该长度的碱基。
-d int
    忽略read后面该长度的碱基。
-g int
    设置使用该数目的碱基来获取 k-mers,默认是使用所有的碱基来获取 k-mer。

kmer_freq_hash 的主要结果文件为 species.freq.stat。该文件有 2 列:第1列是kmer重复的次数,第二列是kmer的种类数。该文件有255行,第225行表示kmer重复次数>=255的kmer的总的种类数。该文件作为 gce 的输入文件。
kmer_freq_hash 的输出到屏幕上的信息结果保存到文件 kmer_freq.log 文件中。该文件中有粗略估计基因组的大小。其中的 Kmer_individual_num 数据作为 gce 的输入参数。

4. gce 的使用

gce 的常用例子:

$ /opt/biosoft/gce-1.0.0/gce \
  -f species.freq.stat -c 85 -g 4112118028 -m 1 -D 8 -b 1 > species.table 2> species.log

gce 的常用参数:

-f string
    kmer depth frequency file
-c int
    kmer depth frequency 的主峰对应的 depth。gce 会在该值附近找主峰。
-g int
    总共的 kmer 数。一定要设定该值,否则 gce 会直接使用 -f 指定的文件计算 kmer 的总数。由于默认下该文件中最大的 depth 为 255,因此,软件自己计算的值比真实的值偏小。同时注意该值包含低覆盖度的 kmer。
-M int
    支持最大的 depth 值,默认为 256 。
-m int
    估算模型的选择,离散型(0),连续型(1)。默认为 0,对真实数据推荐选择 1 。
-D int
    precision of expect value,默认为 1。如果选择了 -m 1,推荐设置该值为 8。
-H int
    使用杂合模式(1),不使用杂合模式(0)。默认值为 0。只有明显存在杂合峰的时候,才选择该值为 1 。
-b int
    数据是(1)否(0)有 bias。当 K > 19时,需要设置 -b 1 。

gce 的结果文件为 species.table 和 species.log 。species.log 文件中的主要内容:

raw_peak	now_node	low_kmer	now_kmer	cvg	genome_size	a[1]	b[1]
84	35834245	22073804	4044916750	84.6637	4.83093e+07	0.928318	0.637648

raw_peak: 覆盖度为 84 的 kmer 的种类数最多,为主峰。
now_node: kmer的种类数。
low_kmer: 低覆盖度的 kmer 数。
now_kmer: 去除低覆盖度的 kmer 数,此值 = (-g 参数指定的总 kmer 数) - low_kmer 。
cvg:估算出的平均覆盖度
genome_size:基因组大小,该值 = now_kmer / cvg 。
a[1]: 在基因组上仅出现 1 次的 kmer 之 种类数比例。
b[1]: 在基因组上仅出现 1 次的 kmer 之 数量比例。该值代表着基因组上拷贝数为 1 的序列比例。

如果使用 -H 1 参数,则会得额外得到如下信息:

for hybrid: a[1/2]=0.223671 a1=0.49108
kmer-species heterozygous ratio is about 0.125918

上面结果中,0.125918 是由 a[1/2] 计算出来的。 0.125918 = a[1/2] / ( 2- a[1/2] ) 。
a[1/2]=0.223671 表示在所有的 uniqe kmer 中,有 0.223671 比例的 kmer 属于杂合 kmer 。

此外,有 a[1/2] 和 b[1/2] 的值在最后的统计结果中。重复序列的含量 = 1 - b[1/2] - b[1] 。

则杂合率 = 0.125918 / kmer_size 。 若计算出的杂合率低于 0.2%,个人认为测序数据应该是纯合的。这时候,应该不使用 -H 1 参数。使用 -H 1 参数会对基因组的大小和重复序列含量估算造成影响。

5. 不同杂合率,有无重复序列的 kmer species 和 kmer individuals 图

下图中 a 和 b 是对理想中无重复的基因组在不同杂合率下的曲线图;
下图中 c 和 d 是对有重复的基因组(human)在不同杂合率下的曲线图。
从下图可以参考不同杂合率下的曲线状况。
kerm_pictures

使用FGAP进行补洞

1. FGAP简介

FGAP利用BLAST将contigs序列比对到基因组草图序列上,寻找重叠到gap区间的最优序列,从而进行补洞。其参考文献:Piro, Vitor C., et al. “FGAP: an automated gap closing tool.” BMC research notes 7.1 (2014): 371.

2. FGAP下载和安装

FGAP官网:http://www.bioinfo.ufpr.br/fgap/

$ wget http://sourceforge.net/projects/fgap/files/MCR_LINUX64b.tar.gz/download
$ tar zxf MCR_LINUX64b.tar.gz
$ cd MCR_LINUX64b
$ ./installMCR.sh /opt/biosoft/MCR

$ wget http://sourceforge.net/projects/fgap/files/FGAP_1_7_LINUX64b.tar.gz/download
$ tar zxf FGAP_1_7_LINUX64b.tar.gz -C /opt/biosoft/

2. FGAP的使用

FGAP的简单使用示例:

$ ln -s /opt/biosoft/FGAP_1_7_LINUX64b/sample_data/* .
$ /opt/biosoft/FGAP_1_7_LINUX64b/run_fgap.sh /opt/biosoft/MCR/v717/ \
 -d DRAFT_ecoli_hiseq454.fasta \
 -a "DATASET_ecoli_454.fasta,DATASET_ecoli_hiseq.fasta"

FGAP的使用参数:

-d /--draft-file	Draft genome file [fasta format - Ex: "draft.fasta"]
-a /--datasets-files	List of datasets files to close gaps [fasta format - Ex: "dataset1.fasta,dataset2.fasta"]

-s /--min-score		Min Score (raw) to return results from BLAST (integer) - Default: 25
-e /--max-evalue	Max E-Value to return results from BLAST (float) - Default: 1e-7
-i /--min-identity	Min identity (%) to return results from BLAST (integer [0-100]) - Default: 70

-C /--contig-end-length	Length (bp) of contig ends to perform BLAST alignment (integer) - Default: 300
-T /--edge-trim-length	Length of ignored bases (bp) upstream and downstrem of the gap (integer) - Default: 0
-R /--max-remove-length	Max number of bases (bp) that can be removed (integer) - Default: 500
-I /--max-insert-length	Max number of bases (bp) that can be inserted (integer) - Default: 500

-p /--positive-gap	Enable closing of positive gaps (with insertion) (integer [0-1]) - Default: 1
-z /--zero-gap		Enable closing of zero gaps (without insert any base) (integer [0-1]) - Default: 0
-g /--negative-gap	Enable closing of negative gaps (overlapping contig ends) (integer [0-1]) - Default: 0

-c /--gap-char				Base that represents the gap (char) - Default: "N"
-b /--blast-path			Blast+ package path (only makeblastdb and blastn are needed, version 2.2.28+ or higher) - Default: ""
-l /--blast-alignment-parameters	BLAST alignment parameters (opengap,extendgap,match,mismatch,wordsize) - Default: "1,1,1,-3,15"
-r /--blast-max-results			Max results from BLAST for each query (integer) - Default: 200
-t /--threads				Number of threads (integer) - Default: 1

-m /--more-output	More output files with gap regions after and before gap closing (integer [0-1]) - Default: 0
-o /--output-prefix	Output prefix [File or folder - Ex: "out" or "out/" ] - Default: "output_fgap"
-h /--help		This help message

使用ABACAS进行有参考基因组的contigs连接

1. ABACAS简介

ABACAS 用于在有参考序列的情况下,对contigs序列进行连接。软件官网:http://abacas.sourceforge.net/

2. ABACAS 下载和安装

ABACAS是一支perl程序,其运行需要安装有MUMmer。

$ wget http://sourceforge.net/projects/abacas/files/abacas.1.3.1.pl
$ chmod 755 abacas.1.3.1.pl
$ mv abacas.1.3.1.pl ~/bin/
$ ln -s ~/bin/abacas.1.3.1.pl ~/bin/abacas.pl
$ perl -p -i -e 's#/usr/local/bin/perl#/usr/bin/perl#' ~/bin/abacas.1.3.1.pl

3. ABACAS 的使用

ABACAS 的简单例子:

$ abacas.pl -r ref.fasta -q query.fasta -p nucmer -o out_prefix

ABACAS 的使用参数:

-r string
    输入参考序列。参考序列为fasta格式并仅有1条序列。
-q string
    输入query序列。query序列为fasta格式。
-p nucmer|promer
    该参数可以设为nucmer或promer,表示使用相应的程序进行contig ordering。
-o string
    设定输出文件的前缀。默认为query.fasta_ref.fastq。
-c
    参考序列是一个环。
-m
    使用该参数后,则将定位并定向的query序列按顺序打印到文件 out_prefix.MULTIFASTA.fa 中。
-b
    使用该参数后,则将未能定位的 congtig 序列内容打印到文件 ou_prefix.contigsInbin.fas 中。
-N
    使用该参数后,则将由定位定向的conitg序列连接而成的pseudomolecule的序列打印到文件 out_prefix.NoNs.fasta 中。默认下,输出文件 out_prefix.fasta 中的 pseudomolecule 序列中 congtig 之间的 gap 使用 N 表示,有重叠则加 100个N;该参数则是生成额外的fasta文件,其中contig之间序列无N。
-t
    使用该参数后,则对 ou_prefix.contigsInbin.fas 中的序列运行 blastall,额外生成 blast.out 文件。使用该参数则必须有 -b 参数。
-g out_prefix
    使用该参数后,则将 pseudomolecule 的gap部分在ref上的序列答应到文件 out_prefix.Gaps_onRef 中。
-P
    使用该参数后,则设计用于 close gaps 的引物。需要安装 prime3_core,但我个人加该参数运行失败。
-R
    使用该参数,则不需要运行 mumer,适合于多次运行程序,节约时间。但我个人加该参数运行失败。
-e
    不进行 contig ordering,有利于直接进行引物设计。
-i int
    允许的最小的identity的百分比。默认为 40。
-v int
    允许的最小的 contig 覆盖率。默认为 40。
-V 0|1
    contig比对到ref序列上的位置的数目。默认为1。使用 -V 0 表示将contig序列随机放到1个比对上的位点。
-l int
    不使用低于此长度的contig序列。默认为 100。 
-d
    ABACAS在调用MUmer的时候,默认下使用了--maxmatch参数来增加匹配,使用 -d 参数则会使用MUmer的默认参数,即不会使用--maxmatch参数。

使用 REAPR 进行基因组错误评估

1. REAPR 简介

REAPR(Recognition of Errors in Assemblies using Paired Reads)能利用成对的reads来识别基因组序列中的错误。从而,能将基因组序列从错误的gap处断开或将错误序列使用 Ns 代替。同时,对错误信息进行统计。

2. REAPR下载和安装

REAPR官网:http://www.sanger.ac.uk/resources/software/reapr/
安装 REAPR 需要先安装 R 和 Perl 模块: File::Basename, File::Copy, File::Spec, File::Spec::Link, Getopt::Long, List::Util。

$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.tar.gz
$ tar zxf Reapr_1.0.17.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/Reapr_1.0.17
$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.manual.pdf
$ ./install.sh
$ echo 'PATH=$PATH:/opt/biosoft/Reapr_1.0.17' >> ~/.bashrc
$ source ~/.bashrc

测试Reapr能否正常运行:
$ wget ftp://ftp.sanger.ac.uk/pub4/resources/software/reapr/Reapr_1.0.17.test_data.tar.gz
$ tar zxf Reapr_1.0.17.test_data.tar.gz 
$ ./test.sh

3. REAPR 使用

reapr 最少需要输入基因组fasta文件和mate-pair数据(需要将RF方向转变成FR方向)文件。常用的运行步骤如下:

使用 facheck 命令检查基因组fasta文件,以防fasta文件序列名含有怪异字符。
$ reapr facheck genome.fasta
使用下列命令对fasta文件进行修改。
$ reapr facheck genome.fasta new_genome.fasta

使用SMALT将 mate-pair reads比对到基因组上。
$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam

运行reapr pipeline对基因组进行评估。将结果输出到 outdir 文件夹中。
$ reapr pipeline genome.fasta jumping.bam outdir

reapr 也支持输入paired-end数据(可选)。通过将paired-end数据uniquely并perfecly比对到基因组,从而能计算出基因组上error-free的碱基位点并进行统计。该项目分析对error calling没有影响。值得注意的是: 默认设置下call一个error-free位点需要至少5x的覆盖度;同时输入的paired-end数据的方向要为inward(一般paired-end数据方向即为inward),此外reapr仅支持inward方向,因此输入的jumping数据也要先转换为inward方向。

使用paired-end数据的命令行:

$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam
$ reapr perfectmap genome.fasta frag.1.fastq frag.2.fastq insert_size perfect_prefix
$ reapr pipeline genome.fasta jumping.bam outdir perfect_prefix

reapr的运行流程:
reapr_pipeline

4. reapr的结果

REAPR 会对基因组每个位点的 FCD(fragment coverage distribution)进行分析。期望的FCD和实际上的上FCD之差为”FCD error”。 FCD error往往代表不正确的scaffold连接、大的插入缺失或contig序列中不正确的连接。为了能更好地计算FCD,特别是gap位点的FCD,mate-pair数据的insert size越大越好。

REAPR对输入的paired reads数据进行检测,并将结果放入到文件夹 outdir/00.Sample/ 下。

gc_vs_cov.lowess.pdf: GC偏好性检测图
insert.in.pdf: FR(inner/inward)类型插入片段的长度检测图
insert.out.pdf: RF类型插入片段的长度检测图
insert.stats.txt: 插入片段长度统计

从错误位点打断后的基因组序列: outdir/04.break.broken assembly.fa 。
如果FCD error区有gap,则从gap处将基因组序列打断;若FCD error区没有gap,则将其用Ns代替。同时将被代替的序列写入到文件 outdir/04.break.broken assembly bin.fa中。

基因组错误统计文件: outdir/05.summary.report.txt。该文件统计位于 contig 中的 FCD error 的数目,和含有 gap 的 FCD error 的数目。