Accounting for hidden confounding factors in eQTL analysis
1
0
Entering edit mode
9.3 years ago

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?

eQTL RNA-Seq peer Matrix eQTL R • 5.0k views
ADD COMMENT
1
Entering edit mode
9.2 years ago

PEER generates multiple output data frames after:

PEER_update(model)

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

PEER_getX(model) 

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

PEER_getResiduals(model)

As your "gene" variable within:

Matrix_eQTL_main()

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

model = PEER()
    PEER_setNk(model,k)
    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)
    PEER_update(model)
    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,
        min.pv.by.genesnp = FALSE,
        noFDRsaveMemory = TRUE);
ADD COMMENT

Login before adding your answer.

Traffic: 1665 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6