View Single Post
Old 03-17-2014, 02:11 AM   #31
shocker8786
Member
 
Location: Urbana Illinois

Join Date: Jan 2013
Posts: 28
Default

I retried the code and was still not successful with one of my datasets. I am sure that the code I am using is the same for both datasets (I have been copying and pasting the code for both datasets). Below is the code and the results for the commands you asked me to enter for the dataset that did not work properly.

> cuff.res=read.delim(file="gene_exp.diff", sep="\t")
> cuff.fc=cuff.res$log2.fold_change.
> gnames=cuff.res$gene
> sel=gnames!="-"
> gnames=as.character(gnames[sel])
> cuff.fc=cuff.fc[sel]
> names(cuff.fc)=gnames
> gnames.eg=pathview::id2eg(gnames, category ="symbol", org="Ss")
Loading required package: org.Ss.eg.db

> sel2=gnames.eg[,2]>""
> cuff.fc=cuff.fc[sel2]
> names(cuff.fc)=gnames.eg[sel2,2]
> range(cuff.fc)
[1] -Inf Inf
> cuff.fc[cuff.fc>10]=10
> cuff.fc[cuff.fc< -10]=-10
> exp.fc=cuff.fc
> out.suffix="cuff"
> require(gage)
> kg.ssc=kegg.gsets("ssc")
> kegg.gs=kg.ssc$kg.sets[kg.ssc$sigmet.idx]
> fc.kegg.p <- gage(exp.fc, gsets = kegg.gs, ref = NULL, samp = NULL)
> sel <- fc.kegg.p$greater[, "q.val"] < 0.1 &
+ !is.na(fc.kegg.p$greater[, "q.val"])
> path.ids <- rownames(fc.kegg.p$greater)[sel]
> sel.l <- fc.kegg.p$less[, "q.val"] < 0.1 &
+ !is.na(fc.kegg.p$less[, "q.val"])
> path.ids.l <- rownames(fc.kegg.p$less)[sel.l]
> path.ids2 <- substr(c(path.ids, path.ids.l), 1, 8)
> require(pathview)
> pv.out.list <- sapply(path.ids2[1:3], function(pid) pathview(
+ gene.data = exp.fc, pathway.id = pid,
+ species = "ssc", out.suffix=out.suffix))
Start tag expected, '<' not found
Parsing ./sscNA.xml file failed, please check the file!
Start tag expected, '<' not found
Parsing ./sscNA.xml file failed, please check the file!
Start tag expected, '<' not found
Parsing ./sscNA.xml file failed, please check the file!
> str(cuff.fc)
Named num [1:7434] 0.6134 0 0.3767 0.0177 -2.1022 ...
- attr(*, "names")= chr [1:7434] "100620481" "100512748" "100513321" "100154202" ...
> head(cuff.fc)
100620481 100512748 100513321 100154202 733660 100153485
0.6133590 0.0000000 0.3766890 0.0177489 -2.1022400 -0.2612780
> lapply(kegg.gs[1:3], head, 3)
$`ssc00970 Aminoacyl-tRNA biosynthesis`
[1] "100144457" "100153423" "100155115"

$`ssc02010 ABC transporters`
[1] "100048963" "100127449" "100144586"

$`ssc03008 Ribosome biogenesis in eukaryotes`
[1] "100151786" "100151831" "100152658"

> head(fc.kegg.p$greater)
p.geomean stat.mean p.val
ssc04064 NF-kappa B signaling pathway 0.02469658 1.990928 0.02469658
ssc04620 Toll-like receptor signaling pathway 0.03526210 1.828538 0.03526210
ssc04066 HIF-1 signaling pathway 0.04658115 1.693048 0.04658115
ssc04640 Hematopoietic cell lineage 0.04929811 1.670237 0.04929811
ssc04670 Leukocyte transendothelial migration 0.04995310 1.657766 0.04995310
ssc04610 Complement and coagulation cascades 0.07305357 1.467898 0.07305357
q.val set.size exp1
ssc04064 NF-kappa B signaling pathway 0.8234987 53 0.02469658
ssc04620 Toll-like receptor signaling pathway 0.8234987 58 0.03526210
ssc04066 HIF-1 signaling pathway 0.8234987 59 0.04658115
ssc04640 Hematopoietic cell lineage 0.8234987 49 0.04929811
ssc04670 Leukocyte transendothelial migration 0.8234987 63 0.04995310
ssc04610 Complement and coagulation cascades 0.8234987 44 0.07305357

Below are the results for the dataset that did work properly:

Working in directory /srv/mds01/shared/Epigenetics/Rnaseq/kyle_results/TJT/tophat_expression/cuffdiff
Writing image file ssc04151.cuff.png
Working in directory /srv/mds01/shared/Epigenetics/Rnaseq/kyle_results/TJT/tophat_expression/cuffdiff
Writing image file ssc04010.cuff.png
Working in directory /srv/mds01/shared/Epigenetics/Rnaseq/kyle_results/TJT/tophat_expression/cuffdiff
Writing image file ssc04080.cuff.png
> str(cuff.fc)
Named num [1:344025] 1.279 0 0.155 0.806 -0.958 ...
- attr(*, "names")= chr [1:344025] "100620481" "100512748" "100513321" "100154202" ...
> head(cuff.fc)
100620481 100512748 100513321 100154202 100156632 733660
1.278660 0.000000 0.154690 0.806045 -0.958121 -0.103950
> lapply(kegg.gs[1:3], head, 3)
$`ssc00970 Aminoacyl-tRNA biosynthesis`
[1] "100144457" "100153423" "100155115"

$`ssc02010 ABC transporters`
[1] "100048963" "100127449" "100144586"

$`ssc03008 Ribosome biogenesis in eukaryotes`
[1] "100151786" "100151831" "100152658"

> head(fc.kegg.p$greater)
p.geomean stat.mean p.val
ssc00770 Pantothenate and CoA biosynthesis 0.9995680 -3.988628 0.9995680
ssc00630 Glyoxylate and dicarboxylate metabolism 0.9995962 -3.942065 0.9995962
ssc00052 Galactose metabolism 0.9999012 -4.281375 0.9999012
ssc04740 Olfactory transduction 0.9999254 -4.576407 0.9999254
ssc00100 Steroid biosynthesis 0.9999490 -4.852319 0.9999490
ssc00970 Aminoacyl-tRNA biosynthesis 0.9999580 -4.829070 0.9999580
q.val set.size exp1
ssc00770 Pantothenate and CoA biosynthesis 1 10 0.9995680
ssc00630 Glyoxylate and dicarboxylate metabolism 1 11 0.9995962
ssc00052 Galactose metabolism 1 15 0.9999012
ssc04740 Olfactory transduction 1 12 0.9999254
ssc00100 Steroid biosynthesis 1 11 0.9999490
ssc00970 Aminoacyl-tRNA biosynthesis 1 12 0.9999580
shocker8786 is offline   Reply With Quote