创建 BioProject 号和 BioSample 号

对某一个物种进行了基因组测序，则申请 BioProject 和 BioSample 号各一个。

使用 tbl2asn 准备后缀为 .sqn 的 ASN.1 文件

在 windows 下可以使用 Sequin 来制作 .sqn 文件。该文件是下面所述的 3 个文件的信息的综合体。tbl2asn 是命令行的工具，适合大基因组数据的 .sqn 文件生成。

1. 生成包含作者信息的 .sbt 模板文件(Submission Template)

推荐使用网页http://www.ncbi.nlm.nih.gov/WebSub/template.cgi，填入数据生成 template.sbt 文件，并下载到本地。当然，此文件也可以使用 Sequin 生成。
填写信息时，可填入 BioProject 和 BioSample 号。

2. 准备后缀为 .fsa 的fasta文件

fasta 文件的的 header 要求如下：

1. ">" 和 第一个空格之间的内容是序列名。
2. header部分可以加入其它因素，比如：
organism [organism=Saccharomyces cerevisiae]
strain [strain=S288C]
isolate [isolate=CWS1]  # 代表在什么个体上获得的样品
chromosome [chromosome=XVI]
topology [topology=circular]
location [location=mitochondrion]
molecule [moltype=mRNA] (DNA is the default)
technique [tech=wgs]
protein name [protein=helicase]
genetic code [gcode=4]

3. 准备后缀为 .tbl 的表格格式的基因组注释信息文件

此文件有 5 列，每列用 tab 分割，称为 feature table。
此文件是最为关键的一步。该文件必须包含：编码基因的结构注释信息、非编码基因的结构注释信息 和 基因的功能注释信息。一旦做不好，NCBI的工作人员就会发email反馈修改意见。

feature table 格式的要点如下：

1. 对每条序列的所有注释之前，有一行额外的内容，例如：
>Feature scaffold_1
该行内容后面的所有注释信息属于序列 scaffold_1 ，一定不能遗漏 Feature 这个单词，Feature 和 scaffold_1 用空格分隔。
2. 每个 feature 使用 5 行内容进行阐述，并分成 2 个部分。
第 1 部分是 feature 在序列上的结构信息。有 3 列，分别为该 feature 的起始位点、结束位点和 feature 名。若 feature 在正义链上，则起始位点 < 结束位点，若在负义链上，则起始位点 > 结束位点。若 feature 为断裂基因的 CDS 或 exon 等信息时，则有多行数据，但仅在其首行的第 3 列上显示 feature 名。
第 2 部分是 feature 的功能注释信息。使用第 4、5 列，前面有 3 个 tab 键。第 4 列对应 feature 的 qualifier，第 5 列是 qualifier 的值。 qualifier 是对 feature 的描述标签。如果有多个 qualifier 及其值，则用多行进行表示。
3. feature 和 qualifier 的具体标签名称参考http://www.insdc.org/documents/feature_table.html。
4. 常用的 feature 名称有：gene, mRNA, CDS, exon, 5'UTR, 3'UTR, tRNA, rRNA, ncRNA 等。其中 ncRNA 是指除了 tRNA 和 rRNA 以外的其余 ncRNA。
5. gene 的 qualifier 标签一般是 gene, 第 5 列使用 NCBI 提供的 locus_tag + 数字编号。 mRNA 和 CDS 的 qualifier 标签一般使用 product，第 5 列是 Nr 注释的结果。exon 的 qualifier 标签一般使用 number，其值为 1,2,3... 。 UTR 的 qualifier 标签可以使用 note。 tRNA 的 qualifier 标签一般使用 product，第 5 列是 tRNA 名称, 例如 tRNA-Lys。rRNA 的 qualifier 标签要有 gene 和 product，第 5 列中product 是 "16S ribosomal RNA", "23S ribosomal RNA", "5S ribosomal RNA", 相应的 gene 的值可以是 rrsA, rrlA, rrfA ... , rr 表示是 rRNA， s l f 分别对应 16 23 5, A是一个编号，下面的编号是 B, C D... 此外，每个 rRNA 区域要有个 gene 的 feature。 ncRNA 的 qualifier 标签中必须有 ncRNA_class，第 5 列则是 ncRNA 的类别，比如 miRNA, siRNA, scRNA 等。此外，可以使用 note 作为 qualifier 的标签，其值可随意标示。
6. mRNA 和 CDS 的 product 的取值，使用 Nr 注释的最优结果。最优结果如果包含 "hypothetical protein" 、 "predicted protein" 、 "unknown" 、 "partial"  或 "homolog" 时，则需要取其它注释结果，或采取一定的措施了。

