This is a note to self: I can never remember how to install an R package on the Sanger computer system. Here goes:
1. Change directory to ~alc/R/library.
2. Start R.
2. Download the R library eg. for the 'bblme2' package, by typing in R:
> system("wget cran.r-project.org/src/contrib/bbmle_1.0.5.2.tar.gz")
3. Install the R package by typing in R:
> install.packages("bbmle_1.0.5.2.tar.gz", repos=NULL, lib="~alc/R/library")
4. Now you should be able to use the package in R:
> library(bbmle, lib="~alc/R/library")
[Note to self: Magdalena has installed R packages compiled for R3-0.0 here:
/nfs/users/nfs_m/mz3/bin/R-packages ]
Monday, 15 April 2013
Wednesday, 10 April 2013
Making maximum likelihood estimates of parameters using R
Making a maximum likelihood estimate of a binomial probability
To make a maximum likelihood estimate of a binomial probability you can use the mle2() function in the 'bbmle' R package.
For example, if you roll a particular (six-sided) die 10 times, and observe seven '5's, then you can estimate the probability of getting a '5' when you roll that die. Here our data is simply that we observe seven '5's:
> mydata <- c(7)
We just have one observation in our data set here.
Next we need to construct a negative log-likelihood function, as the mle2() R function (which we will use to calculate the maximum likelihood estimate, see below) requires a negative log-likelihood function as input, rather than a likelihood function.
Our negative log likelihood function will be minus the log of the probability of observing seven '5's out of 10 throws, according to a binomial distribution B(10, p) with a particular value of the binomial probability, p='prob':
> myfunc <- function(size,prob) { -sum(dbinom(mydata,size,prob,log=TRUE)) }
Then we can call the 'mle2' function in the 'bbmle' package to get the maximum likelihood estimate of the probability of observing a '5' with our die, given that we observed seven '5's on 10 throws of the die:
> library('bbmle')
> mle2(myfunc, start=list(prob=0.5), data=list(size=10))
Coefficients:
prob
0.7000032
Log-likelihood: -1.32
This tells us that the maximum likelihood estimate of the binomial probability is 0.7. That is, this is our estimate of the probability of observing a '5' with our die. This makes sense, as we did observe '5's on 7 out of 10 throws.
The log likelihood is given as -1.32. This is the log of the probability of observing seven '5's on 10 throws of a die that has p=0.7:
> log(dbinom(7, size=10, prob=0.7))
[1] -1.321151
We can also make a plot of the likelihood function L(p) plotted against the possible values of the binomial probability p (Note: the 'curve' function requires the function in terms of x, so we use x here in R instead of p):
> curve(dbinom(7, 10, x))
We see that the peak is at about 0.7.
Note I have plotted the binomial probability mass function (p.m.f.) at each possible value of p, this is the likeihood function. Note that in 'myfunc' above, we used the negative log likeilhood function (because mle2() requires a negative log likelihood function rather than a likelihood function), which is -log of the binomial probability mass function:
> curve(-log(dbinom(7,10,x)))
The negative log likelihood function has its minimum where the likelihood function has its maximum, that is, at p = 0.7.
Making a maximum likelihood estimate of a geometric parameter
A geometric distribution is indexed by one parameter, p. Say we observe three values, 3, 4, 8, that come from a geometric distribution with unknown value of the parameter p. These are the numbers of tries up to and including the first success.
The dgeom function in R assumes the data is in the form of the number of failures before the first success. Thus, the values given to dgeom should be 2, 3, and 7.
We can make a maximum likelihood estimate of the value of p by typing:
> mydata <- c(2, 3, 7)
> myfunc <- function(prob) { -sum(dgeom(mydata,prob,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(prob=0.5))
Coefficients:
prob
0.2000023
Log-likelihood: -7.51
This tells us that the maximum likelihood estimate of p is about 0.2.
We can also make a plot of the likelihood function L(p) plotted against the possible values of the geometric parameter p:
> xvalues <- seq(0.01,0.99,by=0.01)
> myfunc2 <- function(prob) { prod(dgeom(mydata,prob)) }
> yvalues <- sapply(xvalues, myfunc2)
> plot(xvalues,yvalues,type="l")
The peak of the plot is at 0.2. This makes sense as the average of 3, 4, and 8 is 5, and 1/5 is 0.2.
Note I have plotted the products of the geometric probability mass function (p.m.f.) for the data at each possible value of p, which is the likelihood function. In 'myfunc' above, we used the negative log likelihood function, which is just the negative of the sum of the logs of the geometric probability mass functions.
Making a maximum likelihood estimate of a Poisson parameter
A Poisson distribution is indexed by one parameter, lambda. Say we have observed counts of some event with the following frequencies:
Count Frequency
0 58
1 25
2 13
3 2
4 2
5 1
6 1
7 0
8 1
>=9 0
Suppose the counts are observations from a Poisson distribution with mean lambda. We can make a maximum likelihood estimate of the value of lambda by typing:
> mydata <- c(rep(0,58),rep(1,25),rep(2,13),rep(3,2),rep(4,2),rep(5,1),rep(6,1),rep(8,1))
> myfunc <- function(lambda) { -sum(dpois(mydata,lambda,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(lambda=0.5))
lambda
0.8155337
Log-likelihood: -142.05
That is, the maximum likelihood estimate of lambda is 0.8155337, that is, about 0.816. This makes sense, as the sample mean is about 0.816:
> mean(mydata)
[1] 0.815534
Making a maximum likelihood estimate of an exponential parameter
An exponential distribution has one parameter, theta. Say we observe the waiting times between successive occurrences of some event (eg. a town flooding) to be 5 years, 8 years, and 7 years. Suppose the time between occurrences of the event are observations from an exponential distribution with parameter theta. We can make a maximum likelihood estimate of theta by typing:
> mydata <- c(5,8,7)
> myfunc <- function(x) { -sum(dexp(mydata,x,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(x=0.5))
Coefficients:
x
0.15
Log-likelihood: -8.69
That is, a maximum likelihood estimate of the value of theta is 0.15.
Making a maximum likelihood estimate of a Rayleigh parameter
The probability density function (p.d.f.) of a Rayleigh distribution is:
f(x; theta) = (x/theta^2) * exp(-x^2/(2*theta^2)), where theta > 0
That is, the distribution has one parameter, theta. Say we observe a random sample of six observations of X: 22.2, 2.8, 4.0, 13.9, 11.7 and 8.3. We can make a maximum likelihood estimate of theta as follows:
> mypdf <- function(x, theta)
{
ifelse(theta > 0, (x/(theta^2)) * exp(-(x*x)/(2*theta*theta)), NA)
}
> mydata <- c(22.2, 2.8, 4.0, 13.9, 11.7, 8.3)
> myfunc <- function(theta)
{
mylen <- length(mydata)
myprod <- mypdf(mydata[1], theta)
for (i in 2:mylen) { myprod <- myprod * mypdf(mydata[i], theta) }
mylogL <- -log(myprod)
return(mylogL)
}
Here the function mypdf() is for calculating the p.d.f., and the function myfunc() is for calculating the negative of the log likelihood. Now calculate the maximum likelihood estimate of theta:
> library('bbmle')
> mle2(myfunc, start=list(theta=8.5))
8.735013
Log-likelihood: -19.28
This tells us that the maximum likelihood estimate of theta is about 8.735.
Making a maximum likelihood estimate of a parameter in a model
Say we observe four different types of objects, 187 of object 1, 35 of object 2, 37 of object 3, and 31 of object 4.
And say we have some theory (model) that says that the four types of objects should occur with probabilities (9/16) + p, (3/16) - p, (3/16) - p, and (1/16) + p, respectively, where the value of parameter p is unknown.
Therefore, the likelihood function is:
L(p) = (((9/16) + p)^187) * (((3/16) - p)^35) * (((3/16) - p)^37) * (((1/16) + p)^31)
= (((9/16) + p)^187) * (((3/16) - p)^72) * (((1/16) + p)^31)
We can plot this function in R by typing (note that we use x instead of p in R, as the curve() function expects a function in terms of x):
> curve( (((9/16)+x)^187) * (((3/16)-x)^72) * (((1/16)+x)^31) )
We see that the function has a vertical asymptote at p=1. However, suppose we have some prior information that the value of p that we are interested in is just between 0 and 0.1. We can plot the function for this range:
> curve( (((9/16)+x)^187) * (((3/16)-x)^72) * (((1/16)+x)^31) , from = 0, to = 0.1)
We see that there is a local maximum at about 0.058 or 0.059.
To estimate the maximum likelihood value of p, we need to define the negative log likelihood function for just the range 0 to 0.1:
> myfunc <- function(p)
{
ifelse((0 <= p & p <= 0.1),-log((((9/16)+p)^187) * (((3/16)-p)^72) * (((1/16)+p)^31)),NA)
}
> mle2(myfunc, start=list(p=0.05))
Coefficients:
p
0.05838241
Log-likelihood: -302.01
This tells us that the maximum likelihood estimate of p is 0.05838241, that is, about 0.0584.
Note that if we try to estimate the maximum likelihood estimate of p by using the function:
-log((((9/16)+p)^187) * (((3/16)-p)^72) * (((1/16)+p)^31))
without restricting the range of p to 0-0.1, the mle2() function will give us an error message, because of the vertical asymptote at p=1.
Furthermore, if we give mle2() a starting value that is too far from 0.05838241, it will give an error message, eg.:
> mle2(myfunc, start=list(p=0))
Error in optim(par = 0, fn = function (p) :
non-finite finite-difference value [1]
To make a maximum likelihood estimate of a binomial probability you can use the mle2() function in the 'bbmle' R package.
For example, if you roll a particular (six-sided) die 10 times, and observe seven '5's, then you can estimate the probability of getting a '5' when you roll that die. Here our data is simply that we observe seven '5's:
> mydata <- c(7)
We just have one observation in our data set here.
Next we need to construct a negative log-likelihood function, as the mle2() R function (which we will use to calculate the maximum likelihood estimate, see below) requires a negative log-likelihood function as input, rather than a likelihood function.
Our negative log likelihood function will be minus the log of the probability of observing seven '5's out of 10 throws, according to a binomial distribution B(10, p) with a particular value of the binomial probability, p='prob':
> myfunc <- function(size,prob) { -sum(dbinom(mydata,size,prob,log=TRUE)) }
Then we can call the 'mle2' function in the 'bbmle' package to get the maximum likelihood estimate of the probability of observing a '5' with our die, given that we observed seven '5's on 10 throws of the die:
> library('bbmle')
> mle2(myfunc, start=list(prob=0.5), data=list(size=10))
Coefficients:
prob
0.7000032
Log-likelihood: -1.32
This tells us that the maximum likelihood estimate of the binomial probability is 0.7. That is, this is our estimate of the probability of observing a '5' with our die. This makes sense, as we did observe '5's on 7 out of 10 throws.
The log likelihood is given as -1.32. This is the log of the probability of observing seven '5's on 10 throws of a die that has p=0.7:
> log(dbinom(7, size=10, prob=0.7))
[1] -1.321151
We can also make a plot of the likelihood function L(p) plotted against the possible values of the binomial probability p (Note: the 'curve' function requires the function in terms of x, so we use x here in R instead of p):
> curve(dbinom(7, 10, x))
We see that the peak is at about 0.7.
Note I have plotted the binomial probability mass function (p.m.f.) at each possible value of p, this is the likeihood function. Note that in 'myfunc' above, we used the negative log likeilhood function (because mle2() requires a negative log likelihood function rather than a likelihood function), which is -log of the binomial probability mass function:
> curve(-log(dbinom(7,10,x)))
The negative log likelihood function has its minimum where the likelihood function has its maximum, that is, at p = 0.7.
Making a maximum likelihood estimate of a geometric parameter
A geometric distribution is indexed by one parameter, p. Say we observe three values, 3, 4, 8, that come from a geometric distribution with unknown value of the parameter p. These are the numbers of tries up to and including the first success.
The dgeom function in R assumes the data is in the form of the number of failures before the first success. Thus, the values given to dgeom should be 2, 3, and 7.
We can make a maximum likelihood estimate of the value of p by typing:
> mydata <- c(2, 3, 7)
> myfunc <- function(prob) { -sum(dgeom(mydata,prob,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(prob=0.5))
Coefficients:
prob
0.2000023
Log-likelihood: -7.51
This tells us that the maximum likelihood estimate of p is about 0.2.
We can also make a plot of the likelihood function L(p) plotted against the possible values of the geometric parameter p:
> xvalues <- seq(0.01,0.99,by=0.01)
> myfunc2 <- function(prob) { prod(dgeom(mydata,prob)) }
> yvalues <- sapply(xvalues, myfunc2)
> plot(xvalues,yvalues,type="l")
The peak of the plot is at 0.2. This makes sense as the average of 3, 4, and 8 is 5, and 1/5 is 0.2.
Note I have plotted the products of the geometric probability mass function (p.m.f.) for the data at each possible value of p, which is the likelihood function. In 'myfunc' above, we used the negative log likelihood function, which is just the negative of the sum of the logs of the geometric probability mass functions.
Making a maximum likelihood estimate of a Poisson parameter
A Poisson distribution is indexed by one parameter, lambda. Say we have observed counts of some event with the following frequencies:
Count Frequency
0 58
1 25
2 13
3 2
4 2
5 1
6 1
7 0
8 1
>=9 0
Suppose the counts are observations from a Poisson distribution with mean lambda. We can make a maximum likelihood estimate of the value of lambda by typing:
> mydata <- c(rep(0,58),rep(1,25),rep(2,13),rep(3,2),rep(4,2),rep(5,1),rep(6,1),rep(8,1))
> myfunc <- function(lambda) { -sum(dpois(mydata,lambda,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(lambda=0.5))
lambda
0.8155337
Log-likelihood: -142.05
That is, the maximum likelihood estimate of lambda is 0.8155337, that is, about 0.816. This makes sense, as the sample mean is about 0.816:
> mean(mydata)
[1] 0.815534
Making a maximum likelihood estimate of an exponential parameter
An exponential distribution has one parameter, theta. Say we observe the waiting times between successive occurrences of some event (eg. a town flooding) to be 5 years, 8 years, and 7 years. Suppose the time between occurrences of the event are observations from an exponential distribution with parameter theta. We can make a maximum likelihood estimate of theta by typing:
> mydata <- c(5,8,7)
> myfunc <- function(x) { -sum(dexp(mydata,x,log=TRUE)) }
> library('bbmle')
> mle2(myfunc, start=list(x=0.5))
Coefficients:
x
0.15
Log-likelihood: -8.69
That is, a maximum likelihood estimate of the value of theta is 0.15.
Making a maximum likelihood estimate of a Rayleigh parameter
The probability density function (p.d.f.) of a Rayleigh distribution is:
f(x; theta) = (x/theta^2) * exp(-x^2/(2*theta^2)), where theta > 0
That is, the distribution has one parameter, theta. Say we observe a random sample of six observations of X: 22.2, 2.8, 4.0, 13.9, 11.7 and 8.3. We can make a maximum likelihood estimate of theta as follows:
> mypdf <- function(x, theta)
{
ifelse(theta > 0, (x/(theta^2)) * exp(-(x*x)/(2*theta*theta)), NA)
}
> mydata <- c(22.2, 2.8, 4.0, 13.9, 11.7, 8.3)
> myfunc <- function(theta)
{
mylen <- length(mydata)
myprod <- mypdf(mydata[1], theta)
for (i in 2:mylen) { myprod <- myprod * mypdf(mydata[i], theta) }
mylogL <- -log(myprod)
return(mylogL)
}
Here the function mypdf() is for calculating the p.d.f., and the function myfunc() is for calculating the negative of the log likelihood. Now calculate the maximum likelihood estimate of theta:
> library('bbmle')
> mle2(myfunc, start=list(theta=8.5))
8.735013
Log-likelihood: -19.28
This tells us that the maximum likelihood estimate of theta is about 8.735.
Making a maximum likelihood estimate of a parameter in a model
Say we observe four different types of objects, 187 of object 1, 35 of object 2, 37 of object 3, and 31 of object 4.
And say we have some theory (model) that says that the four types of objects should occur with probabilities (9/16) + p, (3/16) - p, (3/16) - p, and (1/16) + p, respectively, where the value of parameter p is unknown.
Therefore, the likelihood function is:
L(p) = (((9/16) + p)^187) * (((3/16) - p)^35) * (((3/16) - p)^37) * (((1/16) + p)^31)
= (((9/16) + p)^187) * (((3/16) - p)^72) * (((1/16) + p)^31)
We can plot this function in R by typing (note that we use x instead of p in R, as the curve() function expects a function in terms of x):
> curve( (((9/16)+x)^187) * (((3/16)-x)^72) * (((1/16)+x)^31) )
We see that the function has a vertical asymptote at p=1. However, suppose we have some prior information that the value of p that we are interested in is just between 0 and 0.1. We can plot the function for this range:
> curve( (((9/16)+x)^187) * (((3/16)-x)^72) * (((1/16)+x)^31) , from = 0, to = 0.1)
We see that there is a local maximum at about 0.058 or 0.059.
To estimate the maximum likelihood value of p, we need to define the negative log likelihood function for just the range 0 to 0.1:
> myfunc <- function(p)
{
ifelse((0 <= p & p <= 0.1),-log((((9/16)+p)^187) * (((3/16)-p)^72) * (((1/16)+p)^31)),NA)
}
> mle2(myfunc, start=list(p=0.05))
Coefficients:
p
0.05838241
Log-likelihood: -302.01
This tells us that the maximum likelihood estimate of p is 0.05838241, that is, about 0.0584.
Note that if we try to estimate the maximum likelihood estimate of p by using the function:
-log((((9/16)+p)^187) * (((3/16)-p)^72) * (((1/16)+p)^31))
without restricting the range of p to 0-0.1, the mle2() function will give us an error message, because of the vertical asymptote at p=1.
Furthermore, if we give mle2() a starting value that is too far from 0.05838241, it will give an error message, eg.:
> mle2(myfunc, start=list(p=0))
Error in optim(par = 0, fn = function (p) :
non-finite finite-difference value [1]
Friday, 5 April 2013
Running Interproscan on a gene set
A useful way to annotate a genome is to run InterProScan on your gene set.
[Note: this is really a note-to-self, as it's only useful to Sanger users.]
Note: some of the below has recently become obsolete, so I've greyed it out.
Using Andrew Page's annotate_eukaryotes script:
At Sanger, you can run interproscan by using the annotate_eukaryotes script by Andrew Page, for example:
% annotate_eukaryotes -a Sratti_v5.2.genes_prot.fa -o Sratti_v5.interpro.gff
where Sratti_v5.2.genes_prot.fa is your input fasta file of proteins, Sratti_v5.interpro.gff is your interproscan output file.
Note: the Interproscan results give GO terms associated with the proteins. My colleague Bernardo Foth has a script for making a file of the GO terms for proteins (/nfs/users/nfs_b/bf3/bin/interproscanResults_parse_GOterms_01.pl).
Using iprscan_chado:
At Sanger, you can do this using the iprscan_chado:
1) Log into a screen window, eg. on farm2-head4:
% screen -RD
2) Edit the .iprscan file in your home directory (eg. /nfs/users/nfs_a/alc/.iprscan) so that it has your team for submitting farm jobs (eg. team133) eg.
LSB_DEFAULTPROJECT=team133
export LSB_DEFAULTPROJECT
3) Start the InterProScan job using:
% iprscan_chado -f <fasta>
where <fasta> is your fasta file
eg.
% iprscan_chado -f Sratti_v4.genes_prot.fa
This will submit your jobs to the farm, there is no need to use 'bsub'.
The results will appear in file <fasta>.interpro, eg. Sratti_v4.genes_prot.fa.interpro
Note: I found that the iprscan_chado program crashes if you have '|' symbols in your gene names, so you will need to rename your genes, and try to run it again if this is the case.
Using Magdalena's script to break up the input file before running InterProScan:
1) Log into pcs4
2) Figure out how many bytes are in your input fasta file (using 'ls -al'). If it is a soft-link, make sure you take the size of the original file.
3) Then work out how many bytes you need per file in order to get ~1000 sequences per file:
(bytes_full_file * 1000)/seqs_full_file
where bytes_full_file is the number of bytes in the input fasta file (from step 2), seqs_full file is the number of sequences in the input fasta file.
This will give us a number of bytes_per_file.
4) Use Magdalena's script to break up the input fasta file into files of about bytes_per_file bytes each:
% perl ~mz3/bin/perl/interpro_scan_splitter.pl input_fasta prefix bytes_per_file
where input_fasta is your input fasta file, prefix is the prefix to give to each of the smaller fasta files, bytes_per_file is the number of bytes each of the smaller fasta files should be.
eg.
% perl ~mz3/bin/perl/interpro_scan_splitter.pl Sratti_v5.genes_prot.fa SrattiV5 490000
This will make many smaller files (usually ~10-20) with some of the sequences from your input fasta file in each smaller file. The smaller files will be called prefix_1, prefix_2, etc., eg. SrattiV5_1, SrattiV5_2, etc.
5) Open a screen session. If the smaller files made in step (4) are called SrattiV5_1, SrattiV5_2, etc. then in the first screen window type:
% iprscan_chado -f SrattiV5_1
In the second screen window type:
% iprscan_chado -f SrattiV5_2
In the third screen window type:
% iprscan_chado -f SrattiV5_3
and so on, until you have started iprscan_chado running for each of the smaller files made in step (4).
Using Magdalena's script to run pfamscan:
% perl -w /nfs/users/nfs_m/mz3/bin/perl/pfamscan_splitter.pl Sratti_v5.genes_prot.fa srattipfam 490000
This will make many smaller files (usually ~10-20) with some of the sequences from your input fasta file in each smaller file. The smaller files will be called prefix_1, prefix_2, etc., eg. srattipfam_1, srattipfam_2, etc.
The pfam jobs are submitted to the farm automatically.
Note added 29-Apr-2016: This doesn't seem to work anymore. See my blog on pfam_scan.pl instead.
Difference between interproscan and pfamscan
Interproscan is more comprehensive as it includes scanning pfam with other things, and it will tell you things such as how many domains of a certain type there are in a protein.
Thanks to Magdalena Zarowiecki for this.
[Note: this is really a note-to-self, as it's only useful to Sanger users.]
Note: some of the below has recently become obsolete, so I've greyed it out.
Using Andrew Page's annotate_eukaryotes script:
At Sanger, you can run interproscan by using the annotate_eukaryotes script by Andrew Page, for example:
% annotate_eukaryotes -a Sratti_v5.2.genes_prot.fa -o Sratti_v5.interpro.gff
where Sratti_v5.2.genes_prot.fa is your input fasta file of proteins, Sratti_v5.interpro.gff is your interproscan output file.
Note: the Interproscan results give GO terms associated with the proteins. My colleague Bernardo Foth has a script for making a file of the GO terms for proteins (/nfs/users/nfs_b/bf3/bin/interproscanResults_parse_GOterms_01.pl).
Using iprscan_chado:
At Sanger, you can do this using the iprscan_chado:
1) Log into a screen window, eg. on farm2-head4:
% screen -RD
2) Edit the .iprscan file in your home directory (eg. /nfs/users/nfs_a/alc/.iprscan) so that it has your team for submitting farm jobs (eg. team133) eg.
LSB_DEFAULTPROJECT=team133
export LSB_DEFAULTPROJECT
3) Start the InterProScan job using:
% iprscan_chado -f <fasta>
where <fasta> is your fasta file
eg.
% iprscan_chado -f Sratti_v4.genes_prot.fa
This will submit your jobs to the farm, there is no need to use 'bsub'.
The results will appear in file <fasta>.interpro, eg. Sratti_v4.genes_prot.fa.interpro
Note: I found that the iprscan_chado program crashes if you have '|' symbols in your gene names, so you will need to rename your genes, and try to run it again if this is the case.
Using Magdalena's script to break up the input file before running InterProScan:
1) Log into pcs4
2) Figure out how many bytes are in your input fasta file (using 'ls -al'). If it is a soft-link, make sure you take the size of the original file.
3) Then work out how many bytes you need per file in order to get ~1000 sequences per file:
(bytes_full_file * 1000)/seqs_full_file
where bytes_full_file is the number of bytes in the input fasta file (from step 2), seqs_full file is the number of sequences in the input fasta file.
This will give us a number of bytes_per_file.
4) Use Magdalena's script to break up the input fasta file into files of about bytes_per_file bytes each:
% perl ~mz3/bin/perl/interpro_scan_splitter.pl input_fasta prefix bytes_per_file
where input_fasta is your input fasta file, prefix is the prefix to give to each of the smaller fasta files, bytes_per_file is the number of bytes each of the smaller fasta files should be.
eg.
% perl ~mz3/bin/perl/interpro_scan_splitter.pl Sratti_v5.genes_prot.fa SrattiV5 490000
This will make many smaller files (usually ~10-20) with some of the sequences from your input fasta file in each smaller file. The smaller files will be called prefix_1, prefix_2, etc., eg. SrattiV5_1, SrattiV5_2, etc.
5) Open a screen session. If the smaller files made in step (4) are called SrattiV5_1, SrattiV5_2, etc. then in the first screen window type:
% iprscan_chado -f SrattiV5_1
In the second screen window type:
% iprscan_chado -f SrattiV5_2
In the third screen window type:
% iprscan_chado -f SrattiV5_3
and so on, until you have started iprscan_chado running for each of the smaller files made in step (4).
Using Magdalena's script to run pfamscan:
% perl -w /nfs/users/nfs_m/mz3/bin/perl/pfamscan_splitter.pl Sratti_v5.genes_prot.fa srattipfam 490000
This will make many smaller files (usually ~10-20) with some of the sequences from your input fasta file in each smaller file. The smaller files will be called prefix_1, prefix_2, etc., eg. srattipfam_1, srattipfam_2, etc.
The pfam jobs are submitted to the farm automatically.
Note added 29-Apr-2016: This doesn't seem to work anymore. See my blog on pfam_scan.pl instead.
Difference between interproscan and pfamscan
Interproscan is more comprehensive as it includes scanning pfam with other things, and it will tell you things such as how many domains of a certain type there are in a protein.
Thanks to Magdalena Zarowiecki for this.
Friday, 29 March 2013
Selecting a probability model for your data
A common task is to choose a statistical model given some data. Here are some general ideas on how to start to think about what is the most appropriate statistical model:
Continuous data eg. measurements, time periods, very large numbers of counts
Examples of continuous distributions are exponential, continuous uniform, and Normal distributions.
Measurement data eg. heights of people; weights of people: a Normal distribution might be appropriate. Characteristics of data for which a Normal model may be appropriate are: the distribution has range from -Infinity to Infinity; is symmetric around a single mode that coincides with the mean; and values that are far from the mean are unlikely.
Very large counts of objects observed eg. number of lottery tickets sold per week; number of fish caught each week by European fishermen: a Normal distribution might be appropriate. You can use a Normal probability plot to help you decide whether a Normal model is appropriate.
Waiting times between successive events eg. time intervals between buses passing a particular bus stop; time intervals between earthquakes: an exponential distribution might be appropriate. An assumption of an exponential model is that events occur at random in time. Characteristics of data for which an exponential distribution may be appropriate are: the distribution has range from 0 to Infinity; is right-skewed (has skewness of 2); has its mode at 0; and the mean and standard deviation of an exponential distribution are equal. You can use an exponential probability plot to help you decide whether an exponential model is appropriate.
Which particular value occurs, out of all possible values in a particular interval (a, b): a continuous uniform distribution might be appropriate. An assumption of a continuous uniform distribution is that each value in the interval is believed to be equally likely. Characteristics of data for which a continuous uniform distribution may be appropriate are: the distribution has finite range from a to b, and has no mode.
Discrete data eg. small numbers of counts
Examples of discrete distributions are Bernoulli (really a special case of the binomial with n = 1), binomial, discrete uniform, genometric and Poisson distributions.
Counts of objects or events observed in a fixed interval of time/space eg. number of goals in each of 100 different soccer matches; number of fish caught in each of 100 different nets: a Poisson distribution might be appropriate. An assumption of a Poisson(mu) distribution is that events occur at random. Characteristics of data for which a Poisson model may be appropriate are: the distribution has an unbounded range {0,1,2...}; is right-skewed (if skewed at all); has one mode; and the mean and variance of a Poisson distribution are equal.
Number of successes in n trials, eg. number of sixes that you get in 100 throws of a die; number of faulty light-bulbs out of 10,000 produced by a factory; number of 1000 people suffering a particular disease who recover when given a particular drug: a Binomial distribution might be appropriate. Assumptions of a B(n, p) model are that the n trials are independent, and that there is a constant probability p of success at each trial. Characteristics of data for which a Binomial model may be appropriate are: the distribution has a finite range {0...n}; and has just one mode.
Number of trials up to and including the first success, eg. number of times you throw a die before you get the first six: a Geometric distribution might be appropriate. Assumptions of a G(p) distribution are that you have a sequence of independent trials, and that there is a constant probability of success p between trials. Characteristics of data for which a Geometric model may be appropriate are: the distribution has an unbounded range {1,2,3...}; is right skewed; and has one mode at 1.
Which particular outcome happens out of a set of several equally likely outcomes, eg. getting 1, 2, 3, 4, 5, or 6 when you throw a die: a discrete uniform distribution might be appropriate. Assumptions of a discrete uniform distribution are that there is a finite set of outcomes possible, and that every outcome is believed to be equally likely. Characteristics of data for which a discrete uniform model may be appropriate are: the distribution has a finite range, and no mode.
Continuous data eg. measurements, time periods, very large numbers of counts
Examples of continuous distributions are exponential, continuous uniform, and Normal distributions.
Measurement data eg. heights of people; weights of people: a Normal distribution might be appropriate. Characteristics of data for which a Normal model may be appropriate are: the distribution has range from -Infinity to Infinity; is symmetric around a single mode that coincides with the mean; and values that are far from the mean are unlikely.
Very large counts of objects observed eg. number of lottery tickets sold per week; number of fish caught each week by European fishermen: a Normal distribution might be appropriate. You can use a Normal probability plot to help you decide whether a Normal model is appropriate.
Waiting times between successive events eg. time intervals between buses passing a particular bus stop; time intervals between earthquakes: an exponential distribution might be appropriate. An assumption of an exponential model is that events occur at random in time. Characteristics of data for which an exponential distribution may be appropriate are: the distribution has range from 0 to Infinity; is right-skewed (has skewness of 2); has its mode at 0; and the mean and standard deviation of an exponential distribution are equal. You can use an exponential probability plot to help you decide whether an exponential model is appropriate.
Which particular value occurs, out of all possible values in a particular interval (a, b): a continuous uniform distribution might be appropriate. An assumption of a continuous uniform distribution is that each value in the interval is believed to be equally likely. Characteristics of data for which a continuous uniform distribution may be appropriate are: the distribution has finite range from a to b, and has no mode.
Discrete data eg. small numbers of counts
Examples of discrete distributions are Bernoulli (really a special case of the binomial with n = 1), binomial, discrete uniform, genometric and Poisson distributions.
Counts of objects or events observed in a fixed interval of time/space eg. number of goals in each of 100 different soccer matches; number of fish caught in each of 100 different nets: a Poisson distribution might be appropriate. An assumption of a Poisson(mu) distribution is that events occur at random. Characteristics of data for which a Poisson model may be appropriate are: the distribution has an unbounded range {0,1,2...}; is right-skewed (if skewed at all); has one mode; and the mean and variance of a Poisson distribution are equal.
Number of successes in n trials, eg. number of sixes that you get in 100 throws of a die; number of faulty light-bulbs out of 10,000 produced by a factory; number of 1000 people suffering a particular disease who recover when given a particular drug: a Binomial distribution might be appropriate. Assumptions of a B(n, p) model are that the n trials are independent, and that there is a constant probability p of success at each trial. Characteristics of data for which a Binomial model may be appropriate are: the distribution has a finite range {0...n}; and has just one mode.
Number of trials up to and including the first success, eg. number of times you throw a die before you get the first six: a Geometric distribution might be appropriate. Assumptions of a G(p) distribution are that you have a sequence of independent trials, and that there is a constant probability of success p between trials. Characteristics of data for which a Geometric model may be appropriate are: the distribution has an unbounded range {1,2,3...}; is right skewed; and has one mode at 1.
Which particular outcome happens out of a set of several equally likely outcomes, eg. getting 1, 2, 3, 4, 5, or 6 when you throw a die: a discrete uniform distribution might be appropriate. Assumptions of a discrete uniform distribution are that there is a finite set of outcomes possible, and that every outcome is believed to be equally likely. Characteristics of data for which a discrete uniform model may be appropriate are: the distribution has a finite range, and no mode.
Wednesday, 27 March 2013
Speeding up BLAST jobs
I have discussed with colleagues lately how to speed up some large BLASTX jobs. Here are several ideas that we came up with:
1.) Reduce the database you are searching:
Just take a selection of representative species for the taxa of interest (eg. Bacteria), rather than all species. This will probably not miss many proteins. However, some species have species-specific genes or even strain-specific genes, especially in bacteria.
2.) Search a non-redundant database:
You can download the whole NCBI nr database from the NCBI ftp site NCBI ftp site. However, it is pretty big (14 Gbyte unzipped). Also, it has all species in it, so it can be slow to extract a subset of sequences (ie. idea (1) above). Instead you can download RefSeq sequences for a particular species or set of species, by going to the NCBI Protein website, and setting 'Limits' to Field=Organism and Source Database=RefSeq, and then searching for a particular species (eg. 'Escherchia coli') in the search box.
3). Use a larger BLAST word size:
The word size can be set in BLAST using the -W option. By default it is set to 3 for proteins. If you set it to 5, it will speed up your BLAST search, although with some loss of sensitivity.
4). Stripe your files in a /lustre filesystem:
If you are using a /lustre filesystem, it will speed up your BLAST search if you stripe the files in the directory. You can do this by typing:
% lfs setstripe <dir> -c -1
where <dir> is the directory containing your BLAST database.
This is a good idea as you will probably be running many BLAST jobs against this single database at the same time. (By striping a large database file across all OSTs, it maximises the IO bandwidth. As the file is large, the overhead in opening the striped file is negligible compared to the time it takes to read the data.)
5). Use megablast instead of BLASTN:
If you are trying to speed up BLAST searches between very similar sequences, Megablast is faster than normal BLASTN.
1.) Reduce the database you are searching:
Just take a selection of representative species for the taxa of interest (eg. Bacteria), rather than all species. This will probably not miss many proteins. However, some species have species-specific genes or even strain-specific genes, especially in bacteria.
2.) Search a non-redundant database:
You can download the whole NCBI nr database from the NCBI ftp site NCBI ftp site. However, it is pretty big (14 Gbyte unzipped). Also, it has all species in it, so it can be slow to extract a subset of sequences (ie. idea (1) above). Instead you can download RefSeq sequences for a particular species or set of species, by going to the NCBI Protein website, and setting 'Limits' to Field=Organism and Source Database=RefSeq, and then searching for a particular species (eg. 'Escherchia coli') in the search box.
3). Use a larger BLAST word size:
The word size can be set in BLAST using the -W option. By default it is set to 3 for proteins. If you set it to 5, it will speed up your BLAST search, although with some loss of sensitivity.
4). Stripe your files in a /lustre filesystem:
If you are using a /lustre filesystem, it will speed up your BLAST search if you stripe the files in the directory. You can do this by typing:
% lfs setstripe <dir> -c -1
where <dir> is the directory containing your BLAST database.
This is a good idea as you will probably be running many BLAST jobs against this single database at the same time. (By striping a large database file across all OSTs, it maximises the IO bandwidth. As the file is large, the overhead in opening the striped file is negligible compared to the time it takes to read the data.)
5). Use megablast instead of BLASTN:
If you are trying to speed up BLAST searches between very similar sequences, Megablast is faster than normal BLASTN.
Monday, 25 March 2013
Using CEGMA to assess genome assemblies
The CEGMA software (Parra, Bradnam & Korf, 2007) can be used to assess the completeness of genome assemblies for eukaryotic species. CEGMA defines a set of 458 conserved protein families that occur in a wide range of eukaryotes, based on a subset of the KOGs database.
This subset consists of KOGs families that contain at least one protein from six selected species (human, Drosophila melanogaster, Arabidopsis thaliana, Caenorhabditis elegans, Saccharomyces cerevisiae, Schizosaccharomyces pombe), which when re-aligned using t-coffee give an alignment that passes several criteria (all proteins must cover at least 75% of the length of the alignment; each protein must have only 5 internal gaps longer than 10 amino acids; the average percent identity over all rows of the alignment must be >10%).
They use an approach combining geneid, GeneWise and TBLASTN searches to predict genes for each of these 458 protein families in the genome assembly of interest. First, TBLASTN is used to identify candidate regions in the genome assembly. Only the 5 best candidate regions are considered further.
Next, a Hidden Markov Model (HMM) is built for each of the 458 protein families using HMMER, and gene predictions are made with GeneWise using these HMMs in the corresponding candidate regions of the genome. The GeneWise alignments are then given to the geneid program to aid geneid in making gene predictions in those candidate regions. They find that the combined use of geneid and GenWise in this way produces more accurate predictions compared to using GeneWise alone.
A filter is applied to the geneid predictions to determine if they are similar enough to the rest of the KOG protein family to be considered orthologs. To do this, the predicted proteins (from the gene predictions) are aligned against the HMM defived from the corresponding protein family. Only predictions that have a strong match to the HMM of the gene family are kept (compared to the match that existing members of the KOG family have to the HMM).
As mentioned above, 5 candidate regions are considered for each of the 458 conserved KOG families. Thus, it is possible at this stage that we have geneid predictions in more than two candidate regions for a particular KOG family. The geneid prediction with the highest scoring match to the HMM for the family is conserved the true ortholog of the KOG family, and any other predictions are assumed to be paralogs.
The last step is to train geneid using the initial set of predictions, so that geneid will produce species-specific coding and splice site models. The version of geneid that has now been trained for your particular assembly is now re-run on the (up to) five candidate regions for each of the 458 KOG families, to make a final set of gene predictions.
To assess how complete your genome is, you can count the number of the 458 KOG families that CEGMA manages to make gene predictions for. You would expect that if the assembly is fairly complete, it should manage to make gene predictions for all 458 KOG families.
The CEGMA paper reports that CEGMA should find gene predictions for all 458 families if they are indeed present in the assembly, but in some rare cases it happens that the gene is present, but CEGMA fails to find it because it is too diverged.
Intermediate files produced by CEGMA steps
(i) step 1: running TBLAST to find the top 5 BLAST matches to the CEGMA family in the genome:
I think this produces files genome.blast and genome.blast.gff
(ii) step 2: uses a HMM for the CEGMA family to make GeneWise predictions in each of the 5 BLAST hit regions:
I think this produces file local.genewise.gff
(iii) step 3: the GeneWise predictions are given to Geneid to make gene predictions in each of the 5 BLAST hit regions:
I think this produces file local.geneid.gff
Note: the Geneid predictions are filtered, and only those that are very similar to the rest of the CEGMA protein family are kept. Only Geneid predictions that have a strong match to the HMM of the CEGMA family are kept, compared to the match that existing members of the CEGMA family have to the HMM.
(iv) step 4: among the Geneid predictions made in each of the top 5 BLAST hit region (multiple predictions per region), only the Geneid prediction (in a region) with the highest scoring match to the HMM of the CEGMA family is kept.
I think this produces file local.geneid.selected.gff [alternatively it might be local_self.geneid.gff, I'm not sure].
(v) step 5: Geneid is trained using the initial set of predictions for all CEGMA families (to train splice sites, etc.), and then run again in the 5 top BLAST hit regions for each family, to make the final CEGMA gene set for the assembly.
I think this produces output.cegma.gff.
How to run CEGMA
To run CEGMA, you can type:
% /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma --tmp --ext -g <genome>
where /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma is where you have installed CEGMA,
--tmp means that temporary files are not deleted,
--ext gives you extended output files with more information,
-g <genome> specifies your assembly fasta file of scaffolds/chromosomes.
eg.
% /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma --tmp --ext -g 03.SS.velvet.scaffolds.fa
where '03.SS.velvet.scaffolds.fa' is the file of scaffolds.
I found that CEGMA dies and gives an error if the names of the sequences in your fasta file of scaffolds contain '|', so you will need to rename the sequences if this is the case.
To run CEGMA, you probably need about 150 Mb of memory for a ~40 Mbase assembly file.
Running CEGMA on farm2 and farm3
[Sanger users only] CEGMA v2 is installed on farm2. To run CEGMA on farm2, you need to tell CEGMA which BLAST, geneid, genewise, etc. to use. It is easiest to do this using Eleanor Stanley's wrappe
% ~es9/bin/runcegma_kog.sh -g genome.fa
where genome.fa is your genome fasta file. You have to bsub this command to the farm. It might need ~2000 Mbyte of RAM.
CEGMA v2.4 is installed on farm3. To run CEGMA on farm3, it's easiest to use a wrapper:
% run_cegma_v2.sh -g genome.fa
CEGMA and HMMER2/HMMER3
Note: in our group we use CEGMA version 2.0, which uses HMMER2 format HMMs. There is a newer CEGMA release, version 2.4, which uses HMMER3 format HMMs. We found that the results for CEGMA 2.0 seemed slightly better than those for v. 2.4, and it was a little easier to use if you want to give it your own file of HMMs.
It is possible to convert HMMER2 HMMs to HMMER3 HMMs using the 'hmmconvert' program. Also, if you want to search CEGMA HMMs for protein matches, in HMMER3 hmmpfam has become hmmscan (searches a protein sequences against a database of HMMs).
CEGMA output
CEGMA gives you an output report that is called 'output.completeness_report'. This contains a summary of which of the subset of the 248 most highly-conserved CEGMA KOGs are present (either partially or completely). On this page, it explains that the 248 are 'likely to be found in low number of inparalogs in a wide range of species' as they are 'cases when only one ortholog in at least four of [the] six species'.
It will say something like this:
#Prots %Completeness - #Total Average %Ortho
Complete 223 89.92 - 237 1.06 5.83
Partial 247 99.60 - 275 1.11 10.93
This means that, for 248 highly conserved CEGMA KOGs, 237 (89.92%) are found as full predictions. Ideally, you would like a genome assembly to have all (100%) of these highly-conserved CEGMA KOGs.
Similarly, for 248 highly conserved CEGMA KOGs, 247 (99.6%) are found as full or partial predictions. This is of course higher than the percent that are found as full predictions.
The 'Average' column gives the 'average CEG gene number'. This tells us the average number of full (or partial) predictions made per CEGMA KOG. Here there are 1.06 full predictions made per CEGMA KOG, and 1.11 partial predictions.
Alternatives to CEGMA
An alternative to CEGMA for assessing genome completeness is BUSCOS (Benchmarking of Universal Single Copy Orthologs), which is mentioned briefly at the end of the OrthoDB paper. The data sets for BUSCOS are available here.
Finding out which CEGMA genes were found
To find out which of the 248 'CEGs' were found (and so count towards %completeness in the CEGMA output report), we can use the following:
% perl /nfs/users/nfs_m/mz3/bin/cegma_v2/src/completeness.pl -m hmm_select.aln /nfs/users/nfs_m/mz3/bin/cegma_v2/data/completeness_cutoff.tbl > cegma_MISSING.txt
(where hmm_select.aln is produced by CEGMA)
Thanks to James Cotton for telling me this!
On the CEGMA webpage, it suggests that you can test whether CEGMA has run properly on your computer system by using the sample.dna and sample.prot files that come with CEGMA. It says that you should get this output.
I have tried to test CEGMA on the Sanger farm2 and farm3 compute farms.
farm2:
farm2 has CEGMA v2 installed, which (according to the CEGMA README) requires wu_blast (which we have), hmmer 2.3.2 (which we have), geneid v1.2 (which we have), and GeneWise 2.2.3-rc7 (which we have). To run the test run of CEGMA v2, I typed:
% export CEGMA=/nfs/users/nfs_m/mz3/bin/cegma_v2
% export PATH=/software/pubseq/bin/wu_blast:$PATH
% export PATH=/nfs/users/nfs_m/mz3/bin/hmmer-2.3.2/bin:$PATH
% export PATH=/nfs/users/nfs_m/mz3/bin/geneid1.2/bin:$PATH
% export PERL5LIB=/nfs/users/nfs_m/mz3/bin/cegma_v2/lib:$PERL5LIB
% export PERL5LIB=/nfs/users/nfs_m/mz3/bin/cegma_v2/src:$PERL5LIB
% export CEGMATMP=$CEGMA_dir/temp
% $CEGMA/bin/cegma --tmp --ext -g /software/pathogen/external/apps/usr/local/cegma_v2.1.251109/sample/sample.dna /software/pathogen/external/apps/usr/local/cegma_v2.1.251109/sample/sample.prot
My output.completeness.report file says:
#Prots %Completeness - #Total Average %Ortho
Complete 6 2.42 - 6 1.00 0.00
Partial 6 2.42 - 6 1.00 0.00
farm3: [Note: this section is greyed out as I need to update it]
farm3 has CEGMA v2.4 installed, which (according to the CEGMA README), requires ncbi_blast 2.2.5 (which we have), hmmer 3.0 (which we have), genewise 2.2.3-rc7 (we have 2.4.1, but the CEGMA webpage says 2.4.1 is fine too), and geneid 2.4 (which we have). To run the test of CEGMA v2.4, I typed:
% export CEGMA=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312
% export PERL5LIB=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312/lib:$PERL5LIB
% export PERL5LIB=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312/src:$PERL5LIB
% export CEGMATMP=$CEGMA_dir/temp
% $CEGMA/bin/cegma --tmp --ext -g /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.dna /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.prot
Note I submitted this to the farm3 with 2000 Mbyte of memory, otherwise sometimes it ran out of memory and blast failed as a result (with a segmentation fault).
My run gave slightly different results than expected according to the CEGMA webpage, I got:
Found 2079 candidate regions in /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.dna
NOTE: created 87 geneid predictions
NOTE: Found 20 geneid predictions with scores above threshold value
DATA COLLECTED: 20 Coding sequences containing 57 introns
NOTE: created 46 geneid predictions
NOTE: Foud 20 geneid predictions with scores above threshold value
However, according to the CEGMA webpage, this is what is expected:
Found 86 candidate regions in /path/to/CEGMA/sample/sample.dna
NOTE: created 23 geneid predictions
NOTE: Found 15 geneid predictions with scores above threshold value
DATA COLLECTED: 15 Coding sequences containing 48 introns
NOTE: created 21 geneid predictions
NOTE: Foud 15 geneid predictions with scores above threshold value
My output.completeness.report file says:
#Prots %Completeness - #Total Average %Ortho
Complete 8 3.23 - 8 1.00 0.00
Partial 10 4.03 - 10 1.00 0.00
So it looks like a higher completeness was found using this test set than for CEGMA v2 on farm2 (see above).
I'm not sure why I get different results than on the CEGMA webpage. One difference is that I am using GeneWise 2.4.1 (not 2.2.3-rc7, the recommended one). However, a strange finding is that my output says that 2079 candidate regions were found using blast, while the CEGMA webpage says that 86 candidate regions were found, even though the blast versions are the same (blast+ 2.2.28).
I'm not sure about what difference using GeneWise 2.2.3-rc7 would make. I've tried to install this on farm3, but no luck so far compiling it..
Alternatives to CEGMA
An alternatives
Thanks
Thanks to Eleanor Stanley and Alan Tracey for help with CEGMA.
This subset consists of KOGs families that contain at least one protein from six selected species (human, Drosophila melanogaster, Arabidopsis thaliana, Caenorhabditis elegans, Saccharomyces cerevisiae, Schizosaccharomyces pombe), which when re-aligned using t-coffee give an alignment that passes several criteria (all proteins must cover at least 75% of the length of the alignment; each protein must have only 5 internal gaps longer than 10 amino acids; the average percent identity over all rows of the alignment must be >10%).
They use an approach combining geneid, GeneWise and TBLASTN searches to predict genes for each of these 458 protein families in the genome assembly of interest. First, TBLASTN is used to identify candidate regions in the genome assembly. Only the 5 best candidate regions are considered further.
Next, a Hidden Markov Model (HMM) is built for each of the 458 protein families using HMMER, and gene predictions are made with GeneWise using these HMMs in the corresponding candidate regions of the genome. The GeneWise alignments are then given to the geneid program to aid geneid in making gene predictions in those candidate regions. They find that the combined use of geneid and GenWise in this way produces more accurate predictions compared to using GeneWise alone.
A filter is applied to the geneid predictions to determine if they are similar enough to the rest of the KOG protein family to be considered orthologs. To do this, the predicted proteins (from the gene predictions) are aligned against the HMM defived from the corresponding protein family. Only predictions that have a strong match to the HMM of the gene family are kept (compared to the match that existing members of the KOG family have to the HMM).
As mentioned above, 5 candidate regions are considered for each of the 458 conserved KOG families. Thus, it is possible at this stage that we have geneid predictions in more than two candidate regions for a particular KOG family. The geneid prediction with the highest scoring match to the HMM for the family is conserved the true ortholog of the KOG family, and any other predictions are assumed to be paralogs.
The last step is to train geneid using the initial set of predictions, so that geneid will produce species-specific coding and splice site models. The version of geneid that has now been trained for your particular assembly is now re-run on the (up to) five candidate regions for each of the 458 KOG families, to make a final set of gene predictions.
To assess how complete your genome is, you can count the number of the 458 KOG families that CEGMA manages to make gene predictions for. You would expect that if the assembly is fairly complete, it should manage to make gene predictions for all 458 KOG families.
The CEGMA paper reports that CEGMA should find gene predictions for all 458 families if they are indeed present in the assembly, but in some rare cases it happens that the gene is present, but CEGMA fails to find it because it is too diverged.
Intermediate files produced by CEGMA steps
(i) step 1: running TBLAST to find the top 5 BLAST matches to the CEGMA family in the genome:
I think this produces files genome.blast and genome.blast.gff
(ii) step 2: uses a HMM for the CEGMA family to make GeneWise predictions in each of the 5 BLAST hit regions:
I think this produces file local.genewise.gff
(iii) step 3: the GeneWise predictions are given to Geneid to make gene predictions in each of the 5 BLAST hit regions:
I think this produces file local.geneid.gff
Note: the Geneid predictions are filtered, and only those that are very similar to the rest of the CEGMA protein family are kept. Only Geneid predictions that have a strong match to the HMM of the CEGMA family are kept, compared to the match that existing members of the CEGMA family have to the HMM.
(iv) step 4: among the Geneid predictions made in each of the top 5 BLAST hit region (multiple predictions per region), only the Geneid prediction (in a region) with the highest scoring match to the HMM of the CEGMA family is kept.
I think this produces file local.geneid.selected.gff [alternatively it might be local_self.geneid.gff, I'm not sure].
(v) step 5: Geneid is trained using the initial set of predictions for all CEGMA families (to train splice sites, etc.), and then run again in the 5 top BLAST hit regions for each family, to make the final CEGMA gene set for the assembly.
I think this produces output.cegma.gff.
How to run CEGMA
To run CEGMA, you can type:
% /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma --tmp --ext -g <genome>
where /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma is where you have installed CEGMA,
--tmp means that temporary files are not deleted,
--ext gives you extended output files with more information,
-g <genome> specifies your assembly fasta file of scaffolds/chromosomes.
eg.
% /nfs/users/nfs_m/mz3/bin/cegma_v2/bin/cegma --tmp --ext -g 03.SS.velvet.scaffolds.fa
where '03.SS.velvet.scaffolds.fa' is the file of scaffolds.
I found that CEGMA dies and gives an error if the names of the sequences in your fasta file of scaffolds contain '|', so you will need to rename the sequences if this is the case.
To run CEGMA, you probably need about 150 Mb of memory for a ~40 Mbase assembly file.
Running CEGMA on farm2 and farm3
[Sanger users only] CEGMA v2 is installed on farm2. To run CEGMA on farm2, you need to tell CEGMA which BLAST, geneid, genewise, etc. to use. It is easiest to do this using Eleanor Stanley's wrappe
% ~es9/bin/runcegma_kog.sh -g genome.fa
where genome.fa is your genome fasta file. You have to bsub this command to the farm. It might need ~2000 Mbyte of RAM.
CEGMA v2.4 is installed on farm3. To run CEGMA on farm3, it's easiest to use a wrapper:
% run_cegma_v2.sh -g genome.fa
CEGMA and HMMER2/HMMER3
Note: in our group we use CEGMA version 2.0, which uses HMMER2 format HMMs. There is a newer CEGMA release, version 2.4, which uses HMMER3 format HMMs. We found that the results for CEGMA 2.0 seemed slightly better than those for v. 2.4, and it was a little easier to use if you want to give it your own file of HMMs.
It is possible to convert HMMER2 HMMs to HMMER3 HMMs using the 'hmmconvert' program. Also, if you want to search CEGMA HMMs for protein matches, in HMMER3 hmmpfam has become hmmscan (searches a protein sequences against a database of HMMs).
CEGMA output
CEGMA gives you an output report that is called 'output.completeness_report'. This contains a summary of which of the subset of the 248 most highly-conserved CEGMA KOGs are present (either partially or completely). On this page, it explains that the 248 are 'likely to be found in low number of inparalogs in a wide range of species' as they are 'cases when only one ortholog in at least four of [the] six species'.
It will say something like this:
#Prots %Completeness - #Total Average %Ortho
Complete 223 89.92 - 237 1.06 5.83
Partial 247 99.60 - 275 1.11 10.93
This means that, for 248 highly conserved CEGMA KOGs, 237 (89.92%) are found as full predictions. Ideally, you would like a genome assembly to have all (100%) of these highly-conserved CEGMA KOGs.
Similarly, for 248 highly conserved CEGMA KOGs, 247 (99.6%) are found as full or partial predictions. This is of course higher than the percent that are found as full predictions.
The 'Average' column gives the 'average CEG gene number'. This tells us the average number of full (or partial) predictions made per CEGMA KOG. Here there are 1.06 full predictions made per CEGMA KOG, and 1.11 partial predictions.
Alternatives to CEGMA
An alternative to CEGMA for assessing genome completeness is BUSCOS (Benchmarking of Universal Single Copy Orthologs), which is mentioned briefly at the end of the OrthoDB paper. The data sets for BUSCOS are available here.
Finding out which CEGMA genes were found
To find out which of the 248 'CEGs' were found (and so count towards %completeness in the CEGMA output report), we can use the following:
% perl /nfs/users/nfs_m/mz3/bin/cegma_v2/src/completeness.pl -m hmm_select.aln /nfs/users/nfs_m/mz3/bin/cegma_v2/data/completeness_cutoff.tbl > cegma_MISSING.txt
(where hmm_select.aln is produced by CEGMA)
Thanks to James Cotton for telling me this!
Testing CEGMA on farm2 and farm3
On the CEGMA webpage, it suggests that you can test whether CEGMA has run properly on your computer system by using the sample.dna and sample.prot files that come with CEGMA. It says that you should get this output.
I have tried to test CEGMA on the Sanger farm2 and farm3 compute farms.
farm2:
farm2 has CEGMA v2 installed, which (according to the CEGMA README) requires wu_blast (which we have), hmmer 2.3.2 (which we have), geneid v1.2 (which we have), and GeneWise 2.2.3-rc7 (which we have). To run the test run of CEGMA v2, I typed:
% export CEGMA=/nfs/users/nfs_m/mz3/bin/cegma_v2
% export PATH=/software/pubseq/bin/wu_blast:$PATH
% export PATH=/nfs/users/nfs_m/mz3/bin/hmmer-2.3.2/bin:$PATH
% export PATH=/nfs/users/nfs_m/mz3/bin/geneid1.2/bin:$PATH
% export PERL5LIB=/nfs/users/nfs_m/mz3/bin/cegma_v2/lib:$PERL5LIB
% export PERL5LIB=/nfs/users/nfs_m/mz3/bin/cegma_v2/src:$PERL5LIB
% export CEGMATMP=$CEGMA_dir/temp
% $CEGMA/bin/cegma --tmp --ext -g /software/pathogen/external/apps/usr/local/cegma_v2.1.251109/sample/sample.dna /software/pathogen/external/apps/usr/local/cegma_v2.1.251109/sample/sample.prot
My output.completeness.report file says:
#Prots %Completeness - #Total Average %Ortho
Complete 6 2.42 - 6 1.00 0.00
Partial 6 2.42 - 6 1.00 0.00
farm3: [Note: this section is greyed out as I need to update it]
farm3 has CEGMA v2.4 installed, which (according to the CEGMA README), requires ncbi_blast 2.2.5 (which we have), hmmer 3.0 (which we have), genewise 2.2.3-rc7 (we have 2.4.1, but the CEGMA webpage says 2.4.1 is fine too), and geneid 2.4 (which we have). To run the test of CEGMA v2.4, I typed:
% export CEGMA=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312
% export PERL5LIB=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312/lib:$PERL5LIB
% export PERL5LIB=/software/pathogen/external/apps/usr/local/cegma_v2.4.010312/src:$PERL5LIB
% export CEGMATMP=$CEGMA_dir/temp
% $CEGMA/bin/cegma --tmp --ext -g /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.dna /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.prot
Note I submitted this to the farm3 with 2000 Mbyte of memory, otherwise sometimes it ran out of memory and blast failed as a result (with a segmentation fault).
My run gave slightly different results than expected according to the CEGMA webpage, I got:
Found 2079 candidate regions in /software/pathogen/external/apps/usr/local/cegma_v2.4.010312/sample/sample.dna
NOTE: created 87 geneid predictions
NOTE: Found 20 geneid predictions with scores above threshold value
DATA COLLECTED: 20 Coding sequences containing 57 introns
NOTE: created 46 geneid predictions
NOTE: Foud 20 geneid predictions with scores above threshold value
However, according to the CEGMA webpage, this is what is expected:
Found 86 candidate regions in /path/to/CEGMA/sample/sample.dna
NOTE: created 23 geneid predictions
NOTE: Found 15 geneid predictions with scores above threshold value
DATA COLLECTED: 15 Coding sequences containing 48 introns
NOTE: created 21 geneid predictions
NOTE: Foud 15 geneid predictions with scores above threshold value
My output.completeness.report file says:
#Prots %Completeness - #Total Average %Ortho
Complete 8 3.23 - 8 1.00 0.00
Partial 10 4.03 - 10 1.00 0.00
So it looks like a higher completeness was found using this test set than for CEGMA v2 on farm2 (see above).
I'm not sure why I get different results than on the CEGMA webpage. One difference is that I am using GeneWise 2.4.1 (not 2.2.3-rc7, the recommended one). However, a strange finding is that my output says that 2079 candidate regions were found using blast, while the CEGMA webpage says that 86 candidate regions were found, even though the blast versions are the same (blast+ 2.2.28).
I'm not sure about what difference using GeneWise 2.2.3-rc7 would make. I've tried to install this on farm3, but no luck so far compiling it..
Alternatives to CEGMA
An alternatives
Thanks
Thanks to Eleanor Stanley and Alan Tracey for help with CEGMA.
Tuesday, 19 March 2013
Making quantile-quantile plots (probability plots) in R
To check whether your sample of data is likely to have come from a population with a particular underlying probability distribution, it is useful to make a quantile-quantile plot (probability plot) for your data.
Making a Normal probability plot
To investigate whether a Normal distribution is a appropriate for modelling your data, you can make a Normal probability plot. To make a Normal probability plot of your data, you can use the function NormalProbPlot():
> NormalProbPlot <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
plot(x[oo], normvals, xlab="Data", ylab="yi", pch=20)
}
For example, we can use it as follows:
> x <- c(4, 0, -12, -18, 4, 12, -6, -16)
> NormalProbPlot(x)
For the model to be a good fit for the data, the points on a Normal probability plot should lie close to a line (it should in theory pass through the origin, but in practice can pass nearby). In this case, the data don't really lie on a straight line, so it is not very convincing that a Normal distribution is appropriate for modelling these data. However, the sample size is small, and the evidence against a Normal distribution isn't very strong either.
Let's try using a random sample of 5000 drawn from a Normal distribution:
> x <- rnorm(5000)
> NormalProbPlot(x)
We see a nice straight line.
Now let's try using a random sample of 5000 drawn from a continuous uniform distribution:
> x <- runif(5000)
> NormalProbPlot(x)
We see that the plot differs from a straight line. The pattern is characteristic of a distribution with tails that are too 'light' compared to a Normal distribution.
Let's try sampling 5000 points from an exponential distribution:
> x <- rexp(5000)
> NormalProbPlot(x)
This is clearly not a straight line. In fact, this curve is typical of what you see when you make a Normal probability plot for a very right-skewed data sample, like one originating from an exponential distribution.
Note that another way of making a Normal probability plot in R is to use the qqnorm() and qqline() functions:
> qqnorm(x)
> qqline(x)
Note that this plot shows the quantiles of the sample data on the y-axis and the quantiles of a theoretical Normal distribution on the x-axis, which is the opposite of the plot above, although it is the exact same data.
In fact, people often make their plot this way; you can also do it using this function:
> NormalProbPlot2 <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
sortedx <- x[oo]
plot(normvals, sortedx, xlab="yi", ylab="Data", pch=20)
}
> NormalProbPlot2(x)
A half-Normal plot
Another type of Normal plot is a 'half-Normal plot', which consists of the negative half of the Normal probability plot superimposed on the positive half:
> HalfNormalProbPlot <- function(x)
{
x <- c(abs(x),-abs(x))
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
sortedx <- x[oo]
plot(normvals[normvals>0], sortedx[normvals>0], xlab="yi", ylab="Data", pch=20)
}
> HalfNormalProbPlot(x)
Making an exponential probability plot
Similarly, to investigate whether an exponential distribution is appropriate for modelling your data, you can make an exponential probability plot. This can be done using the ExpProbPlot function:
> ExpProbPlot <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
expvals <- qexp(quantiles)
plot(x[oo], expvals, xlab="Data", ylab="yi", pch=20)
}
For example, we can use it as follows:
> x <- c(841, 158, 146, 45, 34, 122, 151, 281, 435, 737, 585, 888, 264, 1902,
696, 295, 563, 722, 77, 711, 47, 403, 195, 760, 320, 461, 41, 1337, 336, 1355,
455, 37, 668, 41, 557, 100, 305, 377, 568, 140, 781, 204, 437, 31, 385, 130, 10,
210, 600, 84, 833, 329, 247, 1618, 639, 938, 736, 39, 366, 93, 83, 221)
> ExpProbPlot(x)
For the model to be a good fit for the data, the points on an exponential probability plot should lie close to a line through the origin. In this case the data do lie approximately along a straight line through the origin, so it seems that an exponential distribution is a plausible model for the data.
Making a Normal probability plot
To investigate whether a Normal distribution is a appropriate for modelling your data, you can make a Normal probability plot. To make a Normal probability plot of your data, you can use the function NormalProbPlot():
> NormalProbPlot <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
plot(x[oo], normvals, xlab="Data", ylab="yi", pch=20)
}
For example, we can use it as follows:
> x <- c(4, 0, -12, -18, 4, 12, -6, -16)
> NormalProbPlot(x)
For the model to be a good fit for the data, the points on a Normal probability plot should lie close to a line (it should in theory pass through the origin, but in practice can pass nearby). In this case, the data don't really lie on a straight line, so it is not very convincing that a Normal distribution is appropriate for modelling these data. However, the sample size is small, and the evidence against a Normal distribution isn't very strong either.
Let's try using a random sample of 5000 drawn from a Normal distribution:
> x <- rnorm(5000)
> NormalProbPlot(x)
We see a nice straight line.
Now let's try using a random sample of 5000 drawn from a continuous uniform distribution:
> x <- runif(5000)
> NormalProbPlot(x)
We see that the plot differs from a straight line. The pattern is characteristic of a distribution with tails that are too 'light' compared to a Normal distribution.
Let's try sampling 5000 points from an exponential distribution:
> x <- rexp(5000)
> NormalProbPlot(x)
This is clearly not a straight line. In fact, this curve is typical of what you see when you make a Normal probability plot for a very right-skewed data sample, like one originating from an exponential distribution.
Note that another way of making a Normal probability plot in R is to use the qqnorm() and qqline() functions:
> qqnorm(x)
> qqline(x)
Note that this plot shows the quantiles of the sample data on the y-axis and the quantiles of a theoretical Normal distribution on the x-axis, which is the opposite of the plot above, although it is the exact same data.
In fact, people often make their plot this way; you can also do it using this function:
> NormalProbPlot2 <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
sortedx <- x[oo]
plot(normvals, sortedx, xlab="yi", ylab="Data", pch=20)
}
> NormalProbPlot2(x)
A half-Normal plot
Another type of Normal plot is a 'half-Normal plot', which consists of the negative half of the Normal probability plot superimposed on the positive half:
> HalfNormalProbPlot <- function(x)
{
x <- c(abs(x),-abs(x))
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
normvals <- qnorm(quantiles)
sortedx <- x[oo]
plot(normvals[normvals>0], sortedx[normvals>0], xlab="yi", ylab="Data", pch=20)
}
> HalfNormalProbPlot(x)
Making an exponential probability plot
Similarly, to investigate whether an exponential distribution is appropriate for modelling your data, you can make an exponential probability plot. This can be done using the ExpProbPlot function:
> ExpProbPlot <- function(x)
{
oo <- order(x)
length <- length(x)
quantiles <- seq(1/(1+length),1-(1/(1+length)),1/(1+length))
expvals <- qexp(quantiles)
plot(x[oo], expvals, xlab="Data", ylab="yi", pch=20)
}
For example, we can use it as follows:
> x <- c(841, 158, 146, 45, 34, 122, 151, 281, 435, 737, 585, 888, 264, 1902,
696, 295, 563, 722, 77, 711, 47, 403, 195, 760, 320, 461, 41, 1337, 336, 1355,
455, 37, 668, 41, 557, 100, 305, 377, 568, 140, 781, 204, 437, 31, 385, 130, 10,
210, 600, 84, 833, 329, 247, 1618, 639, 938, 736, 39, 366, 93, 83, 221)
> ExpProbPlot(x)
For the model to be a good fit for the data, the points on an exponential probability plot should lie close to a line through the origin. In this case the data do lie approximately along a straight line through the origin, so it seems that an exponential distribution is a plausible model for the data.
Subscribe to:
Posts (Atom)