Off topic:FW: EdgeR GLM
0
0
Entering edit mode
Shona Wood ▴ 70
@shona-wood-4559
Last seen 9.6 years ago
Hi, I am trying to look at gene expression changes with increasing age. I have three time points young, middle age and old. I have been using the following, ordered factors and coefficient in order to get a fold change and FDR for a linear change in gene expression with increasing age. Can you look over my code and see if I have done this right? Is age.L the fold change in gene expression with increasing age? > library(edgeR) > library(limma) > targets<-read.delim (file="coding_targets.txt") > targets$age<-factor(ordered (targets$age)) > targets$sample<-factor(targets$sample) > targets X files age sample 1 p1 p1_coding_counts.txt a-young 1 2 p2 p2_coding_counts.txt a-young 2 3 p3 p3_coding_counts.txt a-young 3 4 p7 p7_coding_counts.txt b-middle 7 5 p8 p8_coding_counts.txt b-middle 8 6 p9 p9_coding_counts.txt b-middle 9 7 p4 p4_coding_counts.txt c-old 4 8 p5 p5_coding_counts.txt c-old 5 9 p6 p6_coding_counts.txt c-old 6 > d<-readDGE(targets, skip=1, comment.char='#') > colnames(d)<-row.names(targets) > d<- calcNormFactors(d) > d An object of class "DGEList" $samples X files age sample lib.size norm.factors 1 p1 p1_coding_counts.txt a-young 1 3445622 0.9655724 2 p2 p2_coding_counts.txt a-young 2 2696547 0.9902573 3 p3 p3_coding_counts.txt a-young 3 3308099 0.9787044 4 p7 p7_coding_counts.txt b-middle 7 2503479 1.0584660 5 p8 p8_coding_counts.txt b-middle 8 2639127 1.0477893 6 p9 p9_coding_counts.txt b-middle 9 2696547 0.9902573 7 p4 p4_coding_counts.txt c-old 4 3037440 0.9778553 8 p5 p5_coding_counts.txt c-old 5 2647915 1.0144109 9 p6 p6_coding_counts.txt c-old 6 2475370 0.9809077 $counts 1 2 3 4 5 6 7 8 9 ENSRNOG00000000007 1287 1285 1041 788 968 1285 1092 1009 960 ENSRNOG00000000008 0 0 0 0 1 0 0 1 0 ENSRNOG00000000009 0 0 0 3 0 0 1 0 0 ENSRNOG00000000010 38 405 50 18 105 405 372 42 282 ENSRNOG00000000012 0 0 0 0 0 0 0 0 0 22932 more rows ... > d<-d[rowSums(d$counts)>9,] > design<- model.matrix(~ age, data = targets) > design (Intercept) age.L age.Q 1 1 -7.071068e-01 0.4082483 2 1 -7.071068e-01 0.4082483 3 1 -7.071068e-01 0.4082483 4 1 -7.850462e-17 -0.8164966 5 1 -7.850462e-17 -0.8164966 6 1 -7.850462e-17 -0.8164966 7 1 7.071068e-01 0.4082483 8 1 7.071068e-01 0.4082483 9 1 7.071068e-01 0.4082483 attr(,"assign") [1] 0 1 1 attr(,"contrasts") attr(,"contrasts")$age [1] "contr.poly" > d<-estimateGLMCommonDisp(d, design) > names(d) [1] "samples" "counts" "common.dispersion" > glmfit<- glmFit(d, design, dispersion=d$common.dispersion) > results<- glmLRT(d, glmfit, coef=c(2,3)) > topTags(results) Coefficient: age.L age.Q logConc age.L age.Q LR P.Value ENSRNOG00000011672 -11.27122 -3.2445886 -1.9642177 109.53641 1.638591e-24 Many Thanks [[alternative HTML version deleted]]
• 920 views
ADD COMMENT
This thread is not open. No new answers may be added
Traffic: 761 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