结论
LDhat:不能用来计算大样品的重组率。18个样品,200位点耗时大概一小时;52样品,200个位点,耗时遥遥无期。
CPLDhat:改进后的LDhat,可以使用mpirun进行多线程计算,大大节省运行时间;52样品,200位点已运行5天
Phase:单次只能计算500个位点,速度慢;52样品,500个位点,耗时11天。
Mlrho:不能计算两两位点之间的重组率,只能计算全基因组的重组率,运行速度很快。
LDhat
命令行示例
$LDhat_DIR/interval -seq chrA01_1.site -loc chrA01_1.locus -lk Bra.lk -its 1000000 -bpen 5 -samp 2000
input file
- LOCUS文件
- SITE文件
- Lookup文件
CPLDhat
命令行示例
$DIR/mpirun -n 30 $DIR/interval -seq chrA01_3.site -loc chrA01_3.locus -lk Bra.lk -convergence
input file
- CPLDhat 输入文件格式与LDhat完全一致
- CPLDhat 不需要设置参数
- CPLDhat 只是LDhat的基础上稍微修改,使其能用mpirun多线程运行,进而提升节省运行时间
Phase
Usage
$DIR/PHASE -MR -X10 chrA01_1.inp chrA01_1.out
*input file *
Phase只有一个输入文件,其格式如下:
3 ## 3个样品
5 ## 5个位点
P 300 1313 1500 2023 5635 ##P固定格式,后跟5个位点的具体位置
MSSSM ##M代表该位点为多等位基因,S代表该位点二等位基因。1、5 为多等位位点,2、3、4为二等位位点
#sample1 ##样品1名字
12 1 0 1 3 ##For multiallelic loci a positive integermust be used for each allele (representing the number of repeats atmicrosatellite loci),and data for each locus should be separated by aspace
11 0 1 0 3 ##For biallelic loci, any two characters (e.g.A/C, G/T, 0/1) can be used to represent the two alleles, and theydo not need to be separated by a space
#sample2
12 1 1 1 2
12 0 0 0 3
#sample3
-1 ? 0 0 2 ##Missing alleles at multiallelic loci should be representedby -1.
-1 ? 1 1 13 ##Missing alleles at SNPloci should be entered as ?
out put
计算重组率的结果在_recom文件里
_freqs Estimates of the sample haplotype frequencies (which can alsobe used as estimates of the population haplotype frequencies).The first two columns are self-explanatory. The next 2 columnsgive i) estimates (posterior means) of haplotype frequencies forthe whole sample, and ii) estimated standard deviations (squareroot of the variance of the posterior distribution) for these frequencies. If the -c option is used (see below) then additionalpairs of columns give estimated haplotype frequencies and standard deviations for each group specified in the input file.
_pairs List of the most likely pairs of haplotypes for each individual,together with their probability. Only pairs whose probabilityexceeds the threshold given by the -O flag are listed.
_recom Contains estimates of recombination parameters across the region, using the general model for varying recombination ratefrom Li and Stephens (2003). (Note: these estimates are produced only if the recombination model is used: see the -M option.Also, the estimates have a straightforward interpretion only ifthe positions of the loci are specified in the input file.) Thefirst line of the file gives the positions of the loci in the inputfile (to facilitate subsequent analyses using results in this file).Each subsequent line of the file gives a sample from the posterior distribution of the recombination parameters. The firstcolumn gives estimates of the background recombination rate(more precisely, the value of the population genetics parameterρ, here denoted ρ ¯). Column 2 gives the factor by which the recombination rate between locus 1 and 2 exceeds the background rate, and, similarly, column i gives the factor by which the ratebetween locus i − 1 and locus i exceeds the background rate.To get point estimates of these quantities I suggest taking themedian of the results for each column. If you are planning tomake use of the results from this file, then it is advisable to usethe -X option described later (eg -X10 or -X100), to obtain moreaccurate estimates.
_hotspot Contains estimates of recombination parameters if one of thehotspot options is used: see the -M option. The estimates have astraightforward interpretion only if the positions of the loci arespecified in the input file. Each line of the file gives a samplefrom the posterior distribution of the recombination parameters.The first column gives estimates of the background recombination rate (more precisely, the value of the population geneticsparameter ρ ¯). Column 2 gives the estimated left hand edge ofthe first hotspot, column 3 gives the right-hand edge, and column 4 gives the estimated intensity (the factor by which therecombination rate in the hotspot exceeds the background rate.)If this last column is 1 then this corresponds to no increase inrecombination rate in the “hotspot”. To get point estimates ofthese quantities I suggest examining both the mean and the median of the results for each column. (Histograms of each columncan also be helpful in giving you an idea of the uncertainty andskewness in the posterior.) If you are planning to make use ofthe results from this file, then it is advisable to use the -X option described later (eg -X10 or -X100) to obtain more accurateestimates.
_signif Gives the p-value for a permutation test of the null hypothesisthat the cases and controls are random draws from a commonset of haplotype frequencies (only if the case-control status isspecified for each individual; see the -c option below)).
_monitor The monitor file for the goodness of fit of the estimated haplotypes to the underlying model. It has two columns: the first isthe pseudo-likelihood from Stephens and Donnelly (2003); thesecond is the PAC-B likelihood from Li and Stephens (2003).When comparing behaviour of different runs from different starting points I recommend using the first column, as the secondcolumn can be very dependent on the order of the individuals, which can vary a lot across runs. This file replaces the‘‘temp.monitor’’ file produced by previous versions of PHASE
Mlrho
usage
#format
/gpfs01/duliuwen/BioSoft/MlRho_2.8/Auxiliary_Software/FormatPro_0.5/formatPro Bra.pro
#Compute the genome-wide population mutation rate
/gpfs01/duliuwen/BioSoft/MlRho_2.8/mlRho -M 0 -I
#The -I switch in the last command instructed mlRho to write the likelihoods to the binary file profileDb.lik,see likelihood file
/gpfs01/duliuwen/BioSoft/MlRho_2.8/Auxiliary_Software/InspectPro_0.3/inspectPro profileDb.lik |head -n 10
#To compute ∆ and ρ as a function of ∆, use
/gpfs01/duliuwen/BioSoft/MlRho_2.8/mlRho -m 100 -M 200 >rho
input file
Bra.pro格式示例
>Contig1 ##contig1名字
1 0 0 0 5 ##第一列位点位置,第二列A的个数,第二列C的个数,第三列G的个数,第四列T的个数
2 1 10 0 0 ## 此文件由bam文件转换而来,官网有格式转换的辅助软件
>Contig2
5 2 0 7 0
6 12 0 0 0
8 0 10 1 0