tRNAscan-SE 的使用

1. tRNAscan-SE 简介

tRNAscan-SE 能在基因组水平上进行 tRNA 扫描。该软件实际上是一个 perl 脚本,整合了 tRNAscan、 EufindRNA 和 Cove 这 3 个独立的 tRNA 检测软件。tRNAscan-SE 首先调用 tRNAscan 和 EufindRNA 鉴定基因组序列中的 tRNA 区域,然后调用 Cove 进行验证。这样既保证了前者的 sensitivities, 又保证了后者较低的假阳性概率,同时在搜索速度上提升了很多。
有关 tRNAscan-SE 的详细说明,参考其本地化软件包中的 man 文档。
tRNAscan-SE 的网页版:http://lowelab.ucsc.edu/tRNAscan-SE/。但一次最多只能进行 5M bp 序列的 tRNA 预测。

2. tRNAscan-SE 本地安装

$ wget http://lowelab.ucsc.edu/software/tRNAscan-SE.tar.gz
$ tar zxf tRNAscan-SE.tar.gz
$ cd tRNAscan-SE-1.3.1
$ perl -p -i -e 's#\$\(HOME\)#/opt/biosoft/tRNAscan-SE-1.3.1#' Makefile
$ make && make install
$ make testrun
$ echo 'PATH=$PATH:/opt/biosoft/tRNAscan-SE-1.3.1/bin/' >> ~/.bashrc
$ echo 'PERL5LIB=$PERL5LIB:/opt/biosoft/tRNAscan-SE-1.3.1/bin/' >> ~/.bashrc
$ source ~/.bashrc

3. tRNAscan-SE 的使用

常用例子与主要参数:

$ tRNAscan-SE -o tRNA.out -f rRNA.ss -m tRNA.stats genome.fasta

-A 适合于古细菌。该参数选择了古细菌特异性的covariance model(cm),同时稍微放宽了 EufindtRNA 的 cutoffs。
-B 适合于细菌。默认情况下,不选择,-A -B -G 或 -O 参数,则适合于真核生物。
-G 适合于古细菌,细菌和真核生物的混合序列。该参数使用 general tRNA covariance model。
-O 适合于线粒体和叶绿体。选择该参数,则仅使用 Cove 进行分析,搜索速度会很慢,同时也不能给出 pseudogenes 检测。
-i 使用 Infernal cm analysis only。该参数设置后,需要 cmsearch 命令,但是 tRNAscan-SE 软件包中貌似没有该程序,最终无法运行。
-C 仅使用 Cove 进行 tRNA 分析。虽然从一定程度上提高了准确性,但是会极慢,当然不建议了。
-o <file> 将结果保存到文件。
-f <file> 将 tRNA 的二级结构结果保存到文件
-m <file> 将统计结果保存到文件。

4. tRNAscan-SE 的结果说明

在真核生物中,tRNA 由 RNA 聚合酶III 在核内转录生成 pre-tRNA, 再进行加工生成有功能的 tRNA 分子(特别是一些 tRNA 序列还含有内含子)。若 tRNA 存在内含子,则结果文件中第 7 8 列会给出内含子区间,否则其值为 0 。
tRNAscan-SE 的结果中, 如果 begin 比 end 的值大,则表示 tRNA 在负义链上。
有些结果中第 5 列为 pseudogene, 这表示其一级或二级结构比较差。
最后一列是 Cove Score,该分值最低阈值为 20 。该值是一个 log ratio值。ratio 是符合 tRNA covariance model概率与随机序列模型概率的比值。
当然,最后最好是将表格格式结果转换为 GFF3 结果,以利于在基因组上的可视化。

上传基因组数据到NCBI

创建 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 的制作依然需要自己手工制作,再上传。

GENEWISE 的使用

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程序能进行该项处理(我没有用过该程序,我自己写了想要的程序)。

使用 ICORN 进行基因组核苷酸的修正

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 为止。 http://icorn.sourceforge.net/workflow.jpg

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 对参考基因组进行修正,即得到了对应其品种的基因组序列。

使用 Edena 进行基因组组装

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 。

使用 CISA 进行多个 genome Assemblies 的整合

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步进行,如下图所示:
CISA算法图

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 的整合的经验

进行多个 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步的时候,对于较大的基因组,使用此软件几乎是不可能的事情。推荐使用此软件。

使用 L_RNA_scaffolder 进行转录组序列对基因组序列的 rescaffolding

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

使用 musket 进行 reads correction

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 个输入文件。