Installing QIIME-1.9.1 on CentOS 6.5 (By Yue Zheng)

此方法是由郑越同学提供的。

QIIME consists of native Python 2 code and additionally wraps many external applications. As a consequence of this pipeline architecture, QIIME has a lot of dependencies and can be very challenging to install.

1. Setting up qiime-deploy on CentOS

1.1 sudo vim /etc/yum.repos.d/zeromq.repo

Paste the following into that file:

[home_fengshuo_zeromq]
name=The latest stable of zeromq builds (CentOS_CentOS-6)
type=rpm-md
baseurl=http://download.opensuse.org/repositories/home:/fengshuo:/zeromq/CentOS_CentOS-6/
gpgcheck=1
gpgkey=http://download.opensuse.org/repositories/home:/fengshuo:/zeromq/CentOS_CentOS-6/repodata/repomd.xml.key
enabled=1

Save and exit that file

1.2 Install the qiime-deploy dependencies on your machine

sudo yum groupinstall -y "development tools"
sudo yum install -y ant compat-gcc-34-g77 java-1.6.0-openjdk java-1.6.0-openjdk-devel freetype freetype-devel zlib-devel mpich2 readline-devel zeromq zeromq-devel gsl gsl-devel libxslt libpng libpng-devel libgfortran mysql mysql-devel libXt libXt-devel libX11-devel mpich2 mpich2-devel libxml2 xorg-x11-server-Xorg dejavu* python-devel sqlite-devel tcl-devel tk-devel R R-devel ghc

2. Installing requisite Python and R packages

# Installing sqlite-devel
sudo yum install sqlite-devel –y

# Installing Python 2.7
wget https://www.python.org/ftp/python/2.7.8/Python-2.7.8.tgz
tar xf Python-2.7.8.tgz
cd Python-2.7.8
./configure --prefix=/usr
make && make install

# Install setuptools & pip
# First get the setup script for Setuptools:
wget https://bitbucket.org/pypa/setuptools/raw/bootstrap/ez_setup.py
# Then install it for Python 2.7 :
sudo python2.7 ez_setup.py
# Now install pip using the newly installed setuptools:
sudo easy_install-2.7 pip
# With pip installed you can now do things like this:
pip2.7 install [packagename]

# Install virtualenv for Python 2.7
sudo pip2.7 install virtualenv

# Check the system Python interpreter version
python --version
# This will show Python 2.7.8

# Maybe you will found yum can not be used this moment, because yum is associated with python2.6. Thus, we modified the yum conf files to use python2.6
sudo vim /usr/bin/yum
# Replace “#!/usr/bin/python” by “#!/usr/bin/python2.6”
# Installing R packages
# Run R and execute the following commands
install.packages(c('ape', 'biom', 'optparse', 'RColorBrewer', 'randomForest', 'vegan'))
source('http://bioconductor.org/biocLite.R')
biocLite(c('DESeq2', 'metagenomeSeq'))
q()

3. Install the latest QIIME release and its base dependencies is with pip

sudo pip2.7 install numpy
sudo pip2.7 install qiime -v
# For Chines user, you may find the suspend of pip, as the limitation of network. For example, If FastTree cannot be download, you can download it by another port of internet, and then post the install package into your local address. Next step, Downloading the qiime-1.9.1.tar.gz and changing the description of FastTree in setpu.py. After you modified the qiime-1.9.1.tar.gz you can post it into your local address. Finally, run sudo pip2.7 install qiime –v –i [local address]

# Installing QIIME 1.9.0's dependencies
# Downloading the zip packages of ‘qiime deploy’ and ‘qiime deploy conf’ from Github
cd
unzip qiime-deploy-master.zip qiime-deploy-conf-master.zip
mkdir ~/qiime_software
cd qiime-deploy-master
sudo python2.7 qiime-deploy.py ~/qiime_software/ -f ~/qiime-deploy-conf/qiime/qiime-1.9.1/qiime.conf --force-remove-failed-dirs
# After this step, it will display the list including ‘Packages deployed successfully’, ‘Packages skipped’ and ‘Packages failed to deply’
source ~/qiime_software/active.sh
print_qiime_config.py –tf