原核生物 feature table 的说明：http://www.ncbi.nlm.nih.gov/genbank/genomesubmit/#prepare_table。
原核生物的 feature table 的要点：

1. 必须包含的 feature 是 gene, CDS, rRNA 和 tRNA。不要有 mRNA 。并且，每个 CDS，rRNA 或 tRNA 都属于一个 gene。
2. Gene 的 qualifier 标签必须有 locus_tag， 也可以有 gene。 gene 的值为 gene 的名称，其名称有相应的标准，以 3 个小写字母开始的。
3. CDS 的 qualifier 标签必须有 product 和 protein_id 。product 的值也有相应的标准。protein_id 的值一般为 gnl|xxxx|string，其中 xxxx 推荐是实验室的名字， string 是 protein_id 的标示，可以使用 locus_tag 。
4. rRNA, tRNA, misc_RNA 和 ncRNA 都必须有相应的 gene feature， 其 qualifier 必须有 product 。与 RNA 相应的 gene 的 qualifier 中，其 gene 的值：5S rRNA => rrfA, 16s rRNA => rrsA, 23s rRNA => rrlA, tRNA-Lys => tRNAK, tRNA-Thr => tRNAT 。 rRNA 和 tRNA 的注释必须要有。

4.

4. tbl2asn 命令生成 .sqn 文件

在当前目录下生成了 3 个文件: species.sbt, species.fsa, specis.tbl 。
运行 tbl2asn 生成目标文件 species.sqn 。

tbl2asn -t C001.sbt -p ./ -a s -V vb
# -a s 一个fasta文件有多条序列时，使用此参数配置。
# -V vb v表示对输入的数据进行验证，生成 2 个 .val 的文件；-b 生成GeneBank格式的文本文件，以 .gbf 为后缀。
# 运行完毕后需要查看 val 文件，其中可能有很多错误与警示信息。 有些蛋白质序列不是以 M 开头，会在此处提示 ERROR 。特别是细菌基因组会出现这样的提示。应该在 .fsa 文件的 header 部分加上 [gcode=11] 来解决。表明其遗传密码子表是 11 号。根据错误，去除所有的错误提示。

使用 GenomeMacroSend 上传 .sqn 文件

在 GenomeMacroSend 网页 http://www.ncbi.nlm.nih.gov/projects/GenomeSubmit/genome_submit.cgi 的最下方的输入框中填写信息上传 .sqn 文件。

全网页方法上传数据

基因组数据上传：Genomes(WGS) submission portal
转录组数据上传：TSA submission portal
使用网页方式上传数据和上述方法基本一致。 feature tab 的制作依然需要自己手工制作，再上传。

使用 latex 命令将 TEX 文件转换为 DVI 文件
使用 dvipdf 命令将 DVI 文件转换为 pdf 文件
例如：

$ latex wise2.tex
$ dvipdf wise2.dvi

1. GeneWise 简介

Genewise主要用于将蛋白质序列和DNA序列进行比对，从而对DNA序列上的编码区进行预测。

2. GeneWise 安装

$ wget http://www.ebi.ac.uk/~birney/wise2/wise2.4.1.tar.gz
$ tar zxf wise2.4.1.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/src

$ yum install *glib*
$ find ./ -name makefile | xargs sed -i 's/glib-config/pkg-config --libs glib-2.0/'
$ export C_INCLUDE_PATH=/usr/include/glib-2.0/:/usr/lib64/glib-2.0/include/:$C_INCLUDE_PATH
$ perl -p -i -e 's/getline/get_line/g' ./HMMer2/sqio.c
$ perl -p -i -e 's/isnumber/isdigit/' models/phasemodel.c
$ make all
$ export WISECONFIGDIR=/opt/biosoft/wise2.4.1/wisecfg/
$ make test
$ echo 'PATH=$PATH:/opt/biosoft/wise2.4.1/src/bin/' >> ~/.bashrc
$ echo 'export WISECONFIGDIR=/opt/biosoft/wise2.4.1/wisecfg/' >> ~/.bashrc 
$ source ~/.bashrc

