使用 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

MAKER的使用

1. MAKER 简介

MAKER 是进行基因组结构注释的一整套流程。它综合了ab initio方法和 evidence 的方法进行基因组注释。使用起来非常方便快捷,一个命令能搞定全套基因组结构注释。

2. MAKER下载与安装

http://yandell.topaz.genetics.utah.edu/cgi-bin/maker_license.cgi页面提交信息并下载最新版本MAKER。

此外,安装 maker 需要先安装一些软件

必须安装:
Perl 5.8.0 or higher
BioPerl 1.6 or higher
WU-BLAST 2.0 or higher or NCBI-BLAST 2.2.X or higher 
RepeatMasker 3.1.6 or higher
    RepeatMasker requires a repeat library, available from Repbase. 
Exonerate 1.4 or higher

可选安装:
SNAP version 2009-02-03 or higher (for eukaryotic genomes).
Augustus 2.0 or higher (for eukaryotic genomes) 
GeneMark-ES (for eukaryotic genomes)
FGENESH 2.6 or higher - requires licence (for eukaryotic genomes)
GeneMarkS (for prokaryotic genomes) 
snoscan
MPICH2 (install MPICH2 with the --enable-sharedlibs flag)
安装perl模块
$ cpan -i BioPerl Bit::Vector DBD::SQLite DBI Error Error::Simple File::NFSLock File::Which forks forks::shared Inline Inline::C IO::All IO::Prompt PerlIO::gzip Perl::Unsafe::Signals Proc::ProcessTable Proc::Signal threads URI::Escape

安装mpich2
$ sudo yum install mpich2*
使用 MPI 运行 Maker 需要先运行服务 mpd。
正常运行 mpd 命令,需要正确的主机名。否则可能会遇到报错。需要修改 /etc/hosts 文件,将主机名的 IP 对应正确。

安装maker
$ wget http://yandell.topaz.genetics.utah.edu/maker_downloads/2D3B/22C6/BE99/41DBB7F9B1D181B2CDCE1E49DFE6/maker-2.31.8.tgz
$ tar zxf maker-2.31.8.tgz 
$ cd maker/src/
$ perl Build.PL 
$ ./Build install
$ echo 'PATH=$PATH:/opt/biosoft/maker/bin/' >> ~/.bashrc
$ source ~/.bashrc 

3. MAKER 的使用

3.1 准备 MAKER 的输入文件

MAKER 需要的输入文件如下:

1. genome.fasta
    基因组序列文件。
2. inchworm.fasta
    转录组组装结果文件。此文件内容代表转录子序列或 EST 序列。为了增加其准确性,我个人选用了 trinity 中 inchworm 的结果文件,而不是转录组最终的组装结果文件。
3. proteins.fasta
    临近物种的蛋白质序列。可以从 NCBI 的 Taxonomy 数据库中下载。
4. consensi.fa.classified
    适合本物种的重复序列数据库。该文件是一个fasta文件,其序列代表着本物种的高重复序列。该文件使用 RepeatModeler 制作。
5. Augustus hmm files [species]
    适合本物种的 Augustus 的 HMM 文件。该文件以文件夹形式存放在 /opt/biosoft/augustus-3.0.3/config/species/ 目录下。将 genome.fasta 和 inchworm.fasta 输入 PASA,得到用于 trainning 的基因,然后使用 Augustus 进行 training,得到 HMM 文件。
6. Snap hmm file [species.hmm]
    适合本物种的 SNAP 的 HMM 文件。该文件是一个单独的文件。将 genome.fasta 和 inchworm.fasta 输入 PASA,得到用于 trainning 的基因,然后使用 SNAP 进行 training,得到 HMM 文件。
7. Genemark_es hmm file [es.mod]
    适合本物种的 Genemark_es 的 HMM 文件。将 genome.fasta 输入 genemark_es 软件进行自我训练得到的 HMM 文件。
8. rRNA.fa
    rRNA 的序列文件。该文件使用 RNAmmer 对基因组进行分析得到。

得到上述输入文件后,使用 MAKER 命令得到初始化的 3 个配置文件。

$ maker -CTL

得到3个配置文件: maker_exe.ctl, maker_bopts.ctl 和 maker_opts.ctl。
maker_exe.ctl: 配置 Maker 需要调用的程序的路径
maker_bopts.ctl: 配置 blast 的各种阈值
maker_opts.ctl: MAKER 的主要配置文件,配置输入文件的路径,以及 MAKER 的使用流程。

