Question: svaseq and combat
gravatar for Folder40g
3.7 years ago by
Folder40g140 wrote:


I have samples (rnaseq) coming from two different experiments. I plotted all samples in a PCA plot and I saw a clear batch effect. Then I used svaseq to plot a PCA without the batch effect and I got something quite reasonable, samples from different experiments mix but without a clear order. Samples from each experiment in the PCA cluster as expected.

But then I run combat, write the output and createa a heatmap with it in morpheus and the clustering that I get is the one I would get without batch effect correction.

Why is this so? What i'm doing wrong?



    count <- read.table("/local/deseq/ALLsamples.txt", header = T, row.names = 1)
    count <- log(count + 1)

    #-- Filtering counts min >=2 genes having >5 counts
    filter <- apply(count, 1, function(x) length(x[x>5])>=2)
    filtered <- count[filter,]
    dat0 <- as.matrix(filtered)

-- Phenotype

    pheno <-"Daf", 10), rep("Tuv",19)), c(rep("Dox",5), rep("Doxless",5), rep("mN",7), rep("mP",6), rep("mT",6))))
    rownames(pheno) <- colnames(count)
    colnames(pheno) <- c("Experiment", "GroupSample")

-- Models

    #-- Full model = adjustment variables + variables of interest 
    #mod = model.matrix(~as.factor(Experiment) + as.factor(GroupSample), data=pheno)
    mod <- model.matrix(~as.factor(Experiment), data=pheno)

    #-- Null model = adjustment variables 
    mod0 <- model.matrix(~1, data=pheno)

    #-- Number of surrogated variables n-sv,mod,method="leek") #-- Value is 1

-- svaseq + PCA plot

    #-- svaseq
    svseq <- svaseq(as.matrix(filtered), mod = mod, mod0 = mod0, = 1)
    svseq_plot <- svaseq(as.matrix(filtered), mod = mod, mod0 = mod0, = 1)$sv

    #-- Plot
    plot(svseq_plot, pch=19, col="blue")

    #-- Clean original input and obtain corrected counts
    cleaningY = function(y, mod, svaobj) {
      X = cbind(mod, svaobj$sv)
      Hat = solve(t(X)%*%X)%*%t(X)
      beta = (Hat%*%t(y))
      P = ncol(mod)
      cleany = y-t(as.matrix(X[,-c(1:P)])%*%beta[-c(1:P),])

    #-- Clean
    counts_sva = cleaningY(count, mod, svseq)
    pca = prcomp(t( counts_sva ), scale=FALSE)

    plot(pca$x[,1], pca$x[,2])
    text(pca$x[,1], pca$x[,2], rownames(pca$x), pos= 3 )

-- Combat

    mod_batch <- mod[,2]
    combat_adj <- ComBat(as.matrix(filtered), batch = mod_batch, mod = mod0, par.prior=TRUE, prior.plots=FALSE)

    write.table(combat_adj, "/local/CombatAdjusted.txt", sep ="\t")
batch effect sva combat • 3.0k views
ADD COMMENTlink modified 22 months ago by liull_hust0 • written 3.7 years ago by Folder40g140

Hi,in the process of the PCA plot,why you choose scale=FALSE ?I have tried my data with the scale TRUE or FALSE ,the results are very different.

ADD REPLYlink written 22 months ago by liull_hust0

Scaling or not will depend on the distribution of your input data.

ADD REPLYlink written 22 months ago by Kevin Blighe69k
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: 873 users visited in the last hour