Skip to content

Commit

Permalink
remove unused functions in snpgdsEIGMIX
Browse files Browse the repository at this point in the history
some functions in the last version of snpgdsEIGMIX() have been moved to
snpgdsGRM()
  • Loading branch information
zhengxwen committed Oct 5, 2015
1 parent 70c6763 commit 9e1ad00
Show file tree
Hide file tree
Showing 6 changed files with 24 additions and 52 deletions.
3 changes: 2 additions & 1 deletion NEWS
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,8 @@ CHANGES IN VERSION 1.3.0-1.3.11

o `snpgdsLDMat()` supports multiple threads and covariance

o `snpgdsPCA()`: non-computed eigenvalues are NaN to avoid misuse of 'eigenval' when `eigen.method="DSPEVX"`
o `snpgdsPCA()`: non-computed eigenvalues are NaN to avoid misuse of
'eigenval' when `eigen.method="DSPEVX"`


CHANGES IN VERSION 1.2.0
Expand Down
33 changes: 9 additions & 24 deletions R/PCA.R
Original file line number Diff line number Diff line change
Expand Up @@ -195,8 +195,8 @@ snpgdsPCASampLoading <- function(loadobj, gdsobj, sample.id=NULL,

snpgdsEIGMIX <- function(gdsobj, sample.id=NULL, snp.id=NULL,
autosome.only=TRUE, remove.monosnp=TRUE, maf=NaN, missing.rate=NaN,
num.thread=1L, eigen.cnt=32, need.ibdmat=FALSE, ibdmat.only=FALSE,
method=c("Coancestry", "GRM"), verbose=TRUE)
num.thread=1L, eigen.cnt=32L, need.ibdmat=FALSE, ibdmat.only=FALSE,
verbose=TRUE)
{
# check and initialize ...
ws <- .InitFile2(cmd="Eigen-analysis on SNP genotypes:",
Expand All @@ -205,31 +205,16 @@ snpgdsEIGMIX <- function(gdsobj, sample.id=NULL, snp.id=NULL,
maf=maf, missing.rate=missing.rate, num.thread=num.thread,
verbose=verbose)

stopifnot(is.numeric(eigen.cnt) & is.vector(eigen.cnt))
stopifnot(length(eigen.cnt) == 1)
eigen.cnt <- as.integer(eigen.cnt)
if (eigen.cnt < 1)
eigen.cnt <- as.integer(ws$n.samp)
stopifnot(is.numeric(eigen.cnt), length(eigen.cnt)==1L)
if (eigen.cnt < 1L)
eigen.cnt <- ws$n.samp

stopifnot(is.logical(need.ibdmat) & is.vector(need.ibdmat))
stopifnot(length(need.ibdmat) == 1)
stopifnot(is.logical(need.ibdmat), length(need.ibdmat)==1L)
stopifnot(is.logical(ibdmat.only), length(ibdmat.only)==1L)

stopifnot(is.logical(ibdmat.only) & is.vector(ibdmat.only))
stopifnot(length(ibdmat.only) == 1)

method <- match.arg(method)
diag.adjustment <- (method == "Coancestry")

in.method <- c("RatioOfAverages", "AverageOfRatios")
in.method <- as.character(in.method)[1]
in.method <- match(in.method, c("RatioOfAverages", "AverageOfRatios"))


#########################################################
# call eigen-analysis

rv <- .Call(gnrEIGMIX, eigen.cnt, ws$num.thread, need.ibdmat,
ibdmat.only, in.method, diag.adjustment, verbose)
rv <- .Call(gnrEIGMIX, eigen.cnt, ws$num.thread, need.ibdmat, ibdmat.only,
verbose)

# return
rv <- list(sample.id = ws$sample.id, snp.id = ws$snp.id,
Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ Development Version:
[http://www.bioconductor.org/packages/devel/bioc/html/SNPRelate.html](http://www.bioconductor.org/packages/devel/bioc/html/SNPRelate.html)


## News (Development Version)
## News

* Supports the [SeqArray](http://bioconductor.org/packages/release/bioc/html/SeqArray.html) GDS format, see [the vignette](http://www.bioconductor.org/packages/devel/bioc/vignettes/SeqArray/inst/doc/AnalysisTutorial.html).

Expand Down
8 changes: 3 additions & 5 deletions man/snpgdsEIGMIX.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,7 @@
\usage{
snpgdsEIGMIX(gdsobj, sample.id=NULL, snp.id=NULL, autosome.only=TRUE,
remove.monosnp=TRUE, maf=NaN, missing.rate=NaN, num.thread=1L,
eigen.cnt=32, need.ibdmat=FALSE, ibdmat.only=FALSE,
method=c("Coancestry", "GRM"), verbose=TRUE)
eigen.cnt=32L, need.ibdmat=FALSE, ibdmat.only=FALSE, verbose=TRUE)
}
\arguments{
\item{gdsobj}{an object of class \code{\link{SNPGDSFileClass}},
Expand All @@ -33,8 +32,6 @@ snpgdsEIGMIX(gdsobj, sample.id=NULL, snp.id=NULL, autosome.only=TRUE,
\item{need.ibdmat}{if TRUE, return the IBD matrix}
\item{ibdmat.only}{return the IBD matrix only, do not compute the
eigenvalues and eigenvectors}
\item{method}{"Coancestry" -- EIGMIX coancestry matrix (by default);
"GRM" -- genetic relationship matrix}
\item{verbose}{if TRUE, show information}
}
\value{
Expand Down Expand Up @@ -105,7 +102,8 @@ legend("topright", legend=levels(tab$pop), pch="o", col=1:4)


# run eigen-analysis
RV <- snpgdsEIGMIX(genofile, sample.id=samp.id[pop_code=="JPT"])
RV <- snpgdsEIGMIX(genofile, sample.id=samp.id[pop_code=="JPT"],
need.ibdmat=TRUE)
z <- RV$ibdmat

mean(c(z))
Expand Down
4 changes: 2 additions & 2 deletions src/SNPRelate.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1136,7 +1136,7 @@ COREARRAY_DLL_EXPORT void R_init_SNPRelate(DllInfo *info)
extern SEXP gnrConvBEDFlag(SEXP, SEXP, SEXP);
extern SEXP gnrConvBED2GDS(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
extern SEXP gnrDiss(SEXP, SEXP);
extern SEXP gnrEIGMIX(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
extern SEXP gnrEIGMIX(SEXP, SEXP, SEXP, SEXP, SEXP);
extern SEXP gnrFst(SEXP, SEXP, SEXP);
extern SEXP gnrHWE();
extern SEXP gnrGRM(SEXP, SEXP, SEXP);
Expand Down Expand Up @@ -1189,7 +1189,7 @@ COREARRAY_DLL_EXPORT void R_init_SNPRelate(DllInfo *info)

CALL(gnrDiss, 2),
// CALL(gnrDistPerm, ),
CALL(gnrEIGMIX, 7), CALL(gnrErrMsg, 0),
CALL(gnrEIGMIX, 5), CALL(gnrErrMsg, 0),
CALL(gnrFst, 3), CALL(gnrHWE, 0),

CALL(gnrGRM, 3),
Expand Down
26 changes: 7 additions & 19 deletions src/genPCA.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1679,8 +1679,7 @@ static bool _EigComp(const TEigPair &i, const TEigPair &j)

/// to compute the eigenvalues and eigenvectors
COREARRAY_DLL_EXPORT SEXP gnrEIGMIX(SEXP _EigenCnt, SEXP _NumThread,
SEXP _NeedIBDMat, SEXP _IBDMatOnly, SEXP _Method, SEXP _DiagAdj,
SEXP _Verbose)
SEXP _NeedIBDMat, SEXP _IBDMatOnly, SEXP _Verbose)
{
const bool verbose = SEXP_Verbose(_Verbose);

Expand All @@ -1700,27 +1699,16 @@ COREARRAY_DLL_EXPORT SEXP gnrEIGMIX(SEXP _EigenCnt, SEXP _NumThread,
// the upper-triangle IBD matrix
CdMatTri<double> IBD(n);

// Calculate the genetic covariace
switch (INTEGER(_Method)[0])
{
case 1:
PCA::DoAdmixCalc_RatioOfAvg(IBD, LOGICAL(_DiagAdj)[0]==TRUE,
INTEGER(_NumThread)[0], verbose);
break;
case 2:
PCA::DoAdmixCalc_AvgOfRatios(IBD, LOGICAL(_DiagAdj)[0]==TRUE,
INTEGER(_NumThread)[0], verbose);
break;
default:
throw "Invalid method.";
}
// calculate the EIGMIX coancestry matrix
PCA::DoAdmixCalc_RatioOfAvg(IBD, true, Rf_asInteger(_NumThread),
verbose);

// ======== The calculation of eigenvectors and eigenvalues ========

int nProtected = 0;
SEXP EigenVal=NULL, EigenVec=NULL, IBDMat=NULL;

if (LOGICAL(_NeedIBDMat)[0] == TRUE)
if (Rf_asLogical(_NeedIBDMat) == TRUE)
{
PROTECT(IBDMat = Rf_allocMatrix(REALSXP, n, n));
nProtected ++;
Expand All @@ -1737,7 +1725,7 @@ COREARRAY_DLL_EXPORT SEXP gnrEIGMIX(SEXP _EigenCnt, SEXP _NumThread,
}
}

if (LOGICAL(_IBDMatOnly)[0] != TRUE)
if (Rf_asLogical(_IBDMatOnly) != TRUE)
{
const size_t NN = n;
vector<double> tmp_Work(NN*3);
Expand All @@ -1750,7 +1738,7 @@ COREARRAY_DLL_EXPORT SEXP gnrEIGMIX(SEXP _EigenCnt, SEXP _NumThread,
NowDateToStr().c_str());
}

int eigencnt = INTEGER(_EigenCnt)[0];
int eigencnt = Rf_asInteger(_EigenCnt);
if (eigencnt > n) eigencnt = n;

PROTECT(EigenVal = NEW_NUMERIC(n));
Expand Down

0 comments on commit 9e1ad00

Please sign in to comment.