由于之前的1年将文章全都写入到其它地方。而现在我的博客完全转移到了阿里云平台,因此更加稳定了。现在将博客内容转移到过来。
日度归档:2015 年 8 月 25 日
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
结果如下:
使用 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)在不同杂合率下的曲线图。
从下图可以参考不同杂合率下的曲线状况。
使用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参数。