3

I am a complete newbie to MATLAB and the first task I have is to create a volcano plot. I have been using this link to understand about it and get started.

I tried to run it on dummy values - 

a=[1 2 3]

b=[4.6 2.7 4.5]

c=[0.05 0.33 0.45]

And then I ran - 

SigStructure = mavolcanoplot(a, b, c)

My understanding is that a represents the gene expression values for condition 1, b for condition 2, and c is the list of p-values for the 3 data points in a and b.

However running this code gives me the error - 

Index exceeds matrix dimensions.

Error in mavolcanoplot (line 127)
appdata.effect = X(paramStruct.goodVals) - Y(paramStruct.goodVals);

Error in volc (line 4)
SigStructure = mavolcanoplot(a, b, c)

Can anyone explain where I am going wrong?

llrs
  • 4,693
  • 1
  • 18
  • 42
user1995
  • 173
  • 2
  • 8

3 Answers3

3

I'm not sure if nowadays matlab is still a good option for informatics. I would rather go for a R script to make a volcano plot. If you don't want to have fancy code but functional one, you can even make the foldchange calculations and log2 transformation in excel and just plot the log2 fold change in x axis and the -log10 transform of the p/q values on the y axis. If you need it, I could provide one of my own scripts.

enter image description here

Update: here is the actual script:

#load data into the 'res' variable
res <- read.table("genes.txt", header=TRUE)
head(res)

# Make a basic volcano plot with grey dots
with(res, plot(log2FoldChange, -log10(qvalue), pch=19, col="#4A4E4D", xlab="Log2 Fold Change", xlim=c(-8,8)))

# Add colored points: "indianred" if padj<0.0001, orange of log2FC>2, green if both)
with(subset(res, qvalue<.0001 ), points(log2FoldChange, -log10(qvalue), pch=19, col="#17B090"))
with(subset(res, abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#F89538"))
with(subset(res, qvalue<.0001 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#FA6367"))
llrs
  • 4,693
  • 1
  • 18
  • 42
Erik
  • 56
  • 3
  • 1
    hey Erik, it would be really helpful for me if you could provide me the scripts you mentioned. I am a real noob in this, and am trying to figure out the best way to go about it – user1995 Oct 04 '17 at 20:19
2

so while I am not so familiar with MATLAB, in the definition of the volcano plot it specifies that PValues is supposed to be a column vector or a DataMatrix object. Your c vector is a row vector, and this confuses MATLAB.

galicae
  • 349
  • 1
  • 5
2

Is there a reason why you are doing this in MATLAB instead of R?

In R with the ggplot2 package you can make a volcano plot pretty easily using something like this (may be errors as it is off the top of my head and untested):

# load libraries
library(ggplot2)
library(data.table)

# create a dataframe/data.table of the data
a = c(1,2,3)
b = c(4.6, 2.7, 4.5)
c = c(0.05, 0.33, 0.45)
dt = data.table(A=a,B=b,pval=c)

# add column to calculate log2 fold change between a/b
dt[,log2Diff := (log2(A)-log2(B))]
# add column to transform pvalue to -log scale
dt[,logPval := -log10(pval)]

#make the plot
ggplot(dt, aes(x=log2Diff,y=logPval))+geom_point()
Reilstein
  • 367
  • 1
  • 14