3. GeneWise的使用

在GeneWise的安装目录下，有一个wise2.pdf文件，阐述了详细的genewise的使用方法。其软件最常用的命令是genewise。该命令的常用示例：

genewise protein.fasta dna.fasta -both -gff

程序输入的蛋白质序列和DNA序列分别是2个fasta文件。这两个fasta文件中仅有第一条序列是有效的，genewise仅对其中的2个第一条序列进行比对。以上示例对dna序列的正负链都进行cds预测，并将gff格式结果文件输出到标准输出。

genewise的常用参数：

-trev
    仅对负义链进行cds预测。
-tfor
    仅对正义链进行cds预测。该参数是默认值。
-both
    对负链都进行cds预测。
-genes
    给出gene结构的结果，非常简单的exon信息结果。默认情况下仅输出适合人类阅读的比对结果。
-gff
    给出gff格式的结果。
-cdna
    给出cdna序列。
-pep
    给出cds翻译出的蛋白质序列。
-splice [model/flat]
    使用的split site是model（默认值）或GT/AG。
-help
    给出帮助信息。
-version
    给出版本信息。
-silent
    标准错误输出不输出messages信息。
-quiet
    标准错误输出不输出report/info信息。

genewise的运行原理简述：

1. genewise的算法：21:93算法是genewise的基础算法。该算法简单讲就是 Match-Insert-Delete，在蛋白质序列和DNA序列比对后能准确划定intron边界。算法将intron分成5部分：5'端splice site、中间intron主体、富含CT区域、连接区、3'端splice site。根据蛋白质序列和DNA序列的比对结果算出Intron部分，从而将DNA序列的CDS区分成了Match、Insert和Delete 3部分，再对这3部分进行蛋白质翻译或移码翻译，从而划定intron边界，得到CDS信息。
2. 6:23算法则是2:93算法的简单版本，也是软件的默认设置。和2:93算法相比，6:23算法的intron没有第3和第4部分（富含CT区域、连接区）。同时，6:23算法更适合于DNA序列中没有屏蔽重复或introns序列比较怪异的情况。使用该算法的时候，-intron参数的值得是tied（也是该参数默认的值），否则会得到错误的很长的intron结果。
3. 若是算法后面带个 L 字样，则表示适用于进行输入的蛋白质序列是 HMM 模型。此外， 还有其它的一些算法，可以参考wise2.pdf文件。
4. genewise对基因进行预测后，有一个得分。该得分 = log2(预测模型的可能性/随机结果的可能性) 。因此，0表示该结果是个随机的结果，不可靠的。根据软件作者的经验，得分高于35的结果是非常可靠的；得分25~35的结果是可信的；得分18~25的结果可能仅适用于某些蛋白质家族；得分低于15的是不可信的。

4. GeneWise的高级使用

用临近物种的protein序列对基因组进行homolog gene预测的时候，需要通过blast将proteins序列和基因组序列进行比对，再提取基因组的目标基因区域和最佳结果protein进行genewise分析。因此，需要自己写一些程序进行并行化的genewise计算，从而达到对全基因组大数据的分析。Genewise软件提供了一支程序/opt/biosoft/wise2.4.1/src/perl/scripts/blastwise.pl程序能进行该项处理（我没有用过该程序，我自己写了想要的程序）。

1. ICORN 简介

ICORN (Iterative Correction Of Reference Nucleotides), 能通过将 reads 比对到基因组，从而修正 SNP 和 小于 3bp 的 INDEL 位点。现在出了新的版本 ICORN2 。
ICORN 官网： http://icorn.sourceforge.net/
ICORN 参考文献： Otto T D, Sanders M, Berriman M, et al. Iterative Correction of Reference Nucleotides (iCORN) using second generation sequencing technology[J]. Bioinformatics, 2010, 26(14): 1704-1707.

2. ICORN 原理

ICORN 的原理如下图所示。ICORN2 调用 SMALT 将 reads map 到基因组序列上；然后 Call SNPs 和 <=3bp INDELs；根据其位点的覆盖度情况决定基因组上该位点的核苷酸类型。最后通过多轮迭代直到不能 call 到新的 variants 为止。

