Question: Accounting for hidden confounding factors in eQTL analysis
gravatar for Alexander Skates
3.7 years ago by
United Kingdom
Alexander Skates320 wrote:

I'm attempting to conduct an eQTL analysis. Broadly speaking, I'm attempting to use Matrix eQTL to perform permutation based tests, including inferred hidden PEER factors as covariates.

I've batch corrected the expression data using ComBat, and I've used PEER to generate the covariate matrix, however I can't use it with Matrix eQTL. PEER generates a matrix with dimensions [number of genes]x[number of factors], whereas Matrix eQTL is expecting one with dimensions [number of factors]x[number of samples]. How do I integrate the two packages?

rna-seq eqtl peer R matrix eqtl • 2.7k views
ADD COMMENTlink modified 3.5 years ago by informatics bot560 • written 3.7 years ago by Alexander Skates320
gravatar for informatics bot
3.5 years ago by
United States
informatics bot560 wrote:

PEER generates multiple output data frames after:


I'm assuming the matrix you're generating comes from:


If you want to combine PEER with MatrixEQTL, I would suggest using the "residual matrix" from:


As your "gene" variable within:


I've attach some code just to show you the basic work-flow:

model = PEER()
    PEER_setPhenoMean(model,as.matrix(expr_cov)) # Give PEER your expression data
    PEER_setAdd_mean(model, TRUE)
    PEER_setCovariates(model, as.matrix(covar4))  # PEER ask NxC matrix, where N=samples and C=covariates
    PEER_setNmax_iterations(model, 100000)
    residuals = PEER_getResiduals(model)  # convert to GxN
    colnames(residuals) = colnames(expr_cov)
    rownames(residuals) = rownames(expr_cov)
    # add mean
    residuals = residuals + apply(expr_cov, 2, mean)
    # convert to SliceData format
    genes = SlicedData$new();
    genes$CreateFromMatrix(residuals); # using your residuals as expression input
    rm(residuals, model)
    # run eQTL with residuals
    me = Matrix_eQTL_main(
        snps = snps,
        gene = genes,
        cvrt = SlicedData$new(),  # or cvrt_snp
        output_file_name = "",
        pvOutputThreshold = 0,  # no tran
        useModel = modelLINEAR,
        errorCovariance = numeric(),
        verbose = FALSE,
        output_file_name.cis = paste0("step1.cis.eQTL.K",k,".xls"),
        pvOutputThreshold.cis = 1e-5,
        snpspos = snp_loc,
        genepos = exons_8_25_2015,
        cisDist = 1e6,
        pvalue.hist = FALSE, = FALSE,
        noFDRsaveMemory = TRUE);
ADD COMMENTlink written 3.5 years ago by informatics bot560
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1199 users visited in the last hour