Search
Question: Multiple samples on the same graph
0
gravatar for Luiz
4 months ago by
Luiz0
Brazil
Luiz0 wrote:

Hi, good morning!

I hope you have a nice weekend.

I'm working on pathview today and i can't put different matrixes on the graph to compare. Only one is working and my version is 1.14.0. Already run the demo code and everything is fine accord with vignette.

I'm trying to do this : https://bioconductor.org/packages/devel/bioc/vignettes/pathview/inst/doc/pathview.pdf pg14

Workflow

setwd("~/R")
require(dplyr)
require(DESeq2)
require(ggplot2)
require(ggrepel)
require(pheatmap)
library(pathview)
library(gage)
library(gageData)
data(kegg.sets.hs)
data(sigmet.idx.hs)
data(gse16873.d)

kegg.sets.hs = kegg.sets.hs[sigmet.idx.hs]
head(kegg.sets.hs, 3)
bckCountTable <- read.table("Matrix.txt", header=TRUE, row.names=1)
samples <- data.frame(row.names=c("R1", "R2", "R3","IR1","IR2", "IR3"), condition=as.factor(c(rep("R",3),rep("IR",3))))
samples
bckCDS <- DESeqDataSetFromMatrix(countData = bckCountTable, colData=samples, design=~condition)
bckCDS <- bckCDS[ rowSums(counts(bckCDS)) > 1, ]
bckCDS_1 <- DESeq(bckCDS)
res <- results(bckCDS_1, alpha=0.05, betaPrior = F, contrast=c("condition", "IR", "R"), pAdjustMethod = "BH")
res$symbol = mapIds(org.Hs.eg.db,
                    keys=row.names(res), 
                    column="SYMBOL",
                    keytype="REFSEQ",
                    multiVals="first")
res$entrez = mapIds(org.Hs.eg.db,
                    keys=row.names(res), 
                    column="ENTREZID",
                    keytype="REFSEQ",
                    multiVals="first")
res$name = mapIds(org.Hs.eg.db,
                    keys=row.names(res), 
                    column="GENENAME",
                    keytype="REFSEQ",
                    multiVals="first")
foldchanges = res$log2FoldChange
names(foldchanges) = res$entrez
head(foldchanges)
head(res,10)
foldchanges = res$log2FoldChange
names(foldchanges) = res$entrez
head(foldchanges)
keggres = gage(foldchanges, gsets=kegg.sets.hs, same.dir=TRUE)
lapply(keggres, head)

bckCountTable2 <- read.table("Matrix2.txt", header=TRUE, row.names=1)
samples2 <- data.frame(row.names=c("CTL1", "CTL2", "CTL3","HT1","HT2", "HT3"), condition=as.factor(c(rep("CTL",3),rep("HT",3))))
samples2
bckCDS2 <- DESeqDataSetFromMatrix(countData = bckCountTable2, colData=samples2, design=~condition)
bckCDS2 <- bckCDS2[ rowSums(counts(bckCDS2)) > 1, ]
bckCDS_12 <- DESeq(bckCDS2)
res2 <- results(bckCDS_12, alpha=0.05, betaPrior = F, contrast=c("condition", "HT", "CTL"), pAdjustMethod = "BH")
res2$symbol = mapIds(org.Hs.eg.db,
                    keys=row.names(res2), 
                    column="SYMBOL",
                    keytype="REFSEQ",
                    multiVals="first")
res2$entrez = mapIds(org.Hs.eg.db,
                    keys=row.names(res2), 
                    column="ENTREZID",
                    keytype="REFSEQ",
                    multiVals="first")
res2$name = mapIds(org.Hs.eg.db,
                    keys=row.names(res2), 
                    column="GENENAME",
                    keytype="REFSEQ",
                    multiVals="first")
foldchanges2 = res2$log2FoldChange
names(foldchanges2) = res2$entrez
head(foldchanges2)
head(res2,10)
foldchanges2 = res2$log2FoldChange
names(foldchanges2) = res2$entrez
head(foldchanges2)
keggres2 = gage(foldchanges2, gsets=kegg.sets.hs, same.dir=TRUE)
lapply(keggres2, head)

Then i proceed to:

pv.out <- pathview(gene.data = foldchanges, cpd.data = foldchanges2, pathway.id = "any", + species = "hsa", out.suffix = "any", keys.align = "y", + kegg.native = T, match.data = F, multi.state = T, same.layer = T) 

But i get a graph with only the gene.data, like cpd.data don't exist.

Thank you for your attention!

 

ADD COMMENTlink modified 4 months ago by Luo Weijun1.4k • written 4 months ago by Luiz0
0
gravatar for Luo Weijun
4 months ago by
Luo Weijun1.4k
United States
Luo Weijun1.4k wrote:

You are not calling pathview the right way. In your analysis, you have two fold change vectors from two analyses, i.e. foldchanges and foldchanges2. You should combined them into a single matrix with two columns, say fc.mat, then pass that into argument gene.data, something like:

pv.out <- pathview(gene.data = fc.mat, pathway.id = "any", species = "hsa", out.suffix = "any")

 

the row names of fc.mat should be gene ID (Entrez in your case).

note cpd.data argument is for compound or metabolite data, not for addition gene data (as you did in your code). For details of pathview, do:

?pathview

 

You may always use pathview web server as an alternative:

https://pathview.uncc.edu/

Pathview Web server provides a intuitive GUI and API to Pathview, a full pathway analysis workflow supporting multiple omics data and integrated analysis, and much more.

 

ADD COMMENTlink written 4 months ago by Luo Weijun1.4k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 318 users visited in the last hour