3. ICORN 的下载和安装

ICORN 的使用会调用第三方的软件： GATK、SMALT、samtools 和 SNP-o-AMTIC 。

$ wget ftp://ftp.sanger.ac.uk/pub4/resources/software/pagit/ICORN2/icorn2.V0.95.tgz
$ tar zxf icorn2.V0.95.tgz -C /opt/biosoft/
$ /opt/biosoft/ICORN2/icorn2.sh --help

4. ICORN 的使用

4.1 使用前准备

准备的输入文件是： readroot_1.fastq、readroot_2.fastq、genome.fasta，将这3个文件放置于工作目录中。
同时，需要给一些环境变量赋值：

设定程序所在的目录
$ export ICORN2_HOME=/opt/biosoft/ICORN2/
设定运行的线程数
$ export ICORN2_THREADS=24
设定输出信息的多或少,对 debug 有用
$ export ICORN2_VERBOSE=2

4.2 运行软件 $ /opt/biosoft/ICORN2/icorn2.sh readroot 350 genome.fasta 1 3

以上命令第 2 个参数是数据的 insert size； 1 和 3 代表迭代的起始和终止，表示迭代 3 次。如果需要继续迭代，则设置起始为 4 。上述程序的结果文件为 ICORN2.Query.contigs.fa.4 。在 ICORN2_3 文件夹中也含有该 fasta 文件，同时程序也生成了一个 gff 文件。作者推荐将这两个文件载入到 artemis 基因组浏览器中进行查看。

5. 思考

有些人在有参考基因组的状况下做了重测序，总是想要得到其重测序的基因组结果。可以使用 ICORN 对参考基因组进行修正，即得到了对应其品种的基因组序列。

1. 简介

Edena 是一个基于 overlaps–graph-based 的 de novo assembler。仅能很好地使用 Illumina 数据进行基因组组装。软件能使用 direct-reverse (paired-ends) 和 reverse-direct (mate-pairs) 的 Illumina 数据。同时要求 reads 的长度必须一致，否则所有 reads 会从 3′ 端开始截短到最短的 reads 的长度。

2. Edena 下载和安装

$ wget http://www.genomic.ch/edena/EdenaV3.131028.tar.gz
$ tar zxf EdenaV3.131028.tar.gz -C　/opt/biosoft
$ cd /opt/biosoft/EdenaV3.131028/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/EdenaV3.131028/bin/' >> ~/.bashrc
$ source ~/.bashrc

3. Edena 的使用

软件直接使用 Edena 命令即可。包含两个步骤： overlapping 和 assembling。

3.1 overlapping

使用例子：

$ edena -nThreads 24 -DRpairs frag1_1.fq frag1_2.fq frag2_1.fq frag2_2.fq\
 -RDpairs jump1_1.fq jump1_2.fq jump2_1.fq jump2_2.fq

最后的结果文件为 out.ovl

3.2 Assembling

使用例子：

$ edena -e out.ovl

最终的结果文件为 out_contigs.fasta 。

1. CISA 简介

CISA (Contig Integrator for Sequence Assembly). 主要用于细菌基因组（小基因组）的整合。该软件使用简单，算法较好。
CISA 软件官网：http://sb.nhri.org.tw/CISA/en/CISA
CISA参考文献：Lin S H, Liao Y C. CISA: contig integrator for sequence assembly of bacterial genomes[J]. PloS one, 2013, 8(3): e60843.

2. CISA 算法

CISA分4步进行，如下图所示：

2.1 鉴定代表性的 contigs 并对其进行延伸

