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

PostGreSQL 数据库的简单使用

1. PostGreSQL 简介

PostgreSQL 是一种非常复杂的对象-关系型数据库管理系统(ORDBMS),也是目前功能最强大,特性最丰富和最复杂的自由软件数据库系统。有些特性甚至连商业数据库都不具备。这个起源于伯克利(BSD)的数据库研究计划目前已经衍生成一项国际开发项目,并且有非常广泛的用户。
有人将 PostgreSQL 和 MySQL 作了比较, PostgreSQL 貌似完胜。

2. PostGreSQL 的安装

# yum install postgres*
# /etc/init.d/postgresql start       启动数据库,第一次启动则会初始化数据库。

3. PostGreSQL 的简单使用

参考自: PostgreSQL 8.1 中文文档

3.1 登录 PostGreSQL 的命令

使用 psql 命令登录 PostGreSQL 数据库。
初始化 PostGreSQL 数据库后,则会有名为 postgres 和 tempalte1 的数据库,有一个名为 postgres 的用户,无密码登录。
例如,使用 postgres 用户登录 postgres 数据库

$ psql postgres postgres
psql (8.4.18)
Type "help" for help.

postgres=#

psql的用法与常用参数:

psql [options] [DBNAME [USERNAME]]

DBNAME
    设置要登录的数据库,默认为当前系统用户名
USERNAME
    设置登录数据库的用户,默认为当前系统用户名
-h | --host=HOSTNAME
    设置连接的服务器名,default: localhost
-p | --port=PORT
    设置端口, default: 5432
-U | --username=USERNAME
    设置登录的用户名,默认为当前系统用户名

3.2 简单常用的 postgresql 命令行

登录 postgresql 数据库后,提示符以登录的数据名加等于号表示。此时,可以输入 postgresql 命令行用于辅助数据库操作。这些命令行以 \ 开头。几个简单常用的命令如下:

\c[onnect] DBNAME
    切换数据库
\l[+]
    列举数据库
\dt
    列举表
\d TABLENAME
    查看表的结构
\q
    退出 postgresql
\cd [DIR]
    切换 linux 系统的工作目录
\! [COMMAD]
    执行 linux 系统命令。比如 \! ls --color 用于打印出当前目录的文件名; 但是 \! cd 命令无效。
\?
    展示所 postgresql 的命令行帮助文档
\h CREATE ROLE
    \h 为查看 SQL 命令的帮助信息,以上命令为查看 CREATE ROLE 命令的详细帮助。

3.3 简单常用的 SQL 命令

使用 SQL 语言对 postgresql 的数据库和表进行操作,初学者需要的 SQL 命令:

CREATE DATABASE chado
    创建名为 chado 的数据库
DROP DATABASE chado
    删除名为 chado 的数据库
CREATE ROLE chenlianfu [PASSWORD '123456']
    创建名为 chenlianfu 的用户
DROP ROLE chenlianfu
    删除名为 chenlianfu 的用户

3.4 PostGres 的 linux 命令

在 linux 终端下,可以使用一些命令来操作 PostGres 数据库。

createdb
    创建数据库
createlang
    给数据库安装语言
createuser
    用户创建用户
dropdb
    删除数据库
droplang
    从数据库删除语言
dropuser
    删除用户
initdb
    初始化数据库
pg_dump
    将数据库转换为文本文件,用于备份

基因预测流程

1. 获取准确的转录子序列

利用 RNA-seq 数据进行转录组 de novo 组装。 推荐使用 trinity 软件进行有基因指导的 de novo 组装,只进行 inchworm 组装,这样能获得比较准确的转录组序列。

2. 获取完整的基因模型

使用 PASA 将上一步骤得到的转录子序列比对到基因组上。根据比对结果,提取完整并且长度较长的基因模型。

3. HMM training

根据上一步骤的基因模型进行 AUGUSTUS 和 SNAP 的 training,得到这 2 个基因预测软件的 HMM 参数文件。

4. 由 RNA-seq 数据得到 hints

对基因组进行重复序列分析,对 DNA 转座子和逆转录转座子进行 hardmask (即对该区域使用 N 屏蔽),对 low-complexity 区域进行 softmask(即将该区域碱基换成小写)。 将 RNA-seq 数据比对到屏蔽了重复序列的基因组,得到 hints 信息。

5. 获取低表达区域来自其它物种的保守蛋白

在屏蔽重复序列的基础上,对上一步骤得到 hints 区域进行hardmask。 从 NCBI 的 Taxonomy 数据库下载蛋白质序列,并将这些序列比对到屏蔽了重复序列和基因表达区域的基因组上,设定阈值筛选出保守蛋白。这些蛋白比对到了表达量低的基因区,有利于这些基因的预测。

6. 使用 MAKER 进行第一遍基因预测

使用 MAKER 进行基因预测,输入文件有:

基因组序列: 基因组装的结果文件
重复序列文库: 使用 ReapeterModeler 得到
HMM 文件: AUGUSTUS, SNPA 和 GeneMark-ES 的 HMM 文件
转录子序列: PASA 的结果文件
蛋白质序列: 筛选出来的蛋白序列

7. 使用 MAKER 进行第二遍基因预测

根据上一步的预测结果,提取 AED 值较低的基因模型,进行第二遍 HMM training,然后再次使用 MAKER 进行基因预测。

8. 进行基因模型的人工校正

samtools 和 bamtools 安装

1. samtools 安装

SAMtools 新的 1.0 版本改进比较大。 SAMtools 包含 3 部分:

Samtools : 对 SAM/BAM/CRAM 格式文件进行读、写、索引和查看操作。
BCFtools : 对 BCF2/VCF/gVCF 格式文件进行读写操作; 对 SNP 和 short indel 进行 calling/filtering/summarising
HTSlib : 用于高通量数据读写操作的 C library

新版的 samtools 将 Samtools 和 BCFtool 分离开,成为 2 个独立的软件包。 这 2 个都运用了 HTSlib,并且其软件包中都包含了 HTSlib 。
新版不的 Samtools 的官方网址更改为了: http://www.htslib.org/

1.1 安装 Samtools

$ wget http://downloads.sourceforge.net/project/samtools/samtools/1.1/samtools-1.1.tar.bz2
$ tar jxf samtools-1.1.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/samtools-1.1/
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/samtools-1.1/' >> ~/.bashrc

1.2 安装 BCFtools

$ wget http://sourceforge.net/projects/samtools/files/samtools/1.1/bcftools-1.1.tar.bz2
$ tar jxf bcftools-1.1.tar.bz2 -C /opt/biosoft/
$ cd /opt/biosoft/bcftools-1.1
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/bcftools-1.1' >> ~/.bashrc

2. Bamtools

Bamtools 是 Samtools 的 c++ 版本的 Samtools, 其官网: https://github.com/pezmaster31/bamtools

$ wget https://github.com/pezmaster31/bamtools/archive/master.zip -O bamtools-master.zip
$ unzip bamtools-master.zip -d /opt/biosoft/
$ cd /opt/biosoft/bamtools-master/
$ mkdir build
$ cd build/
$ cmake ..
$ make -j 8
$ echo 'PATH=$PATH:/opt/biosoft/bamtools-master/bin/' >> ~/.bashrc 
$ source ~/.bashrc