Question: How to recover treated/control count from DESeq2 output
1
gravatar for gundalav
2.3 years ago by
gundalav280
La La Land
gundalav280 wrote:

I have the following DESeq2 code and the corresponding results

suppressMessages(library(DESeq2))
suppressMessages(library(airway))

data(airway)
airway_se <- airway
airway_dds <- DESeqDataSet(airway_se, design = ~cell + dex)
deseq <- DESeq(airway_dds)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
results <- results(deseq)
results
#> log2 fold change (MAP): dex untrt vs trt 
#> Wald test p-value: dex untrt vs trt 
#> DataFrame with 64102 rows and 6 columns
#>                  baseMean log2FoldChange      lfcSE       stat
#>                 <numeric>      <numeric>  <numeric>  <numeric>
#> ENSG00000000003 708.60217     0.37415246 0.09884435  3.7852692
#> ENSG00000000005   0.00000             NA         NA         NA
#> ENSG00000000419 520.29790    -0.20206175 0.10974241 -1.8412367
#> ENSG00000000457 237.16304    -0.03616686 0.13834540 -0.2614244
#> ENSG00000000460  57.93263     0.08445399 0.24990709  0.3379415
#> ...                   ...            ...        ...        ...
#> LRG_94                  0             NA         NA         NA
#> LRG_96                  0             NA         NA         NA
#> LRG_97                  0             NA         NA         NA
#> LRG_98                  0             NA         NA         NA
#> LRG_99                  0             NA         NA         NA
#>                       pvalue        padj
#>                    <numeric>   <numeric>
#> ENSG00000000003 0.0001535423 0.001289269
#> ENSG00000000005           NA          NA
#> ENSG00000000419 0.0655868795 0.197066711
#> ENSG00000000457 0.7937652416 0.913856017
#> ENSG00000000460 0.7354072415 0.884141575
#> ...                      ...         ...
#> LRG_94                    NA          NA
#> LRG_96                    NA          NA
#> LRG_97                    NA          NA
#> LRG_98                    NA          NA
#> LRG_99                    NA          NA

My question is how can I recover the treated and control count for each gene to calculate the fold change in log2FoldChange output above?

dseq2 rna-seq • 1.8k views
ADD COMMENTlink modified 2.3 years ago by Benn7.1k • written 2.3 years ago by gundalav280
7
gravatar for Santosh Anand
2.3 years ago by
Santosh Anand4.9k
Santosh Anand4.9k wrote:
# Un-normalized counts
counts(deseq)

# Normalized counts  (Normalized for size factors)
counts(deseq, normalized = TRUE)

Be aware that the log2FoldChange reported by DESeq2 is shrunk. So it will be usually lesser than what calculated from your normalized counts data. See the DESeq2 paper, especially the Fig.1 https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8

That said, you can also get the MLE (or not shrunken) estimate of Log2FC by using addMLE = T in results()

results(deseq, addMLE = T)

This will add a column named lfcMLE in output, which should be closer to the log2FC calculated from normalized data.

ADD COMMENTlink modified 11 months ago • written 2.3 years ago by Santosh Anand4.9k
3
gravatar for Benn
2.3 years ago by
Benn7.1k
Netherlands
Benn7.1k wrote:

It's in the airway object, but it's a S4 object, so you'll have to put out the right slots.

my_counts <- airway@assays$data$counts

To add the colnames and rownames:

colnames(my_counts)<-airway@colData@listData$SampleName

rownames(my_counts) <- airway@rowData@partitioning@NAMES

head(my_counts)
ADD COMMENTlink written 2.3 years ago by Benn7.1k

Thanks. airway@assays$data$counts is that normalized or not? If not how can I get the normalized one from S4?

ADD REPLYlink written 2.2 years ago by gundalav280

This answer A: How to recover treated/control count from DESeq2 output tells you how to get the counts object (normalized) from a deseq/S4 object

ADD REPLYlink written 2.2 years ago by WouterDeCoster40k
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: 1685 users visited in the last hour