首先，将所有 Assemblies 的 contigs 进行合并，并按长度进行排序。然后调用 NUCmer 软件按从长到短的顺序对序列逐个进行基因组比对。通过此方法来得到序列在各个 Assemblies 之间的重叠状况，并将结果写入到 CISA1/explained.txt(序列重叠信息) 和 CISA1/Extend_info(序列延伸信息) 文件中。
设定的阈值为：非代表性序列与代表性序列的比对结果达到 >95% 覆盖度 和 >95% 的一致性； 或非代表性序列超出了代表性序列的两端时，比对结果 >80% 覆盖度 和 >95% 的一致性。
通过这种方法，最终得到延长的代表性序列，输出到 R1_Contigs.fa 文件中。
如图所示：对 3 个 Assemblies 进行 CISA 分析。首先找到所有 contigs 序列中最长的 contig 序列（蓝色线条表示）进行全基因组比对，得到了第 1 轮结果；然后使用剩下的序列（+号后面的序列）挑出最长的 contig 序列（蓝色线条表示），对剩下的序列使用 NUCmer 进行比对，得到第 2 轮结果； 继续使用剩下的序列挑出最长的 contig 序列（黑色线条表示），最为代表性序列进行比对；最终得到了 4 条（蓝蓝黑红）代表性序列，其中的空心椭圆和实心箭头分别代表延伸的位置和方向。

2.2 对 misassembled 的 contigs 进行去除或截断, 并对 contigs两端不确定区域进行截短

将 R1_Contigs.fa 的序列使用 NUCmer 比对到所有的 contigs 上。通过比对结果，从而对 misassembled 的 contigs 进行去除或截断, 或对 contigs两端不确定区域进行截短。
如图所示：
Phase(2)左图表示 1 条代表性序列比对到了 2 个不同的 contigs 的中间，则表明该代表性序列是错误的组装，则去除该代表性序列，同时引进相应的能比对到该代表性序列的 contigs。这些信息在 CISA2/Remove_Info 文件中。
Phase(2)右图表示 1 条代表性序列中间区域没有对应的 contigs 能匹配，则表明该序列是错误的组装，则将该区域进行去除，截断该代表性序列。有关Gaps 的信息在 CISA2/Gaps 文件中，若 Gap 长度大于其代表性序列长度的 95%，则去除该序列。
此外，对代表性序列首尾不确定性区域的信息在 CISA2/clip_info 和 CISA2/clip_out 文件中。
最终，对 misassembled 的 contigs 进行处理后，得到结果文件 R2_Contigs.fa。该文件的序列总长度偏大。

2.3 通过 Overlap 信息将 contigs 融合

如果 2 个 contigs 有 end-to-end 重叠，并且重叠区对较小 contigs 长度的比例 >30%, 则将其进行融合。使用多轮 blastn 进行迭代来融合 contigs 序列。同时也进行重复序列鉴定。
融合信息写入到每轮文件夹下的 Extend_info 文件中，同时每轮生成 temp.fa 文件，作为下一轮 blastn 的数据库构建的输入文件。
重复序列信息文件写入到每轮文件夹下的 Repeat_Region.txt 文件中。

2.4 通过 small overlap 进行 contigs 融合

上一步骤将 contigs 进行融合需要 >30% 覆盖度的 overlap，同时也鉴定了重复序列。如图所示：Phase3中，黑色与红色线条重叠区比例大，被整合到一起，该序列和其他两条蓝色序列的overlap比较small，不能进行融合。但是鉴定了其重复序列（绿色块表示）后，通过本步骤，若重复区域长度小于其 small overlap，则能将 contigs 进行融合。
此步骤和上一步骤类似，使用 blastn 进行多轮运算。同时将结果写入到每轮文件夹的 Extend_info 和 temp.fa 文件中。
最终的输出文件为 CISA.fa (由 cisa.config 配置文件指定的名称)

3. CISA 的下载与安装

$ wget http://sb.nhri.org.tw/CISA/upload/en/2014/3/CISA_20140304-05194132.tar
$ tar xf CISA_20140304-05194132.tar -C /opt/biosoft

4. CISA 的使用

4.1 创建 merge.config 配置文件

配置文件如下：

count=3 
data=abyss.fasta,title=abyss  
data=velvet.fasta,title=velvet
data=allpathslg.fasta,title=allpathslg
Master_file=merge_contigs.fa
min_length=100 (default:100)
Gap=11

4.2 Merge.py

$ /opt/biosoft/CISA1.3/Merge.py merge.config

根据上面 merge.config 配置文件信息，对 3 个 Assemblies 的 contigs 序列进行提取，并整合到 merge_contigs.fasta 文件中，同时输出各个 Assemblies 的 contigs size。

4.3 创建 cisa.config 配置文件

配置文件内容如下：

