Question: MA plot from my data frame
0
gravatar for bgraphit
3.8 years ago by
bgraphit20
United States
bgraphit20 wrote:

Hello! I am trying to figure out how I can make an MA plot from the following my data frame "dar"

I am currently using DESeq2 but I extracted the information from data frame "dr" which is the results (dds) I saved this as a csv.

Opened another R session and loaded "dr" where I can have the column with gene name. I changed the header of that column to "gene"

then merged "dr" by "gene" it to my list of genes of interest to make my data frame dim.

So dim at the moment isnt seen as a DESeq2 object data frame.

dim(dar)

[1] 13348     6

 

head(dar)
                           baseMean log2FoldChange     lfcSE          stat
A1BG                   7.278775e-02    0.140265235 1.6670434   0.084140124
A1CF                   4.130988e-02   -1.102589688 1.6676373  -0.661168755
AAAS                   1.657182e+01    1.173655704 0.5986687   1.960442661
AACS                   6.373922e+00   -0.653331952 0.8046976  -0.811897466
AADAC                  3.251195e+00    0.736454731 0.7014096   1.049963935

                            pvalue                    padj
A1BG                   9.329450e-01           NA
A1CF                   5.085041e-01           NA
AAAS                   4.994407e-02 1.377416e-01
AACS                   4.168505e-01 6.001605e-01
AADAC                  2.937347e-01 4.793114e-01
 

 

I want to make an MA plot with the plotMA function where x axis is mean expression("baseMean") and y axis is log2FoldChange

plotMA(Merged_H3K27me3_not_inforMAplot_nocols, ylim=c(-5,5))

 

How can I get this data frame back into DESeq2 format? and make my MA plot?

 

plotma R • 2.2k views
ADD COMMENTlink modified 3.8 years ago • written 3.8 years ago by bgraphit20
1
gravatar for Steven Lakin
3.8 years ago by
Steven Lakin1.4k
Fort Collins, CO, USA
Steven Lakin1.4k wrote:

Based on the way your dataframe is formatted, I'm guessing you're using DESeq2.  Make sure DESeq2 was the last package you loaded that has a plotMA function; there are many other packages with a function of the same name, so whichever you loaded last would be the one you're using.  Can you confirm that you're using the DESeq2 plotMA function and report back if you have further problems?

Try the following:

library(DESeq2)
plotMA(dar, main="DESeq2", ylim=c(-2,2))

 

Or:

DESeq2::plotMA(dar, main="DESeq2", ylim=c(-2,2))
ADD COMMENTlink written 3.8 years ago by Steven Lakin1.4k

Hi Steven

I am using DESeq2 but I extracted the information from "dr" which is the results (dds) I saved this as a csv.

Opened another R session and loaded "dr" where I can have the column with gene name. I changed the header of that column to "gene"

then merged "dr" by "gene" it to my list of genes of interest to make my data frame dim.

So dim at the moment isnt seen as a DESeq2 object data frame.

 

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by bgraphit20
1

To use the DESeq plotMA, you need to maintain the original architecture of the dds result object.  Try the following in the same R environment as your processed dds object:

res <- results(dds)
myGenesOfInterest <- c("Gene1", "Gene2")   # place all genes of interest in this vector
dar <- res[rownames(res) %in% myGenesOfInterest, ]
plotMA(dar, main="DESeq2", ylim=c(-2,2))

Keep in mind your gene names should be your row names, and row names is just another column you can manipulate using rownames().

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by Steven Lakin1.4k

My genesofInterest list looks like this:

> Merged_Genes_notin_H3K27
                        gene                    samplex     ......

2                       A1BG                      0.000000
3                       A1CF                      0.000000
4                       AAAS                      0.000000
5                       AACS                      4.247928
6                      AADAC                      0.000000
7                AADACL2-AS1                      0.000000
8                    AADACP1                      6.554589
9                      AADAT                      2.584963

Is there a way to use this instead of the vector your describe?
 

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by bgraphit20
1

Yep, assuming row 1 is also in the data and "gene" is the column name, this should work:

res <- results(dds)
dar <- res[rownames(res) %in% Merged_Genes_notin_H3K27$gene, ]
plotMA(dar, main="DESeq2", ylim=c(-2,2))

 

If row 1 is not needed, then this should work:

res <- results(dds)
genesOfInterest <- as.character(Merged_Genes_notin_H3K27[2:nrow(Merged_Genes_notin_H3K27), 1])

dar <- res[rownames(res) %in% genesOfInterest, ]
plotMA(dar, main="DESeq2", ylim=c(-2,2))
ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by Steven Lakin1.4k

ok so i did

res2 <- results(dds)
> dim(res2)
[1] 25369     7
dar <- res2[rownames(res2) %in% Merged_Genes_notin_H3K9me3$gene, ]

> names(Merged_Genes_in_H3K9me3)
  [1] "gene" 

and  the problem is that dar 

> dar
[1] gene           baseMean       log2FoldChange lfcSE          stat          
[6] pvalue         padj          
<0 rows> (or 0-length row.names)

plotMA(dar, main="DESeq2", ylim=c(-2,2))

 

This is odd because there are genes that match in both lists.

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by bgraphit20

Try recopying the second block of code, where I explicitly converted it to character.  Sometimes these column types can be tricky and convert to factors.

ADD REPLYlink written 3.8 years ago by Steven Lakin1.4k

>genesOfInterest <- as.character(Merged_Genes_in_H3K9me3[1:nrow(Merged_Genes_in_H3K9me3), 1])

> genesOfInterest
  [1] "ABCA11P"       "ABCA13"        "ABI3BP"        "ADCY2"        
  [5] "AMN1"          "ANKMY1"        "ANKRD20A12P"   "ANKRD20A9P"   
  [9] "ANKRD30BP2"    "ANKRD36BP1"    "APMAP"         "ARHGAP39"     
 [1

> dar <- res2[rownames(res2) %in% genesOfInterest, ]

> dar
[1] gene           baseMean       log2FoldChange lfcSE          stat          
[6] pvalue         padj          
<0 rows> (or 0-length row.names)

> dim(dar)
[1] 0 7

no still not working well

 

ADD REPLYlink written 3.8 years ago by bgraphit20

That is strange.  It worked for me on a test dataset from DESeq2 vignette.  Did you make sure to recall the results(dds) from your analyzed dds object?

What is the output of:

str(res2)
ADD REPLYlink written 3.8 years ago by Steven Lakin1.4k

Actually yes it worked. It was my mistake that when I loaded the rdata file that had my list of genes of interest it had a res2 object derived from a .csv file that overwrote over the res2 I needed.

Thanks!!

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by bgraphit20
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: 1094 users visited in the last hour