# If there are some packages were uninstalled, you should install them manually
# For example, usearch and amplicannoise were failed to install.

# Installing usearch manually
# Visting http://www.drive5.com/usearch/download.html to download the USEARCH v5.2.236
# Moving the binary file into /usr/bin and change the name as usearch, then chmod 755 [the binary file] 

# Installing usearch manually
# Downloading the AmpliconNoiseV1.27.tar.gz
tar -xvzf AmpliconNoiseV1.27.tar.gz
cd AmpliconNoiseV1.27
make clean
make
make install
echo "export PATH=$HOME/AmpliconNoiseV1.27/Scripts:$HOME/AmpliconNoiseV1.27/bin:$PATH" >> $HOME/.bashrc
echo "export PYRO_LOOKUP_FILE=$HOME/AmpliconNoiseV1.27/Data/LookUp_E123.dat" >> $HOME/.bashrc
echo "export SEQ_LOOKUP_FILE=$HOME/AmpliconNoiseV1.27/Data/Tran.dat" >> $HOME/.bashrc

# PATH Environment Variable
echo "export PATH=$HOME/bin/:$PATH" >> $HOME/.bashrc
source $HOME/.bashrc

# Finnaly verification
source ~/qiime_software/active.sh
print_qiime_config.py –tf

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)

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

使用 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 的数目。

使用SOPRA进行scaffold连接

1. SOPRA 简介

SOPRA(Statistical Optimization of Paired Read Assembly),主要用于利用成对的reads信息来进行scaffold连接,其准确性非产高。其参考文献:Dayarian, Adel, Todd P. Michael, and Anirvan M. Sengupta. “SOPRA: Scaffolding algorithm for paired reads via statistical optimization.” BMC bioinformatics 11.1 (2010): 345.

2. SOPRA的下载和安装