genome=47814562 (填写Assemblies中最大的 contigs size)
infile=merge_contigs.fa
outfile=cisa.contig.fa
nucmer=/opt/biosoft/MUMmer3.23/nucmer
R2_Gap=0.95 (default:0.95)
CISA=/opt/biosoft/CISA1.3
makeblastdb=/opt/biosoft/ncbi-blast/bin/makeblastdb
blastn=/opt/biosoft/ncbi-blast/bin/blastn

4.4 CISA.py

$ /opt/biosoft/CISA1.3/CISA.py cisa.config

启动主程序进行运算。基因组越大，Assemblies 数目越多，耗时越多（呈几何级数增加）。

5. 思考

CISA 应该要修改下，在第 2 步进行并行化运算，以提高速度。
可以分批进行两两融合，直到最终融合。多个两两融合进行并行化运算。
进行大基因组的融合时，CISA 不适用了。

进行多个 genome Assemblies 整合的软件有很多：GAA, GAM, GAM-NGS, GARM, CISA, minimus2, Reconciliator, MAIA 等。

我对其中的一些进行了尝试，结果下：

1. GAA：软件很快就给出了结果，但是结果明显是错的。我不知道是不是用错了，但是软件运行过程中没有报错。

2. GAM-NGS： 软件安装可能会出现一些问题。有时候能出结果， 将 master 和 slave 反过来运行却会有报错。同时使用 bowtie2 的结果作为输入会报错。该程序安装和运行太过容易出错。 没有使用 MUMer ，虽然说是能节约时间，其实生成 sorted 的 bam 文件也极为耗时的。 最后，对于我的数据，该软件的整合效果太差。

3. GARM：软件的 pipeline 明显没写好，直接报错，连使用其 example data 都狂报错。 一看下载的程序是数年前的，估计作者也没修改其程序了。

4. Reconciliator： 从其发表文章所述的链接，完全下载不到该软件。

5. MAIA： 貌似需要使用 matlab，没装 matlab 切需要用钱购买的软件，这个就麻烦了。

6. minimus2： 使用貌似需要一些其他格式的文件，貌似比较麻烦，同时基因组整合的效果听说比较差。我个人就不去费精力验证了。

7. GAM: 貌似 GAM-NGS 是在该软件基础上做的，那么这个估计也好不到哪儿去。

8. CISA： 该软件结果还不错，使用起来也非常简单。唯一的缺点是运算极度耗时，特别是运行第2步的时候，对于较大的基因组，使用此软件几乎是不可能的事情。推荐使用此软件。

1. L_RNA_scaffolder简介

L_RNA_scaffolder能利用长的转录子序列来辅助基因组的scaffolding拼接。这些长的转录子序列可以是454/sanger/Ion_Torrent的测序结果，也可以是illumina测序的de novo组装结果。
在使用此软件进行scaffolding之前，需要将这些长的转录子序列使用BLAT软件map到基因组上，得到PSL（without heading）文件。将此PSL文件和基因组的初步组装结果文件输入到L_RNA_scaffolder,很快能得到结果文件。

2. L_RNA_scaffolder的安装

$ wget http://www.fishbrowser.org/software/L_RNA_scaffolder/downloads/L_RNA_scaffolder.tar.gz
$ tar zxf L_RNA_scaffolder.tar.gz
$ cd L_RNA_scaffolder
运行软件的主程序，不加参数给出帮助信息
$ ./L_RNA_scaffolder.sh

3. L_RNA_scaffolder的使用

先使用blat，将转录子序列比对到基因组
$ blat genome.fasta transcripts.fasta -noHead out.psl
再使用L_RNA_scaffolder.sh得到最终的基因组文件
$ mkdir L_RNA_scaffolder_output
$ /opt/biosoft/L_RNA_scaffolder/L_RNA_scaffolder.sh -d /opt/biosoft/L_RNA_scaffolder/ -i out.psl -j genome.fasta -o ./L_RNA_scaffolder_output
最终的结果文件为./L_RNA_scaffolder_output/L_RNA_scaffolder.fasta

1. musket 简介

Musket is a an efficient multistage k-mer spectrum based corrector for illumina short read data。
其发表文献为：Liu Y, Schröder J, Schmidt B. Musket: a multistage k-mer spectrum-based error corrector for Illumina sequence data[J]. Bioinformatics, 2013, 29(3): 308-315.
该软件使用简单，能进行多线程运行，速度快。

