How to do DESeq2 for Kallisto
0
0
Entering edit mode
3.6 years ago
Kumar ▴ 170

Hi, I am trying to perform DESeq2 for Kallisto. I ran Kallisto for 36 PE samples and got results files such as abundance.h5. Now I am trying DESeq2 analysis. I am getting everything correct but while making plotPCA it is making the plot of 6 samples. However, I have a total of 36 samples, I tried to find out the issue but not sure why the script taking only 6 samples for the analysis.

Here is the script:

files <- c("/DataAnalysis/270-aligned/abundance.h5", 
 "/DataAnalysis/272-aligned/abundance.h5") ######just displaying 2 samples here, I have total 36 samples 

names(files) <- paste0("sample", 1:length(files))
txi.kallisto <- tximport(files, type = "kallisto", txOut = TRUE)
head(txi.kallisto$counts)

sampleTable <- data.frame(condition = factor(rep(c("mutant", "control"), each = 3)))
rownames(sampleTable) <- colnames(countdata$counts)

dds <- DESeqDataSetFromTximport(countdata, sampleTable, ~condition)
#PCA plot
vsd <- vst(dds)
class(vsd)
head(colData(vsd))
plotPCA(vsd, "condition")

output

> head(colData(vsd))
DataFrame with 6 rows and 1 column
    condition
     <factor>
sample1    mutant
sample2    mutant
sample3    mutant
sample4   control
sample5   control
sample6   control
Kallisto rna-seq DESeq2 • 2.3k views
ADD COMMENT
0
Entering edit mode
sampleTable <- data.frame(condition = factor(rep(c("mutant", "control"), each = 3)))

You only specify 6 samples. What is CountData? Are all samples imported by tximport?

ADD REPLY
0
Entering edit mode

I believe so. All data has been imported by tximport. I changed CountData by txi.kallisto variable but it is showing following error. I am following the process of analysis at https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html. Where do I need to specify samples?

> dds <- DESeqDataSetFromTximport(txi.kallisto, sampleTable, ~condition)
Error in DESeqDataSetFromMatrix(countData = counts, colData = colData,  : 
ncol(countData) == nrow(colData) is not TRUE
ADD REPLY
0
Entering edit mode

Well, you are hand making your ColData there with only 6 samples. If you have count values for 36 samples, how can they possibly be the same? Ar you just copy-pasting that code from some tutorial?

ADD REPLY
0
Entering edit mode

I am not catching entirely. I am following the tutorial I am new in this analysis. I tried to change but not getting improved. If you could help me to improve the code.

sampleTable <- data.frame(condition = 
factor(rep(c("mutant","mutant","mutant","mutant","mutant","mutant","mutant","mutant","mutant",
                                             "mutant","mutant","mutant","mutant","mutant","mutant","mutant","mutant","mutant",
                                               "control","control","control","control","control","control","control","control",
                                               "control","control","control","control","control","control","control","control",
                                              "control","control"), each = 36)))

rownames(sampleTable) <- colnames(txi.kallisto$counts)

dds <- DESeqDataSetFromTximport(txi.kallisto, sampleTable, ~condition)

error:

Error in `.rowNamesDF<-`(x, value = value) : invalid 'row.names' length
ADD REPLY
0
Entering edit mode

Remove the rep(c(....),each=36) part. If you explicitely write down each factor level then you don't have to rep it.

ADD REPLY
0
Entering edit mode

I got it! if I write with rep as following:

sampleTable <- data.frame(condition = factor(rep(c("mutant", "control"), each = 18)))
ADD REPLY
0
Entering edit mode

However, I need one more help in the following code. What should I use here instead of ColData.

vsd <- vst(dds)
class(vsd)
head(ColData(vsd))
plotPCA(vsd, "condition")
ADD REPLY
0
Entering edit mode

Want the values? assay(vsd), read the manual!

ADD REPLY
0
Entering edit mode

SampleTable is 1296 lines long. Are you so confident in your coding skills that you didn't even look at it before trying to use it?

ADD REPLY
0
Entering edit mode

I did not know about it. If you could really want to help, you could have sent the link instead of commenting on me. I mentioned that I am new about the analysis and do not know about the coding either.

ADD REPLY

Login before adding your answer.

Traffic: 2004 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