$ wget http://www.physics.rutgers.edu/~anirvans/SOPRA/SOPRA_v1.4.6.zip
$ mkdir /opt/biosoft/SOPRA_v1.4.6/
$ unzip SOPRA_v1.4.6.zip -d /opt/biosoft/SOPRA_v1.4.6/
安装目录/opt/biosoft/SOPRA_v1.4.6/下有SOPRA的manual文档
$ chmod 755 /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
$ perl -p -i -e 's\s*$/\n/' /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/*.pl
SOPRA的主要运行程序就是这4支perl程序。

3. SOPRA的使用流程

SOPRA 可以支持 Illumina 和 SOLID 的二代测序数据,但此处仅讲解使用 Illumina 的测序数据。SOPRA 也支持同时利用多个paired测序文库的数据,同时支持paired-end和mate-paired数据。

3.1 准备输入数据

不管是paired-end或mate-paired数据,首先需要将数据的方向转换为FR方向。即需要将mate-paired的2端reads都进行反向互补。然后,需要将fatstq文件2端的reads文件合并为一个文件,并转换为fasta文件作为SOPRA的输入数据。
此外,SOPRA还需要contig序列作为输入的基因组序列。
以1个paired-end数据和1个mate-paired数据为例:

先将 jumping 文库的数据进行反向互补
$ fastq_rc.pl jump.1.fastq > 11; mv 11 jump.1.fastq
$ fastq_rc.pl jump.2.fastq > 22; mv 22 jump.2.fastq

然后,合并两端的reads序列
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  frag.1.fastq frag.2.fastq frag.fastq
$ /opt/biosoft/velvet_1.2.10/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  jump.1.fastq jump.1.fastq jump.fastq

再将fastq转换为fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' frag.fastq > frag.fasta
$ perl -e '$num = 1; while (<>) { print ">$num\n"; $_=<>; print; $_=<>; $_=<>; $num ++; }' jump.fastq > jump.fasta

3.2 SOPRA 对contigs序列和测序数据进行格式处理

使用s_prep_contigAseq_v1.4.6.pl程序对基因组的contig序列以及illunia测序数据进行序列名的修改,以利于后续程序运行。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_prep_contigAseq_v1.4.6.pl \
  -contig contig.fasta -mate frag.fasta jump.fasta -a SOPRA_OUT/

程序参数:
-contig string
    输入基因组 contigs 序列。
-mate string
    输入 illumina 数据。如果有多个数据,则多个fasta文件之间使用空格分开。
-a string
    设置输出文件夹目录。

在 SOPRA_OUT 目录中生成了 contigs_sopra.fasta、frag_sopra.fasta 和 jump_sopra.fasta 文件。查看这3个文件内容,发现是对fasta文件中的序列名进行了修改。

3.3 将Illumina数据比对到基因组

可以利用bowite,bwa等序列比对软件将Illumina数据按单端序列比对到基因组序列上。若一条序列比对到多个位置,推荐报告多个结果。因为后续的s_parse_sam_v1.4.6.pl程序会去除掉比对到多个位置的reads。若仅报告最优结果,我不知道程序会如何处理。要是不去除这样的比对结果,可能会对scaffolding的准确性有影响。因此,还是要注意报告多个比对结果,特别是bowtie2默认下仅报告最优结果。
此外,去除reads尾部碱基质量低的碱基后,illumina reads 的匹配率若是提高很多,则推荐去除reads尾部10~20个碱基,再用于比对。由于mate-pair数据进行了反向互补,则其原本的尾部成了首部。

$ cd SOPRA_OUT/
$ bowtie2-build contigs_sopra.fasta contig
$ bowtie2 -p 20 -x contig -k 10 -f -3 15 -U frag_sopra.fasta -S frag_sopra.sam
$ bowtie2 -p 20 -x contig -k 10 -f -5 15 -U jump_sopra.fasta -S jump_sopra.sam

程序参数:
-k 10
    若read比对到多个位点,则最多报告 10 个结果。
-3 15
    去除read尾部15bp碱基后再用于比对。
-5 15
    去除read首部15bp碱基后再用于比对。

3.4 对SAM文件进行分析

使用s_parse_sam_v1.4.6.pl对sam文件进行分析并简化其信息。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_parse_sam_v1.4.6.pl \
  -sam frag_sopra.sam jump_sopra.sam -a ./

程序参数:
-sam string
    输入 sam 文件。多个illumina 文库的sam文件用空格隔开。
-a string
    设置输出文件的路径。

得到sam文件分析后的结果frag_sopra.sam_parsed和jump_sopra.sam_parsed。

3.5 进行序列方向和距离分析

读取上一步简化后的sam文件信息,进行分析。得到contigs序列的覆盖度、library文库的插入片段长度,序列长度等。

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_read_parsed_sam_v1.4.6.pl \
  -parsed frag_sopra.sam_parsed -d 500 -parsed jump_sopra.sam_parsed -d 3000 -a ./

程序参数:
-parsed string
    输入 .sam_parsed 文件。若有多个文件,则用空格分割。
-d int
    illumina数据的插入片段长度。
-c int default: 5
    如果read及其反向互补序列在library中出现的次数>=该值,则不计算其成对信息。
    在运行 s_scaf_v1.4.6.pl 程序后,输入日志中得到如下信息:
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100,则需要考虑降低 -c 参数的值。
    此外,如果数据量比较少,平均覆盖度较低,比如低于 10,则考虑增加 -c 参数的值。

-e int default:0
    如果设置其值为 1, 则使用输入的插入片段长度值。默认下软件会计算插入片段长度值,并使用软件计算出来的值。
-a string
    设置输出文件的路径。

输出文件为 orientdistinfo_c5 。 此外,在 div 文件夹中有contigs序列的覆盖度、library文库的插入片段长度,序列长度等统计信息。

3.6 进行scaffold连接

$ /opt/biosoft/SOPRA_v1.4.6/source_codes_v1.4.6/SOPRA_with_prebuilt_contigs/s_scaf_v1.4.6.pl \
  -o orientdistinfo_c5 -a ./

程序常用参数:
-o string
    输入文件 orientdistinfo_c。
-w int default: 4
    连接2个contigs所需要的最小的pair links。
    在运行 s_scaf_v1.4.6.pl 程序后,输入日志中得到如下信息:
    Average number of links between two contigs using minlength 150 and minlink 2 is 63.38, ...
    Starting cycle 1 of orientation assignment
    Average number of links between two contigs using minlength 150 and minlink 4 is 184.84, ...
    如果该值(184.84)大于100,则需要考虑提高 -w 参数的值。一般 -w 取该值的 4~5%。如果该值较低,比如低于10的时候,则需要降低 -w 参数的值。
-L int default: 150
    用于连接scaffold的最小长度的contigs序列。
-h float default: 2.2
    设置一个系数,用于排除对高覆盖度contigs的scaffold连接。大于 mean_coverage + h * std_coverage 该覆盖度的 contigs ,不能用于scaffold连接
-a string
    设置输出文件的路径。

程序输出文件为 scaffolds_h2.2_L150_w4.fasta 。该文件为最终的结果文件。

使用 snoscan 和 snoGPS 进行 snoRNAs 分析

1. snoscan 和 snoGPS 简介

snoRNAs 主要有 2 种不同类型。其中 C/D box 类型 snoRNAs 用于指导甲基化修饰(2′O-ribose-methylations),该甲基化修饰位于核糖体 2′-OH上; H/ACA box 类型 snoRNAs 用于指导假尿苷酸化修饰(Pseudouridylation),即将尿嘧啶转换为假尿嘧啶。

snoscan用于鉴定 C/D box 类型 snoRNAs; snoGPS用于鉴定 H/ACA box 类型的 snoRNAs。

可以使用Snoscan Server 1.0 网页工具snoGPS Server 1.0 网页工具进行 snoRNAs 预测。网页工具支持的最长序列长度为 100kb 。

2. snoscan 与 snoGPS 下载与安装

snoscan 的下载和安装:

$ wget lowelab.ucsc.edu/software/snoscan.tar.gz
$ tar zxf snoscan.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/snoscan-0.9b/squid-1.5j
$ perl -p -i -e 's/getline/get_line/g' sqio.c
$ cd ..
$ make
$ echo 'PATH=$PATH:/opt/biosoft/snoscan-0.9b' >> ~/.bashrc
$ perl -p -i -e 's#/usr/local/bin/perl#/usr/bin/perl#' sort-snos

snoGPS 的下载和安装:

$ wget http://lowelab.ucsc.edu/software/snoGPS-0.2.tar.gz
$ tar zxf snoGPS-0.2.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/snoGPS-0.2/src
$ perl -p -i -e 's/getline/get_line/g' lib/*.c squid/lib/*.c inc/*.h
$ make
$ cd ..
$ ln -s src/pseudoU_test snoGPS
$ echo 'PATH=$PATH:/opt/biosoft/snoGPS-0.2' >> ~/.bashrc 
$ source ~/.bashrc
$ perl -p -i -e 's#/usr/local/bin/perl#/usr/bin/perl#' /opt/biosoft/snoGPS-0.2/scripts/*.pl
$ echo 'export PERL5LIB=$PERL5LIB:/opt/biosoft/snoGPS-0.2/perlModules/' >> ~/.bashrc
$ sudo cpan -i Class::MethodMaker File::Sort Bio::Tools::Run::Alignment::TCoffee

3. snoscan 的使用

snoscan的常用例子:

$ snoscan rRNA.fa genome.fasta -s -o out

rRNA.fa      : rRNA序列的fasta格式
genome.fasta : 用于搜索 snoRNAs 的基因组序列

snoscan的常用参数:

-h
    打印帮助信息
-m string
    指定包含甲基化位点信息的文件。该文件的格式为:
        >seq_A
        T 5 Verified meth site by TML
        C 12  Predicted by snoscan
        >My-seq-B
        G 1   Partially modified site
        A 6   Mapped by Maden
    与fasta格式类似,头部是和 genome.fasta 序列名相同;内容部分分3列,以空格或制表符分割,第1列是位点的碱基,第2列是1-based的位点,第3列是60个字符以内的描述。
-o string
    将结构输出到指定的文件,默认输出到标准输出
-s
    是否保存 snoRNA 的序列到输出中

4. snoGPS 的使用

4.1 获得可能的假尿苷酸化位点序列

首先,提取 rRNA 序列信息中可能发生假尿苷酸化的位点。一般情况下,认为所有碱基为 U 的位点都可能发生假尿苷酸化修饰。通过程序将 碱基U及其侧翼10bp序列提取出来,提取出来,得到 target 文件。

$ /opt/biosoft/snoGPS-0.2/scripts/getRnaTargets.pl -n rRNA.fa 

-n
    将所有的位点信息写入到一个文件中。
-r
    制作反向 RNA 序列的 targets。
-f string
    输入含有目标序列假尿苷酸化修饰的位点信息的文件。该文件格式例如:
        rRNA1 species 43:46:53
        rRNA2 species 6:7:15:34:37:39:41:43:44:54:58:60:69:72:89:91
    文件内容为用空格分割的3列。第1列是rRNA序列名;第2列是一个单词的描述;第3列是假尿苷酸位点。    
-s string
    比如输入 '656:989:2010',则对指定位点的尿苷酸制作 target 文件。默认下是对所有尿苷酸位点制作 target 文件。

4.2 运行snoGPS

snoGPS 的常用例子:

$ snoGPS -D 0 -S 15 -T snoGPS_tmp/uridin_target_file_for_snoGPS.txt \
  -F snoGPS_hits.fa -L /opt/biosoft/snoGPS-0.2/scoretables/haca2stemv7.tables \
  genome.fasta /opt/biosoft/snoGPS-0.2/desc/haca2stemv4a.desc > snoGPS_out

genome.fasta : 在该基因组序列中鉴定 snoRNAs。
haca2stemv4a.desc : 该文件为 desc 文件,用于决定一些参数。对于human和yeast有不同的desc文件。
程序将主要结果输出到标准输出。

snoGPS 的常用惨素:

-h
    打印帮助文档
-S float
    设置得分阈值。
-T string
    设置 target 文件,即含有可能的假尿苷酸化修饰位点信息的文件,由rRNA序列信息获得。
-F string
    输出的 hits 文件。其中为 snoRNAs 序列。
-t int
    设置从 target 文件中读取的 target 位点数目。
-L string
    载入 score-tables 文件。snoGPS 提供了2个文件,分别适合于 human 和 yeast。该文件位于 /opt/biosoft/snoGPS-0.2/scoretables/ 目录下。
-W
    仅仅对正链进行检测
-C
    仅仅对负链进行检测

5. snoRNA 的一些知识点

Small nucleolar RNA 的 wikipedia

snoRNA 指导的 rRNA 的转录后修饰有 2 种: 甲基化(Methylation) 和 假尿苷酸化(Pseudouridylation)。

snoRNP(small nucleolar ribonucleoprotein): 每个 snoRNA 和至少4个蛋白分子形成复合物,来指导目标RNA的修饰。

snoRNAs 主要有 2 种不同类型。其中 C/D box 类型用于指导甲基化修饰; H/ACA box 类型用于指导假尿苷酸化修饰。

C/D box 类型的 snoRNAs 包含 C(RUGAUGA)环 和 D(CUGA)环,分别在 snoRNA 的 5′ 和 3′ 端。如下图所示:
C/D box

H/ACA box 类型的 snoRNAs 由 2 个发夹结构和 2 个单链区组成,称为 hairpin-hinge-hairpin-tail 结构,如下图所示。结构中包含 H box(ANANNA) 和 ACA box(ACA)。
H/ACA box