3.2 使用 MAKER 进行基因组注释

先修改 maker_opts.ctl 配置文件,需要修改的内容如下:

genome=genome.fasta            # 基因组序列文件
est=inchworm.fasta             # EST 序列文件
protein=proteins.fasta         # 临近物种的蛋白质序列文件
model_org=fungi                # 选择repebase数据库的物种,对真菌则用fungi,也可以使用 all
rmlib=consensi.fa.classified   # 适合本物种的重复序列数据库
snaphmm=species.hmm            # snap traing 的 HMM 文件
gmhmm=es.mod                   # genemark_es 的 HMM 文件
augustus_species=species       # augustus 的 HMM 文件
est2genome=1                   # 使用 EST 序列进行基因预测
protein2genome=1               # 使用蛋白质序列进行基因预测
trna=1                         # 使用 tRNAscan 进行 tRNA 预测
snoscan_rrna=rRNA.fa           # 使用 Snoscan 利用 rRNA 序列进行 C/D box 类型的 snoRNAs 预测,不过我加入了这个信息后,maker运行不成功,不知道为什么。
correct_est_fusion=1           # 进行融合 EST 的校正

并行化运行 maker。由于基因组序列比较大,MAKER 的计算量大,使用并行化能极大节约时间。

$ mpd &
$ maker                  # 不并行化运行
$ mpiexec -n 20 maker    # 使用20个线程并行计算

使用 lftp 进行 FTP 数据下载和上传

1. lftp 简介

适合于 FTP 操作的命令和软件比较多。 linux 下操作命令有 ftp, lftp 和 sftp; 图形化界面非常好用的有 FileZilla。 为了便于在服务器上使用命令行来进行便捷的 FTP 操作,我个人喜欢使用 lftp 。

2. lftp 的使用

lftp 的简单例子:
比如,下载测序公司释放的数据:

从诺和致源公司 FTP 站点 ftpdata.novogene.cn 以用户名 Novo_s4VcJc 密码 dpdATZ77 端口 2300 下载测序释放的数据:
$ lftp -e "get -c data/species.reads1.fq.gz; exit" \
  -u Novo_s4VcJc,dpdATZ77 -p 2300 ftpdata.novogene.cn

比如,上传测序原始数据到 NCBI 的 SRA FTP 站点:

$ lftp -e "put raw_data.reads1.fastq; exit" \
  -u sra,VfOiVJn1 -p 21 ftp-private.ncbi.nih.gov

lftp 的参数:

-f file
    执行该文件中的 FTP 命令,执行完毕后退出。
-c cmd
    将 FTP 命令直接写入到该参数后,执行其中的命令后退出。

-u user[,passwd]
    通过该参数输入用户名,同时可以选择输入密码。
-p port
    设置端口。默认为 21。
-e cmd
    在 lftp 命令后给出 FTP 位点, -u 和 -p 参数登录 FTP 后,再执行该参数后的 FTP 命令。 

3. 常用的 FTP 命令

直接在终端中输入 lftp 命令后,进入了 lftp 的操作界面,从而可以使用 FTP 命令进行操作了。

1. open

用于连接 FTP 站点,用法为:

open [-e cmd] [-u user[,pass]] [-p port] host|url

2. pwd 和 lpwd

进入 FPT 站点后, pwd 命令用于显示服务器端的当前工作目录; lpwd 用于显示本地机器上的当前工作目录。

3. ls 和 !ls

进入 FPT 站点后, ls 命令用于列出服务器端的目录内容; !ls 用于列出本地机器上的目录内容。

4. cd 和 lcd

进入 FPT 站点后,cd 命令用于在服务器端切换路径; lcd 用于在本地机器上切换路径。

5. pget

用于从 FTP 上下载数据。该命令能使用多个连接下载数据,从而加大下载速度,但是会增大服务端和网络的负载。用法和参数为:

pget [OPTS] rfile [-o lfile]
-c
    支持续传
-n int
    设置最大的连接数

6. get 和 mget

用于从 FTP 上下载数据。 get 用于下载一个文件, mget 用于下载多个文件。 用法和参数如下:

get [-E] [-a] [-c] [-O base] rfile [-o lfile] ...
mget [-c] [-d] [-a] [-E] [-O base] files

-c
    支持断点续传。
-E
    成功下载后,删除服务器端的数据
-a
    使用 ascii 模式下载,默认为二进制模式下载
-d
    创建和下载文件名一致的文件夹,并将文件保存到文件夹中
