Question: svaseq and combat
0
gravatar for Folder40g
2.3 years ago by
Folder40g120
Folder40g120 wrote:

Hi

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?

Thanks

    library(sva)
    library(limma)

    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 <- as.data.frame(cbind(c(rep("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
    num.sv(as.matrix(filtered),mod,method="leek") #-- Value is 1

-- svaseq + PCA plot

    #-- svaseq
    svseq <- svaseq(as.matrix(filtered), mod = mod, mod0 = mod0, n.sv = 1)
    svseq_plot <- svaseq(as.matrix(filtered), mod = mod, mod0 = mod0, n.sv = 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),])
      return(cleany)
    }

    #-- 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 • 1.8k views
ADD COMMENTlink modified 6 months ago by liull_hust0 • written 2.3 years ago by Folder40g120

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 6 months ago by liull_hust0

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

ADD REPLYlink written 6 months ago by Kevin Blighe48k
Please log in to add an answer.

Help
Access

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