同一物种得到了多个基因组组装结果，如何将组装结果进行合并，从而使基因组更完整？
Parrish N, Sudakov B, Eskin E. Genome reassembly with high-throughput sequencing data[J]. BMC genomics, 2013, 14(Suppl 1): S8.一文中提到：
A number of software packages have been developed in recent years with the aim of utilizing a set of reference genomes to produce a more optimized scaffolding, or layout, of the contigs produced in de novo assembly. OSLay [20] uses a maximum-weight matching algorithm to identify likely neighboring contigs. Treecat [21] builds a fully connected graph of the contigs, with edges weighted by the distance between syntenic regions in the reference, and attempts to find a minimum-weight Hamiltonian path through the graph using a greedy heuristic. Finally, PGA [22] uses a genetic algorithm to search the space of possible contig orderings. By relying on the contigs produced through de novo assembly, however, these methods may not take full advantage of the reference genome.
从该篇文章开始，寻找方法：

1. AMOScmp

引用文献：Pop M, Phillippy A, Delcher A L, et al. Comparative genome assembly[J]. Briefings in bioinformatics, 2004, 5(3): 237-248.

AMOScmp applies a modified MUMmer algorithm to a newly sequenced genome by mapping it onto a reference genome.

Hawkeye and AMOS are available open source at http://amos.sourceforge.net.

2. Projector 2

引用文献：van Hijum S A F T, Zomer A L, Kuipers O P, et al. Projector 2: contig mapping for efficient gap-closure of prokaryotic genome sequence assemblies[J]. Nucleic acids research, 2005, 33(suppl 2): W560-W566.

文章中写道：
Projector 2 has several distinctive features: a user-friendly web interface, automatic removal of repetitive elements (repeat-masking) and automated primer design for gap-closure purposes. The web interface is freely accessible at http://molgen.biol.rug.nl/websoftware/projector2.

3. OSLay

引用文献：Richter D C, Schuster S C, Huson D H. OSLay: optimal syntenic layout of unfinished assemblies[J]. Bioinformatics, 2007, 23(13): 1573-1579.

该文章中写道：
The underlying algorithm is based on maximum weight matching. The tool provides an interactive visualization of the computed layout and the result can be imported into the assembly editing tool Consed to support the design of primer pairs for gap closure.

OSLay is freely available from: http://www-ab.informatik.unituebingen.de/software/oslay

4. PGA

引用文献：Zhao F, Zhao F, Li T, et al. A new pheromone trail-based genetic algorithm for comparative genome assembly[J]. Nucleic acids research, 2008, 36(10): 3455-3462.

该文章中写道：
A pheromone trail-based genetic algorithm (PGA) was used to search globally for the optimal placement for each contig.An extended version of PGA can predict additional candidate connections for each contig and can thus increase the likelihood of identifying the correct arrangement of each contig. The software and test data sets can be accessed at http://sourceforge.net/projects/pga4genomics/.

5.Treecat

引用文献：Husemann P, Stoye J. Phylogenetic comparative assembly[J]. Algorithms for Molecular Biology, 2010, 5(3).

