Plot profiles for each group - spline analysis in edgeR
1
0
Entering edit mode
boczniak767 ▴ 740
@maciej-jonczyk-3945
Last seen 5 weeks ago
Poland

Hi,

this question relates to post Time-course, comparison of two groups. I have different question related to the data presented there so I've made new question.

I'd like to plot observed and fitted log-CPMs for a selected gene or genes.

I've consulted edgeR user guide (p. 113, the bottom block of code). But there is plot for one time series without replications. I'd like to make a plot for experiment linked above, there are two levels of "line".

I wonder if there is conveient way to plot such data computed by edgeR.

For now I assumed tahat averaging log-CPMs values is fine (???)

Here is my code (continuation of the one in the linked post)


# data arrangement

# id of the top gene
> tab = as.data.frame(topTags(fit2, n=30)) # here I used "fit2" object generated in the linked post
> ids=row.names(tab)[1]

# data for lines for fitted CPMs
> logCPM.fit <- cpm(fit2, log=TRUE) # according to edgeR user guide
> mf1=logCPM.fit[ids,1:12] # data for mutant
> mf2=logCPM.fit[ids,13:24] # data for wt
> mf12=rbind(mf1,mf2) # binding  rows for mutant and wt
> x=as.data.frame(mf12)
# averages or biological reps
> x$sr15=rowMeans(x[,1:3])
> x$sr22=rowMeans(x[,4:6])
> x$sr29=rowMeans(x[,7:9])
> x$sr36=rowMeans(x[,10:12])
> x2=x[,13:16] # selecting only averages
> mf12av=x2

# data for scatterplot of observed CPMs
> logCPM.obs <- cpm(cts.tmm, log=TRUE, prior.count=fit2$prior.count) # according to edgeR user guide
> m1=logCPM.obs[ids,1:12] # data for mutant
> m2=logCPM.obs[ids,13:24] # data for wt
> m12=rbind(m1,m2) # binding data for mutant and wt
> x=m12 # new object, just in case if I corrupt the original
# averages for biological reps
> x$sr15=rowMeans(x[,1:3])
> x$sr22=rowMeans(x[,4:6])
> x$sr29=rowMeans(x[,7:9])
> x$sr36=rowMeans(x[,10:12])
> m12av=x[,13:16] # selecting only averages
> x=t(m12av) # transposing to get data for matplot

# plot
> matplot(c(15, 22, 29, 36), x)
> lines(c(15, 22, 29, 36), mf12av[1,], col="black", lwd=2)
> lines(c(15, 22, 29, 36), mf12av[2,], col="red", lwd=2)

And here is my plot. It is obviously bad way of making such plots. Could anyone suggest me the better method?

timecoursedata edgeR • 523 views
ADD COMMENT

Login before adding your answer.

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