Plotting Microarray Time Course Data

Question

[robjohn7000]

Hi,

I need to plot microarray time course data so that the gene expression profiles of the top 100 genes are shown on the graph, with the curves labeled with the replicate probes (Information below).

I was able to draw the graph frames with the right titles and axes labels, but label on the curves did not have the right names, using the following code:

  for(j in 1:12){
  pch_value = as.character(targets$grp[6*j])
  points(c(0, 6, 24, 48, 72,80), exprs.row[(6*j-1):(6*j)],
  type='b', pch=pch_value)

Contrasts:

mc = makeContrasts("hr6-hr0","hr24-hr0","hr48-hr0","hr72-hr0","hr80-hr0","hr72-hr6", levels = X)



colnames(X) = c('hr0', 'hr24', 'hr48', 'hr6', 'hr72', 'hr80')

Replicate information:

Replicates (hr0 - 6 reps, hr 6 - 6reps, hr24-4reps, hr48-4reps,hr72-6reps, hr80-2reps) are in the following order :
hr0
hr0
hr6
hr6
hr24
hr24
hr48
hr48
hr72
hr72
hr0
hr0
hr6
hr6
hr24
hr24
hr48
hr48
hr72
hr72
hr80
hr80
hr0
hr0
hr6
hr6
hr72
hr72

I have tried maSigPro and timecourse; I was unable to set the correct contrasts using these two programs.

Help will be appreciated.

[Michael]

I am not totally sure what you mean with "label on the curves". The code does not handle plotting "curve labels" explicitly. You are plotting an affy expression set, correct? I recommend you post an image of how your plot looks, and mark in the plot where you want labels to appear. It might be best to control your plot options explicitly, see for example documentation of legend and also axis, par, and text.

[Ying W]

Could you include output for X and mc. If you are having a labeling problem make sure to check your colnames after every step and especially before reassigning them

[robjohn7000]

Thanks everyone. I have added more information as requested. My figure example is here: http://s23.postimg.org/bqcp67jaj/Rplot01_page_1.jpg

It's as if I'm plotting the wrong things. I just want to show how the other time profiles are different from hr0 for each ranked gene and perhaps, I do not really need to label the profile graphs and should just get rid of the 'h' that appears as labels. Any suggestions will be welcome.

mc and X are as follows:

mc

Contrasts
Levels hr6-hr0 hr24-hr0 hr48-hr0 hr72-hr0 hr80-hr0 hr72-hr6
  hr0       -1       -1       -1       -1       -1        0
  hr24       0        1        0        0        0        0
  hr48       0        0        1        0        0        0
  hr6        1        0        0        0        0       -1
  hr72       0        0        0        1        0        1
  hr80       0        0        0        0        1        0

> X
   hr0 hr24 hr48 hr6 hr72 hr80
1    1    0    0   0    0    0
2    1    0    0   0    0    0
3    0    0    0   1    0    0
4    0    0    0   1    0    0
5    0    1    0   0    0    0
6    0    1    0   0    0    0
7    0    0    1   0    0    0
8    0    0    1   0    0    0
9    0    0    0   0    1    0
10   0    0    0   0    1    0
11   1    0    0   0    0    0
12   1    0    0   0    0    0
13   0    0    0   1    0    0
14   0    0    0   1    0    0
15   0    1    0   0    0    0
16   0    1    0   0    0    0
17   0    0    1   0    0    0
18   0    0    1   0    0    0
19   0    0    0   0    1    0
20   0    0    0   0    1    0
21   0    0    0   0    0    1
22   0    0    0   0    0    1
23   1    0    0   0    0    0
24   1    0    0   0    0    0
25   0    0    0   1    0    0
26   0    0    0   1    0    0
27   0    0    0   0    1    0
28   0    0    0   0    1    0
 attr(,"assign")
 [1] 1 1 1 1 1 1
 attr(,"contrasts")
 attr(,"contrasts")$expt_structure
 [1] "contr.treatment"

[Ying W]

Check that your pch_value is changing after every iteration. Also an alternative way could be to create a vector of values ahead of time and then do pch=pch_vector[j]

Answer 1

[Michael]

Your plot character is always "h", the first character of all your replicate names/contrasts. pch can only be a single character, try pch=1:6 instead, if you don't need those. Using different colors (e.g. col=1:6 as argument to plot) or lintypes (lty=1:6) should help differentiate them. Use legend() to add a legend explaining the lines if required.

[robjohn7000]

I appreciate the new suggestions. I was able to remove the 'h', but have 5 plots and I now think it's probably better for me to just have one plot as explained with new information and idea below:

I have reduced the contrasts (removed 80hr-6h) to:

          contrasts <- makeContrasts("hr6-hr0","hr24-hr0","hr48-hr0","hr72-hr0","hr80-hr0",levels = exptDesign)

Hence, the new contrasts:

contrasts

Contrasts
                Levels                  hr6-hr0     hr24-hr0        hr48-hr0        hr72-hr0         hr80-hr0
                     hr0                    -1              -1              -1              -1              -1
                     hr6                     1               0               0               0               0
                     hr24                    0               1               0               0               0
                     hr48                    0               0               1               0               0
                     hr72                    0               0               0               1               0
                     hr80                    0               0               0               0               1

exprs.row has 28 rows and an example corresponding to one ranked gene in the loop (code in the first post) is as follows:

exprs.row

hr0       hr0       hr6       hr6      hr24      hr24      hr48      hr48      hr72      hr72       hr0       hr0       hr6 
8.983808  9.157047  9.128908  9.029095  9.060012  8.962301  9.291968  9.309716 10.226955  9.950621  8.553083  8.907762  9.093225 
hr6      hr24      hr24      hr48      hr48      hr72      hr72      hr80      hr80       hr0       hr0       hr6       hr6 
8.984623  9.015506  8.989505  9.184224  9.019390  9.947579  9.941535  9.995935 10.010917  8.778475  8.847193  8.762086  8.743280 
hr72      hr72 
10.095114  9.926371

According to my replicate information, different time points has the following

0hr(1:2, 11:12, 23:24)
6hr(3:4, 13:14, 25:26)
24hr(5:6, 15:16)
48hr(7:8, 18:19)
72hr(9:10, 19:20, 27:28)
80hr(21:22)

If want to just have a single plot (per gene), I guess I need to average the exprs.row according to my replicate info, and then just plot time points against the averaged exprs.row to show a profile through the different time points (0, 6, 24, 48, 72 and 80h). I'm not sure if this is correct and how to write the codes for this new idea.

Apologies for the initial confusion.