该文章中写道：
Our new algorithm for contig ordering uses sequence similarity as well as phylogenetic information to estimate adjacencies of contigs. An evaluation of our implementation shows that it performs better than recent approaches while being much faster at the same time.
The software is open source (GPL) and available within the Comparative Genomics – Contig Arrangement Toolsuite (cg-cat, http://bibiserv.techfak.uni-bielefeld.de/cg-cat webcite) on the Bielefeld Bioinformatics Server (BiBiServ).

1. PHYLIP简介

PHYLIP,the Phylogeny Inference Package, is a package of programs for inferring phylogenies (evolutionary trees). It can infer phylogenies by parsimony, compatibility, distance matrix methods, and likelihood. It can also compute consensus trees, compute distances between trees, draw trees, resample data sets by bootstrapping or jackknifing, edit trees, and compute distance matrices. It can handle data that are nucleotide sequences, protein sequences, gene frequencies, restriction sites, restriction fragments, distances, discrete characters, and continuous characters.

2. PHYLIP的下载和安装

下载PHYLIP，然后解压缩，得到phylip的安装文件包，其中包含doc、exe和src三个文件夹以及一个phylip.html文件。html文件结合doc中的文件在浏览器中打开是phylip详细的说明文件。当然，官网也提供了phylip 3.6的PDF版的Manual文件。

$ wget http://evolution.gs.washington.edu/phylip/download/phylip-3.695.tar.gz
$ tar zxf phylip-3.695.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/phylip-3.695/src/
$ cp Makefile.unx Makefile
$ make install

1. 重复序列简介和相关软件

参考自文献：Review：Identifying repeats and transposable elements in sequenced genomes: how to find your way through the dense forest of programs。E Lerat。Heredity (2010) 104, 520–533; doi:10.1038/hdy.2009.165; published online 25 November 2009

1.1 Repeats的分类

基因组中的repeats依据其序列特征分成2类：串联重复(tandem repeats) 和 散在分布在基因组中的重复序列(interspersed repeats).其中第二类主要是transposable elements(TEs).

第一类串联重复包含：microsatellites 或 simple sequence repeats（1-6个碱基为一个重复单元） 和 minisatellites（10-60个碱基的长序列为一个重复单元）.

TEs包含2种类型：class-I TEs通过RNA介导的(copy and paste）机制进行转座；class-II TEs通过DNA介导的(cut and paste）机制来转座. 前者称为retroelements，后者称为DNA transposons。

class-I TEs中主要由LTR(long terminal repeat)构成。LTR的部分序列可能具有编码功能。而non-LTR则包含2个子类：LINEs(long interspersed nuclear elements)和SINEs(short interspersed elements），其中前者可能具有编码功能，后者则没有。

class-II TEs中加入了一个子类 MITEs(miniature inverted repeat transposable elements),基于DNA的转座因子，但是确通过”copy and paste”的机制来转座(Wicker et al., 2007)。

1.2 鉴定重复序列的软件

对于不同的重复序列，需要使用不同的软件来进行鉴定。而鉴定的方法可以分为：基于library，基于重复序列的特定结构 或 重头预测。文献中给出的软件很多：http://www.nature.com/hdy/journal/v104/n6/fig_tab/hdy2009165t1.html#figure-title。

1.2.1 基于library-based的软件

library-based的软件，需要构建library，该library中包含很多来自不同物种某一重复序列的一致性序列，然后通过相似性比对来鉴定repeats。这种方法能对所有的种类的重复序列进行鉴定。此方法最经典最流行的软件是RepeatMasker；此方法中CENSOR和MASKERAID两个软件可以用于改良RepeatMasker的结果；此外，用于基因组的重复序列鉴定的还有GREEDIER（Li et al.,2008），该软件在其文章中表明该软件性能还不错，在repeats鉴定的敏感性上稍微比RepeatMasker高一点，但是repeats的鉴定率只有RepeatMasker的一半左右.

1.2.2 基于signature的软件

基于signature的方法主要用于TEs的鉴定。
1. 鉴定LTR逆转座子: LTR_STRUC (McCarthy and McDonald, 2003), LTR_PAR (Kalyanaraman and Aluru, 2006), FIND_LTR (Rho et al., 2007), RETROTECTOR (Sperber et al., 2007), LTR_FINDER (Xu and Wang, 2007) and LTRHARVEST (Ellinghaus et al., 2008)。文献中，这些软件中LTR_STRUC的敏感性最高(96%)，但是LTR的鉴定率只有30%；而LTRharvest的鉴定率最高(42%),敏感性67%.总体上，作者依次推荐的软件是LTRHARVEST和FIND_LTR(敏感性83%，鉴定率37%)。
2. 鉴定non-LTR retrotransposons: TSDFINDER (Szak et al., 2002), SINEDR (Tu et al., 2004) and RTANALYZER (Lucier et al., 2007)。其中，第一个软件用于验证RepeatMasker检测到的L1 insertions；第二个软件用于检测侧翼有TSDs（target site duplications 当重复序列插入到基因组上时，其两侧会带入短核酸序列的重复）的SINEs；第三个软件通过一些特征，比如TSDs，polyA尾和5’端核酸内切酶位点等来通过打分来检测L1逆转座子。
3. 鉴定MITEs:FINDMITE (Tu, 2001), TRANSPO (Santiago et al., 2002), MITE Analysis Kit (MAK) (Yang and Hall, 2003) and MITE Uncovering SysTem (MUST) (Chen et al., 2009)。文献中作者使用第一个软件报错，使用第三个软件却下载不到。第二个软件不能寻找新的MITEs，看来最好是使用最新的第四个软件。
4. 鉴定helitrons: HELITRONFINDER。该软件(Du et al., 2008)用来寻找玉米基因组中的helitrons（在动物和植物中有发现）。

1.2.3 重头预测的软件

1.2.3.1 自我比较的方法

通过BLAST、PALS等方法，将序列和自身进行比较，从而找出重复序列。软件有：REPEAT PATTERN TOOLKIT (Agarwal and States, 1994), RECON (Bao and Eddy, 2002), PILER (Edgar and Myers, 2005) and the BLASTER suite (used in Quesneville et al., 2005).其中RECON软件使用最广泛。

1.2.3.2 k-mer and spaced seed approaches

一定长度的k-mer出现了多次，可以鉴定为重复序列；spaced seed则是k-mer的一种衍生，seed上允许有一定的差异。

软件有：REPUTER (Kurtz and Schleiermacher, 1999), VMATCH (Kurtz, unpublished), REPEAT-MATCH (Delcher et al., 1999), MER-ENGINE (Healy et al., 2003), FORREPEATS (Lefebvre et al., 2003), REAS (Li et al., 2005), REPEATSCOUT (Price et al., 2005), RAP (Campagna et al., 2005), REPSEEK (Achaz et al., 2007), TALLYMER (Kurtz et al., 2008) and P-CLOUDS (Gu et al., 2008).

1.2.4 其它重复序列鉴定软件

其它一些用于鉴定非TEs的软件：Tandem Repeats Finder (TRF) (Benson, 1999), Tandem Repeat Occurrence Locator (TROLL) (Castelo et al., 2002), MREPS (Kolpakov et al., 2003), TRAP (Sobreira et al., 2006) and Optimized Moving Window Spectral Analysis (OMWSA) (Du et al., 2007) have been developed specifically to detect tandem repeats. The Inverted Repeat Finder (IRF) program (Warburton et al., 2004) was designed to search for inverted repeats.

1.3 多个软件整合的pipeline程序

REPEATMODELER pipeline (Smit, unpublished http://www.repeatmasker.org/RepeatModeler.html) includes the programs RECON, REPEATSCOUT, REPEATMASKER and TRF. It uses the output of the RECON and REPEATSCOUT programs to build, refine and classify consensus models of putative interspersed repeats.

当然，在此文献中，也讲述了其它很多专门用途的其它pipeline软件。而REPEATMODELER pipeline是现在运用于基因组的重复序列鉴定最主流的软件。以下将讲述该软件运用。

2. RepeatModeler的安装与使用

2.1 软件的安装

RepeatMasker和RepeatModeler是ISB（Institute for Systems Biology)的软件。ISB is located in the South Lake Union neighborhood。

根据RepeatMasker说明，其安装与使用需要：Perl 5.8.0以上版本，序列比对Engine，TRF和Repeat Database。其中序列比对engine可以安装多个，但每次只能使用其一，可以使用Cross_match,RMBlast,HMMER和ABBlaast/WUBlast等。

根据RepeatMideker说明，其安装与使用需要：Perl 5.8.8以上版本，RepeatMasker,Repeat Database,RECON,RepeatScout,TRF和序列比对engine。其中序列比对engine可以安装多个，但每次只能使用其一，可以使用RMBlast和ABBlaast/WUBlast。

再安装完毕需要的软件后，对RepeatMasker和RepeatModeler进行configure的时候填入相应软件的路径即可安装。

2.2 RepeatModeler的使用

2.2.1 使用RepeatModeler来通过基因组序列构建library

$ $RepeatModelerHome/BuildDatabase -name species \
  -engine ncbi species.genome.fasta
$ $RepeatModelerHome/RepeatModeler -database species

结果生成了一个文件夹，名称为RM_[PID].[DATE] ie. “RM_5098.MonMar141305172005″。该文件夹中的”consensi.fa.classified”即为library，用于RepeatMasker的输入。值得注意的是，可能fasta文件需要序列间不能有空（空格和换行等），否则会程序出错。

2.2.2 使用RepeatMasker来进行重复序列掩盖和重复序列计算

cp RM_5098.MonMar141305172005/consensi.fa.classified .
mkdir Repeat_result
$ $RepeatMaskerHome/RepeatMasker -pa 8 \
  -e ncbi -lib consensi.fa.classified \
  -dir Repeat_result -gff species.genome.fasta

则生成的结果文件位于Repeat_result文件夹中。

1. GapFiller简介

组装出来的基因组草图的scaffold需要进一步进行gaps的关闭。进行这样功能的软件有：SOAPdenovo GapCloser v1.12r6; IMAGE; GapFiller.

GapFiller文章发表在Genome Biology上：Boetzer M,Pirovano W. 2012. Toward almost closed genomes with GapFiller. Genome Biol.13:R56。从文章可以完全明白该软件closing gap的原理。

GapFiller需要输入scaffold序列(FASTA)和NGS paired-read数据(FASTA or FASTAQ)，输出FASTA格式文件。该软件的获得需要填写一些邮箱和单位信息。商业license需要花钱；学术性需要引用其文章。

2. GapFiller安装

下载GapFiller的安装包，解压缩后，里面包含bowtie、bwa和example共3个文件，其最重要的是GapFiller.pl文件，为主程序。还有2个PDF格式的manual文件。

3. GapFiller的使用

直接运行主程序，会给出软件的参数说明，如下：

-l library文件
-s scaffold序列的fasta文件
-m default:29 和gap边缘重叠的最小碱基数,该数值最好设置比reads的长度小一点点的数。比如36bp长度的reads，设置该值为30~35.
-o default:2 在补洞时，延伸一个碱基最小需要的reads数.
-r default:0.7 在补洞时，至少有该比例reads的碱基一致，才能对该碱基位点进行延伸。
-d default:50 gap部分序列的允许的最大差异。填补gap后，若值“填补上的序列长度 - gap长度”大于该阈值，则停止补洞；若小于该阈值，则不进行融合。
-n default:10 在一个scaffold中对邻近的两个contigs进行融合所需要最小重叠的碱基数。
-t default:10 由于gap边缘的碱基大部分是低质量碱基，补洞时需要先将gap边缘该数目的碱基trim掉，作为N处理。
-i default:10 迭代的最大次数。
-g default:1 使用bowtie进行比对的时候允许的最大的gap数，和bowtie中的-v参数一致
-T default:1 运行时使用的线程数
-S 跳过重新读取输入文件
-b 输出文件的basename。

-l 参数所指向的library文件需要先行编辑好。该文件包含7列，每一列之间以空格(space)隔开.其例子和格式如下：

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 bowtie file3.1.fastq file3.2.fastq 4000 0.5 RF

第1列：library名称
第2列：使用的序列比对方法，如果reads长度<50,则使用bowtie；若长度>50并<150,则使用bwa；若长度很大，比如454的reads，则使用bwa。BWA和BWA-sw运行在默认模式下。 第3,4列：双末端测序的fastq文件或fasta文件。 第5,6列：插入片段的长度，以及承认的长度。比如上例子中插入片段长度为400bp，成对的reads的片段长度只有在[400-400*0.25,400+400*0.25]范围内才被承认。 第7列：双端测序reads的方向，有FF，FR，RF和RR几种。

4. 例子

编辑一个libraries.txt文件，内容如下：

Illumina_160bp bwa fragment.reads1.fastq fragment.reads2.fastq 156 0.25 FR
Illumina_6000bp bwa jumping.reads1.fastq jumping.reads2.fastq 6170 0.25 FR

运行GapFiller程序，如下：

$GapFillerHome/GapFiller.pl -l libraries.txt\
  -s genome.fasta -m 90 -T 8 -b species

1. PASA简介

PASA, acronym for Program to Assemble Spliced Alignments, is a eukaryotic genome annotation tool that exploits spliced alignments of expressed transcript sequences to automatically model gene structures, and to maintain gene structure annotation consistent with the most recently available experimental sequence data. PASA also identifies and classifies all splicing variations supported by the transcript alignments.

Note:
Combine genome and Trinity de novo RNA-Seq assemblies to generate a comprehensive transcript database.

2. PASA使用前的准备

2.1 Mysql数据库的准备

创建只读权限用户和所有权限用户各一个。

mysql> GRANT SELECT ON *.* TO 'pasa'@'%' IDENTIFIED BY '123456';
mysql> GRANT ALL ON *.* TO 'chenlianfu'@'%' IDENTIFIED BY '123456';
mysql> FLUSH PRIVILEGES;

2.1 安装perl模块

# cpan
cpan[1]> install DBD::mysql
cpan[1]> install GD

2.3 安装GMAP

$ wget http://research-pub.gene.com/gmap/src/gmap-gsnap-2013-03-31.v5.tar.gz
$ tar zxvf gmap-gsnap-2013-03-31.v5.tar.gz
$ cd gmap-2013-03-31
$ ./configure --prefix=$PWD
$ make -j 8
$ make install

2.4 安装BLAT

$ wget http://hgwdev.cse.ucsc.edu/~kent/src/blatSrc35.zip
$ unzip blatSrc35.zip
$ cd blatSrc
$ MACHTYP=x86_64
$ export MACHTYPE
$ mkdir -p ~/bin/x86_64
$ make -j 8

2.5 安装FASTA

$ wget http://faculty.virginia.edu/wrpearson/fasta/fasta3/CURRENT.tar.gz
$ tar zxvf CURRENT.tar.gz
$ cd fasta-35.4.12
$ cd src
$ make -f ../make/Makefile.linux_sse2 all
$ cd ..
$ ln -s $PWD/bin/fasta35 ~/bin/fasta

2.6 安装PASA

$ wget http://kaz.dl.sourceforge.net/project/pasa/PASA2-r20130425beta.tgz
$ tar zxvf PASA2-r20130425beta.tgz
$ cd PASA2-r20130425beta/
$ make -j 8

2.7 安装GD

安装GD需要先行安装libgd

$ wget https://bitbucket.org/libgd/gd-libgd/get/93368566388c.zip
$ unzip 93368566388c.zip
$ cd libgd-gd-libgd-93368566388c
$ ./bootstrap.sh
$ ./configure
$ make -j 8
$ sudo make install
$ gdlib-config

再安装GD

$ wget http://search.cpan.org/CPAN/authors/id/L/LD/LDS/GD-2.49.tar.gz
$ tar zxvf GD-2.49.tar.gz
$ cd GD-2.49
$ perl Makefile.PL
$ make -j 8
$ sudo make install

安装GD的目的是能通过网页来查看PASA的运行结果。

2.8 配置PASA

2.8.1. 修改PASA的配置文件$PASAHOME/pasa_conf/conf.txt

$ cp $PASAHOME/pasa_conf/pasa.CONFIG.template $PASAHOME/pasa_conf/conf.txt
$ vim $PASAHOME/pasa_conf/conf.txt

2.8.2. 该文件需要修改的地方：

PASA_ADMIN_EMAIL=(your email address)
MYSQLSERVER=(your mysql server name)   此处不能填写IP。
MYSQL_RO_USER=(mysql read-only username)
MYSQL_RO_PASSWORD=(mysql read-only password)
MYSQL_RW_USER=(mysql all privileges username)
MYSQL_RW_PASSWORD=(mysql all privileges password)
BASE_PASA_URL=http://server_name/pasa/cgi-bin/

2.8.3. 修改httpd配置文件，

# vim /etc/httpd/conf/httpd.conf
# /etc/init.d/httpd restart

在/etc/httpd/conf/httpd.conf添加如下几行:

ScriptAlias /pasa "$PASAHOME"
<Directory "$PASAHOME">
        Options MultiViews ExecCGI
        AllowOverride None
        Order allow,deny
        Allow from all
</Directory>

2.9 cleaning the transcript sequences[Optional, requires seqclean to be installed

下载两个污染数据库，为fasta文件。

$ cd $PASAHOME/seqclean
$ tar zxf seqclean.tar.gz
$ cd seqclean
$ wget ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec -O UniVec.fasta
$ wget ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec_Core -O UniVec_Core.fasta

UniVec_Core includes only oligonucleotides and vectors consisting of bacterial, phage, viral, yeast or synthetic sequences. Vectors that include sequences of mammalian origin are excluded.

3. PASA主程序的使用

PASA的主程序是： $PASAHOME/scripts/Launch_PASA_pipeline.pl, 其使用参数如下：

*代表该参数是必须的

-c <filename> *
比对配置文件。可以将$PASAHOME/pasa_conf/pasa.alignAssembly.Template.
txt复制过来，只是将其中的MYSQLDB修改成需要的mysql数据库名。

####################

spliced alignment settings:
--ALIGNERS <string>
比对的软件，可用的软件有gmap和blat。也可以同时选择使用'gmap,blat'

-N <int> default: 1
max number of top scoring alignments

--MAX_INTRON_LENGTH | -I <int>  default: 100000
max intron length parameter passed to GMAP or BLAT

--IMPORT_CUSTOM_ALIGNMENTS_GFF3 <filename>
only using the alignments supplied in the corresponding GFF3 file.

--cufflinks_gtf <filename>
incorporate cufflinks-generated transcripts

####################

actions
-C
    flag, create MYSQL database
-R
    flag, run alignment/assembly pipeline.
-A
    compare to annotated genes.
--ALT_SPLICE
    flag, run alternative splicing analysis

-R 用于比对transcripts ， -A 用于和已有gff3注释文件的比较和更新；这两个参数不
能同时共用，使用不同的参数，则 -C 参数设置不同的参数文件。

####################

input files

-g <filename> *
    genome sequence FASTA file

-t <filename> *
    transcript db

-f <filename>
    file containing a list of fl-cdna accessions.

--TDN <filename>
    file containing a list of accessions corresponding to Trinity
 (full) de novo assemblies (not genome-guided)

####################

polyAdenylation site identification  ** highly recommended **
-T
    flag,transcript db were trimmed using the TGI seqclean tool.
-u <filename>
    value, transcript db containing untrimmed sequences (input to 
seqclean).a filename with a .cln extension should also exist, gen
erated by seqclean.

####################

Jump-starting or prematurely terminating
-x
    flag, print cmds only, don't process anything. (useful to get 
indices for -x or -e opts below)
-s <int>
    pipeline index to start running at (avoid rerunning searches).
-e <int>
    pipeline index where to stop running, and do not execute this 
entry. 

####################

Misc:
--TRANSDECODER
    flag, run transdecoder to identify candidate full-length coding
 transcripts
--CPU <int> default: 2
    multithreading
-d  flag, Debug 
-h  flag, print this option menu and quit

install List::Util FindBin YAML Storable Getopt::Long Algorithm::Diff Data::Dumper Module::Pluggable GraphViz

1. 了解超几何概率分布

参考：The Hypergeometric Distribution
超几何分布的公式为：

p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)
for x = 0, …, k.
其中, m 是桶里面白球的个数, n 是黑球的个数, k 是从桶中随机取出的球数, x 是取出球
中白球的个数.

Fisher‘s Exact Test就是依据超几何概率分布得到的。

2. 超几何分布运算方法

使用R的运算方法：

dhyper(x, m, n, k, log = FALSE)
    计算某一个点的p值
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
    计算一个分布的p值，默认下计算P(X<=x)
qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
    计算某一个p值对应的分位数
rhyper(nn, m, n, k)
    按超几何分布给出nn的可能的模拟结果值

3. 对一组p值进行校正

参考：Adjust P-values for Multiple Comparisons
使用R的运算方法：

p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c(“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,
# “fdr”, “none”)

4. 利用上述原理来编写perl程序，来进行GO或PATHWAY富集分析

输入两个文件：第1个文件是各个基因对应的GO分类或KEGG分类；第2个文件是差异表达基因的名称。以下为KEGG富集分析的perl程序。

 #!/usr/bin/perl

=head1 Name

  kegg_enrich_analysis.pl

=head1 Description

  This program is designed to do enrichment analysis by kegg results.

=head1 Version

  Author: Chen Lianfu, chenllianfu@foxmail.com
  Version: 1.0,  Date: Sat Apr 27 02:06:02 2013

=head1 Usage

  --out         the output result name, default: enrichment_output
  --tmp         to keep the tmp files, default: false
  --fdr         set the false discovery rate, default: 0.05
  --kdn         When a kegg term was perpared for Fisher's exact 
test, the i(the num of this kegg terms have so many differential 
expression genes) should >= this value, default: 3
  --help        output help information to screen

=head1 Example

  perl ../kegg_enrich_analysis.pl annot_keggs genelist.txt

=cut

use strict;
use Getopt::Long;

my ($out,$tmp,$fdr_threshold,$kegg2degs_num_threshold,$Help);

GetOptions(
    "out:s"=>\$out,
    "fdr:s"=>\$fdr_threshold,
    "kdn:s"=>\$kegg2degs_num_threshold,
    "tmp"=>\$tmp,
    "help"=>\$Help
);

$out ||= 'enrichment_output';
$fdr_threshold ||= 0.05;
$kegg2degs_num_threshold ||= 3;

die `pod2text $0` if (@ARGV == 0 || $Help);

my $kegg2gene_file = shift @ARGV;
my $deglist = shift @ARGV;

open kegg2GENE, '<', $kegg2gene_file or die $!;
open DEGLIST, '<', $deglist or die $!;

my %kegg2genes;
my %gene2keggs;
while ( <kegg2GENE> ) {
    next unless /^(\S+).*?(K\d+)/;
    $gene2keggs{$1} .= "$2,";
    $kegg2genes{$2}{"genes"} .= "$1,";
    $kegg2genes{$2}{"num"} += 1;
}

my ( $total_deg_num, $deg_own_kegg_num );
my %deg_kegg2genes;
while ( <DEGLIST> ) {
    chomp;
    my $genename = $_;
    $total_deg_num ++;
    if ( $gene2keggs{$genename} ) {
        $deg_own_kegg_num ++;
        my @keggs = split /,/, $gene2keggs{$genename};

        foreach ( @keggs ) {
            $deg_kegg2genes{$_}{"genes"} .= "$genename,";
            $deg_kegg2genes{$_}{"num"} ++;
        }
    }
}

print "$total_deg_num differential expression genes in total, and 
$deg_own_kegg_num genes own kegg terms.\n";

my @for_enrich_keggs = keys %deg_kegg2genes;
open FOR_R, '>', "tmp_4value" or die $!;

my $kegg_num_for_fisher_test = @for_enrich_keggs;
print "$kegg_num_for_fisher_test kegg terms is perpared for Fisher
's exact test.\n";

foreach my $kegg ( @for_enrich_keggs ) {
    my $n_and_m = keys %gene2keggs;
    my $N = $deg_own_kegg_num;
    my $n = $kegg2genes{$kegg}{"num"};
    my $i = $deg_kegg2genes{$kegg}{"num"};
    my $m = $n_and_m -$n;

    print FOR_R "$kegg $i $n $m $N\n" if $i >= $kegg2degs_num_threshold;
}

open R_SCRIPT, '>', "phyper_adjust.R" or die $!;

print R_SCRIPT 
'phy <- read.table(file="tmp_4value")
pvalue <- phyper(phy$V2,phy$V3,phy$V4,phy$V5,lower.tail=FALSE)
qvalue <- p.adjust(pvalue,method=\'fdr\')
value <- cbind ( pvalue, qvalue )
rownames(value)=phy$V1
value
';

my @fdr = `cat phyper_adjust.R | R --vanilla --slave`;

#print @fdr;

`rm phyper_adjust.R tmp_4value` unless $tmp;

open OUTPUT, '>', $out or die $!;

foreach ( @fdr ) {
    next unless /(\S+)\s+(\S+)\s+(\S+)/;
    my $kegg = $1; my $pvalue = $2; my $fdr = $3;

    if ( $fdr <= $fdr_threshold ) {
        my $n_and_m = keys %gene2keggs;
        my $N = $deg_own_kegg_num;
        my $n = $kegg2genes{$kegg}{"num"};
        my $i = $deg_kegg2genes{$kegg}{"num"};
        my $m = $n_and_m -$n;

        my $genes = $deg_kegg2genes{$kegg}{"genes"};

        print OUTPUT "$kegg\t$i\t$N\t$n\t$n_and_m\t$pvalue\t$fdr\t$genes\n";
    }
}

Reference: The Hypergeometric Distribution

Description

Density, distribution function, quantile function and random generation for the hypergeometric distribution.

Usage

dhyper(x, m, n, k, log = FALSE)
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
rhyper(nn, m, n, k)

dhper    计算某一个点的p值
phyper   计算一个分布的p值，默认下计算P(X<=x)
qhyper   计算某一个p值对应的分位数
rhyper   按超几何分布给出nn的可能的模拟结果值

Arguments

x, q vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.

m the number of white balls in the urn.

n the number of black balls in the urn.

k the number of balls drawn from the urn.

p probability, it must be between 0 and 1.

nn number of observations. If length(nn) > 1, the length is taken to be the number required.

log, log.p logical; if TRUE, probabilities p are given as log(p).

lower.tail logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].

Details

The hypergeometric distribution is used for sampling without replacement. The density of this distribution with parameters m, n and k (named Np, N-Np, and n, respectively in the reference below) is given by

p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)
for x = 0, …, k.

The quantile is defined as the smallest value x such that F(x) ≥ p, where F is the distribution function.

Value

dhyper gives the density, phyper gives the distribution function, qhyper gives the quantile function, and rhyper generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rhyper, and is the maximum of the lengths of the numerical parameters for the other functions.

The numerical parameters other than n are recycled to the length of the result. Only the first elements of the logical parameters are used.

Source

dhyper computes via binomial probabilities, using code contributed by Catherine Loader (see dbinom).

phyper is based on calculating dhyper and phyper(...)/dhyper(...) (as a summation), based on ideas of Ian Smith and Morten Welinder.

qhyper is based on inversion.

rhyper is based on a corrected version of

Kachitvichyanukul, V. and Schmeiser, B. (1985). Computer generation of hypergeometric random variates. Journal of Statistical Computation and Simulation, 22, 127–145.

References

Johnson, N. L., Kotz, S., and Kemp, A. W. (1992) Univariate Discrete Distributions, Second Edition. New York: Wiley.
See Also
Distributions for other standard distributions.

Examples

m x rbind(phyper(x, m, n, k), dhyper(x, m, n, k))
all(phyper(x, m, n, k) == cumsum(dhyper(x, m, n, k))) # FALSE
## but error is very small:
signif(phyper(x, m, n, k) - cumsum(dhyper(x, m, n, k)), digits = 3)