2. musket 下载与安装

$ wget http://cznic.dl.sourceforge.net/project/musket/musket-1.1.tar.bz
$ tar zxf musket-1.1.tar.bz -C /opt/biosoft
$ cd /opt/biosoft/musket-1.1/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/musket-1.1/' >> ~/.bashrc
$ source ~/.bashrc

在软件编译前，可以修改 Makefile 中第二行内容 “MACROS = -DMAX_SEQ_LENGTH=200 -DMAX_KMER_SIZE=28”。此行表示，软件允许的最大 read 长度为 200， 最大的 kmer 为 28。可以根据自己的需求修改此值。

3. musket 参数

-k int unit
此参数后有 2 个值，前者为 k-mer 的大小，默认值为 21。后者是 distinct k-mers 的数目，如果是全基因组测序的数据，则可以理解为基因组的大小。当然，由于 read 中可能存在错误，因此 distinct k-mers 的数目比基因组要大。如果内存够用，则我个人推荐设置为基因组大小的 2 倍。该值的设置主要用于平衡内存的消耗，对结果无影响。
-o str
此参数设置出处文件的文件名。设置此参数后，则不管输入文件有多少，只会输出一个文件。此参数和下一个参数 "-omulti" 是互斥的，仅能设置其中一个。
-omulti str
此参数设置出处文件的文件名。设置此参数后，每个输入文件都会对应一个输出文件。
-p int
程序能使用的线程数。默认值为 2。该值要设置为 >= 2。
-zlib int
设置输出的文件是否压缩，默认值为 0， 即不压缩。
-maxtrim int
允许 trim 掉的最大碱基数。该参数默认值为 0， 即表示不进行 trim。若此参数不为 0， 则程序在对一个 read 进行 error correction 后，再对该 read 进行 trim。trim 的方法为：程序寻找此 read 中没有错误碱基的子序列，然后取出最长的子序列，若最长子序列长度满足此参数的设定，则进行 trim，若不满足参数设定，则不进行 trim。
-inorder
此参数是输出文件的 reads 顺序和输入文件一致。在进行 paired reads 的修正时，一定要加上该参数。
-lowercase
是输出文件中的被修正的碱基为小写字母。默认值为 0 。
-multik bool
设置是否可以使用多个 k-mer sizes 来进行 error correction。默认值为 0。当设置为 1 时， 则可以设置多个 k-mer 值，如 "-k 21 536870912 -k 23 268435456"。

4. 常用的一些例子

$ musket -k 21 536870912 -p 12 reads.fa reads2.fa -o output.fa
$ musket -p 12 reads.fa reads2.fa -o output.fa
$ musket reads.fa reads2.fa -o output.fa -inorder

$ musket lib1_reads1.fastq lib1_reads2.fastq lib2_reads1.fastq lib2_reads2.fastq -omulti output -inorder
上述命令生成 4 个文件 output.0 output.1 output.2 output.3 按顺序分别对应前 4 个输入文件。

SSPACE 能利用 paired reads 的比对结果，将 contigs 或 scaffolds 连接成 scaffolds。其参考文献：Boetzer M, Henkel C V, Jansen H J, et al. Scaffolding pre-assembled contigs using SSPACE[J]. Bioinformatics, 2011, 27(4): 578-579.

1. 安装 SSPACE

软件下载页面：http://www.baseclear.com/lab-products/bioinformatics-tools/sspace-standard/。

$ tar zxf SSPACE-STANDARD-3.0_linux-x86_64.tar.gz
$ ./SSPACE-STANDARD-3.0_linux-x86_64/SSPACE_Standard_v3.0.pl

解压缩软件包后，运行软件文件夹中的 perl 程序即可运行 SSPACE。软件主目录下包含一些软件使用说明和示例等,其中 README 文件描述得非常详细。

2. SSPACE 使用方法

2.1 library 文件

首先要建立一个描述 library 信息的文本文件，例如：

Lib1 bwa file1.1.fasta file1.2.fasta 400 0.25 FR
Lib1 bowtie file2.1.fasta file2.2.fasta 400 0.25 FR
Lib2 bwasw file3.1.fastq file3.2.fastq 4000 0.5 RF
Lib2 TAB file4.tab 4000 0.5 RF
Lib3 TAB file5.tab 10000 0.5 RF
unpaired bowtie unpaired_reads1.fasta
unpaired bwasw unpaired_longreads1.gz