-O string
    指定下载文件存放的路径
-o string
    将 get 命令下载数据的数据保存到此文件中

7. put 和 mput

用于将数据上传到服务器端。put用于上传一个文件,mput用于上传多个文件。用法和参数为:

put [-E] [-a] [-c] [-O base] lfile [-o rfile]
mput [-c] [-d] [-a] [-E] [-O base] files

参数和 get/mget 参数一致。

8. mirror
用于将目标文件夹全部下载到指定目录中。其用法和参数如下:

mirror [OPTS] [source [target]]

mirror 的常用参数:
-c
    支持断点续传
-R
    反向 mirror,即将本地文件上传到服务器端
--parallel=N
    同时并行下载 N 个文件
--use-pget=N
    对每个文件使用 pget 下载,并设置 pget 的连接数。

使用Aspera从NCBI或EBI高速下载数据

1. Aspera简介

Aspera提供了大文件高速传输方案,适合于大数据的传输。客服端的使用是免费的。

2. Aspera下载和安装

Aspera下载网页: http://downloads.asperasoft.com/connect2/

$ wget http://d3gcli72yxqn2z.cloudfront.net/connect/bin/aspera-connect-3.5.1.92523-linux-64.tar.gz
$ tar zxf aspera-connect-3.5.1.92523-linux-64.tar.gz
$ sh aspera-connect-3.5.1.92523-linux-64.sh
$ echo 'PATH=$PATH:~/.aspera/connect/bin/' >> ~/.bashrc
$ source ~/.bashrc
$ ascp --help

软件安装在 ~/.aspera/connect/ 目录下。

3. Aspera 的使用

例子,使用 Aspera 高速下载 NCBI或 EBI 上的数据:

$ ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  --host=ftp-private.ncbi.nlm.nih.gov --user=anonftp --mode=recv \
  /sra/sra-instant/reads/ByRun/sra/ERR/ERR105/ERR105009/ERR105009.sra ./
$ ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/ERR/ERR105/ERR105009/ERR105009.sra ./

$ ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  --host=fasp.sra.ebi.ac.uk --user=era-fasp --mode=recv \
  /vol1/fastq/ERR105/ERR105009/ERR105009_1.fastq.gz ./
$ ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
  era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/ERR105/ERR105009/ERR105009_1.fastq.gz ./

Aspera 的用法和简单参数:

Aspera的用法:
$ ascp [参数] 目标文件 目的地址

Aspera的常用参数:
-T
    不进行加密。若不添加此参数,可能会下载不了。
-i string
    输入私钥,安装 aspera 后有在目录 ~/.aspera/connect/etc/ 下有几个私钥,使用 linux 服务器的时候一般使用 asperaweb_id_dsa.openssh 文件作为私钥。
--host=string
    ftp的host名,NCBI的为ftp-private.ncbi.nlm.nih.gov;EBI的为fasp.sra.ebi.ac.uk。
--user=string
    用户名,NCBI的为anonftp,EBI的为era-fasp。
--mode=string
    选择模式,上传为 send,下载为 recv。
-l string
    设置最大传输速度,比如设置为 200M 则表示最大传输速度为 200m/s。若不设置该参数,则一般可达到10m/s的速度,而设置了,传输速度可以更高。

JBrowse 的下载和安装

1. JBrowse 简介

2. JBrowse 下载和安装

下载 JBrowse 安装包,并安装:

$ wget http://jbrowse.org/releases/JBrowse-1.11.5.zip
$ wget http://jbrowse.org/releases/JBrowse-1.11.5-dev.zip
$ unzip JBrowse-1.11.5.zip -d /opt/biosoft
$ unzip JBrowse-1.11.5-dev.zip -d /opt/biosoft
$ cd /opt/biosoft/JBrowse-1.11.5
$ sudo ./setup.sh
$ echo 'PATH=$PATH:/opt/biosoft/JBrowse-1.11.5' >> ~/.bashrc

配置 apache:

# echo 'Alias /JBrowse "/opt/biosoft/JBrowse-1.11.5/"
<Directory "/opt/biosoft/JBrowse-1.11.5/">
    Options MultiViews ExecCGI Indexes FollowSymlinks
    AllowOverride None
    Order allow,deny
    Allow from all
</Directory>' > /etc/httpd/conf.d/JBrowse.conf
# /etc/init.d/httpd restart

打开软件自带的示例 Volvox JBrowse

3. JBrowse 配置