modelLINEAR_CROSS: Constant for 'Matrix_eQTL_engine'.

In MatrixEQTL: Matrix eQTL: Ultra Fast eQTL Analysis via Large Matrix Operations

Description

Set parameter useModel = modelLINEAR_CROSS in the call of Matrix_eQTL_main to indicate that Matrix eQTL should include the interaction of SNP and the last covariate in the model and test for its significance.

References

The package website: http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/

See Also

See Matrix_eQTL_engine for reference and sample code.

Examples

library('MatrixEQTL')    

# Number of columns (samples)
n = 25;

# Number of covariates
nc = 10;

# Generate the standard deviation of the noise
noise.std = 0.1 + rnorm(n)^2;

# Generate the covariates
cvrt.mat = 2 + matrix(rnorm(n*nc), ncol = nc);

# Generate the vectors with single genotype and expression variables
snps.mat = cvrt.mat %*% rnorm(nc) + rnorm(n);
gene.mat = cvrt.mat %*% rnorm(nc) + rnorm(n) * noise.std + 
            1 + 0.5 * snps.mat + snps.mat * cvrt.mat[,nc];

# Create 3 SlicedData objects for the analysis
snps1 = SlicedData$new( matrix( snps.mat, nrow = 1 ) );
gene1 = SlicedData$new( matrix( gene.mat, nrow = 1 ) );
cvrt1 = SlicedData$new( t(cvrt.mat) );

# name of temporary output file
filename = tempfile();

# Call the main analysis function
me = Matrix_eQTL_main(
    snps = snps1, 
    gene = gene1, 
    cvrt = cvrt1, 
    output_file_name = filename, 
    pvOutputThreshold = 1, 
    useModel = modelLINEAR_CROSS, 
    errorCovariance = diag(noise.std^2), 
    verbose = TRUE,
    pvalue.hist = FALSE );
# remove the output file
unlink( filename );

# Pull Matrix eQTL results - t-statistic and p-value
beta = me$all$eqtls$beta;
tstat = me$all$eqtls$statistic;
pvalue = me$all$eqtls$pvalue;
rez = c(beta = beta, tstat = tstat, pvalue = pvalue)
# And compare to those from the linear regression in R
{
    cat('\n\n Matrix eQTL: \n'); 
    print(rez);
    cat('\n R summary(lm()) output: \n')
    lmodel = lm( gene.mat ~ snps.mat + cvrt.mat + snps.mat*cvrt.mat[,nc], 
                    weights = 1/noise.std^2);
    lmout = tail(summary(lmodel)$coefficients,1)[, c("Estimate", "t value", "Pr(>|t|)")];
    print( tail(lmout) );
}

# Results from Matrix eQTL and 'lm' must agree
stopifnot(all.equal(lmout, rez, check.attributes = FALSE));

Example output

Processing the error covariance matrix 
Task finished in  0.006  seconds
Processing covariates 
Task finished in  0.002  seconds
Processing gene expression data (imputation, residualization, etc.) 
Task finished in  0.001  seconds
Creating output file(s) 
Task finished in  0.005  seconds
Performing eQTL analysis 
100.00% done, 1 eQTLs
Task finished in  0.011  seconds
 


 Matrix eQTL: 
        beta        tstat       pvalue 
1.068443e+00 1.439134e+01 6.235456e-09 

 R summary(lm()) output: 
    Estimate      t value     Pr(>|t|) 
1.068443e+00 1.439134e+01 6.235456e-09 

MatrixEQTL documentation built on Dec. 22, 2019, 5:06 p.m.