此 library 文件由多列组成，列与列之间由 1 个 空格 或 tab 分隔，各列意义如下：

第 1 列： library 名称。程序运行过程中产生的临时文件以此来命名； 多个行可以拥有同一个 library 名称，则其具有相同的 library 设置和不同的数据文件； 同时，libraries 必须按 insert size 来排序，inert size 最小的必须放到第一行，这是因为进行 scaffold 构建时，按此文件提供的 libraries 的顺序来输入数据的； unpaired reads， 则第一列是 ‘unpaired’。
第 2 列： 将 reads 比对到基因组上所使用的软件名， 可以为 bowtie 、 bwa 和 bwasw 等； 如果输入的数据是 reads 比对过后的 tab 格式结果，则此列为 “TAB”。
第 3，4 列： Fasta 或 Fastq 格式的双末端测序文件，并且文件中成对的 paired reads 必须在两个文件中并处于相同的行号上，同时，软件读取数据与序列的 headers 无关。如果是 unpaired reads，则仅需要第 3 列，为 tab 格式的 reads mapping 结果，过后详述。
第 5,6 列：第 5 列为 insert size 的期望值； 第 6 列为 insert size 允许的最小偏差。 比如，这两列值分别为 4000 和 0.5，则 insert size 在 2000-6000 之间的 pairs 才是有效 pairs。
第 7 列：paired-reads 的方向，有 FF，FR，RF 或 RR 几种选项。

2.2 程序参数

-l 输入的 library 文件
-s 输入的 Fasta 文件
-x 是否对 contigs 进行延长。其值可以为 0 或 1。 1 表示进行延伸，0 表示不延伸。默认值为 0。

延伸参数：
-m 进行延伸时，read 和基因组序列最小的 overlap。此值越大，则结果越准确，同时耗内存越少。推荐此值接近最长的 read 的长度。比如，对于 26 bp 长度的 reads， 该值适合设为 32～35。 默认此值为 32 。此值取值范围为 15～50 。软件运行时，将 unmapped reads 全部打断成 m+1 长度的序列，这些序列用于进行 contigs 的延伸。
-o 进行延伸时，延伸 1 个碱基需要的最小 reads 数。此值越大，则结果越准确。默认值为 20 。
-r 进行延伸时，延伸 1 个碱基，此碱基在所有匹配的 reads 中的最小比例。此值越大，则结果越准确。默认值为 0.9 。

Scaffolding 参数：
-k 将两个 contigs 连接成 scaffold 时，需要的最小的 reads pairs 数目。默认值为 5 。
-a 将两个 contigs 连接成 scaffold 时，这两个 contigs 之间的连接数 与 其和其它 contigs 的连接数之间的最小比值。此值越大，则结果越准确。默认值为 0.70
-n 在 scaffold 中，将两个邻近的 contigs 合并到一起需要的最小的 overlap。默认值为 15。
-z 进行 scaffolding 时，允许的最小的 contig 长度。低于此长度的 contig 将不能用于进行 scaffold 组装。默认值为 0 。较长的 contigs 产生的 scaffolds 比较可信； 而小于 100bp 的 contigs 容易是重复序列。

bowtie 比对参数：
-g 使用 bowtie 进行比对时，允许的最大 gaps 数。默认值为 0

其它参数：
-T 设定运行的线程数。默认值为 1。
-b 输出文件夹名及文件夹内的文件前缀。
-S 当程序正在运行时，跳过读取 reads 的阶段。和 -b 参数结合使用，则可以同时运行多个 SSPACE 程序，对每个程序设置不同的参数，这样能较快得到较好的结果。
-v verbose mode
-p 生成可供可视化的 .dot 文件。

2.3 其它工具

SSPACE 提供了一些其它比较有用的小工具：

estimate_insert_size.pl 用于计算 insert size。此程序计算的结果有些问题。
fastq_qualitytrim_pairs.pl 对 reads pairs 进行质量控制的程序。

sam_bam2tab.pl 将 bam sam 文件转换为 tab 格式的程序。