Protein Expression Analysis with Uniprot IDs

#setwd('C:/Users/animeshs/SkyDrive')
#source("http://bioconductor.org/biocLite.R")
#install.packages("pvclust")
#install.packages("knitr")
#rm(list = ls())
#dev.off() 
#reporting
#library("knitr")
#clustering
library(gage)
library(pathview)
## Loading required package: KEGGgraph
## Loading required package: XML
## Loading required package: graph
## 
## Attaching package: 'graph'
## 
## The following object is masked from 'package:XML':
## 
##     addNode
## 
## Loading required package: org.Hs.eg.db
## Loading required package: AnnotationDbi
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## 
## The following object is masked from 'package:stats':
## 
##     xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, as.vector, cbind,
##     colnames, duplicated, eval, evalq, Filter, Find, get,
##     intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rep.int, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unlist
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Loading required package: DBI
## 
## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html.
## 
## The pathview downloads and uses KEGG data. Academic users may freely use the
## KEGG website at http://www.kegg.jp/ or its mirror site at GenomeNet
## http://www.genome.jp/kegg/. Academic users may also freely link to the KEGG
## website. Non-academic users may use the KEGG website as end users for
## non-commercial purposes, but any other use requires a license agreement
## (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################
# with MM20CL14 SS data over Hedgehog signaling pathway
hda=readExpData("L:/Elite/Aida/MM20MQv1p5XcompNewRD.txt",row.names=1)
hda=as.matrix(hda)
#cancer related hsa05200
pathview(hda,pathway.id="hsa05200",gene.idtype="UNIPROT",out.suffix="jak")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa05200.jak.multi.png
## [1] "some node width is different from others, and hence adjusted!"
# http://www.genome.jp/kegg/pathway/map/map00020.html
pathview(hda,pathway.id="hsa00020",gene.idtype="UNIPROT",out.suffix="TCA")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa00020.TCA.multi.png
# http://www.genome.jp/kegg/pathway/map/map00030.html
pathview(hda,pathway.id="hsa00030",gene.idtype="UNIPROT",out.suffix="PPP")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa00030.PPP.multi.png
# http://www.genome.jp/kegg/pathway/map/map00190.html 
pathview(hda,pathway.id="hsa00190",gene.idtype="UNIPROT",out.suffix="OxPhos")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa00190.OxPhos.multi.png
## [1] "some node width is different from others, and hence adjusted!"
# http://www.genome.jp/kegg/pathway/map/map00010.html
pathview(hda,pathway.id="hsa00010",gene.idtype="UNIPROT",out.suffix="Glycolysis")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa00010.Glycolysis.multi.png
data=nrow(hda)
dm=10
km=100
kegg=60
GOmap <-  matrix(c(data-dm, dm,kegg-dm, km), nrow = 2, dimnames = list(DATA = c("mapped", "unmapped"), KEGG = c("mapped", "unmapped")))
fisher.test(GOmap, alternative = "greater")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  GOmap
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  557.3   Inf
## sample estimates:
## odds ratio 
##       1029
library("pvclust")
result <- pvclust(hda, method.dist="cor", method.hclust="average", nboot=10)
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
plot(result)
pvrect(result, alpha=0.95)

plot of chunk cluster bootstrap

hda=readExpData("C:/Users/animeshs/OneDrive/SHBER.txt",row.names=1)
## Warning: cannot open file 'C:/Users/animeshs/OneDrive/SHBER.txt': No such
## file or directory
## Error: cannot open the connection
hda=as.matrix(hda)
hda.d=hda[,1:3]-hda[,4:6]
summary(hda.d)
##       MM1            MM2            MM3      
##  Min.   :-4.9   Min.   :-4.1   Min.   :-4.1  
##  1st Qu.:-0.3   1st Qu.:-0.4   1st Qu.:-0.4  
##  Median : 0.0   Median :-0.1   Median : 0.0  
##  Mean   : 0.0   Mean   :-0.1   Mean   : 0.0  
##  3rd Qu.: 0.4   3rd Qu.: 0.3   3rd Qu.: 0.4  
##  Max.   : 3.8   Max.   : 3.3   Max.   : 4.0  
##  NA's   :2650   NA's   :2550   NA's   :2598
pv.out<-pathview(hda.d,pathway.id="hsa03410",gene.idtype="UNIPROT", limit = list(gene = 5, cpd = 1), out.suffix=proc.time())
## [1] "Downloading xml files for hsa03410, 1/1 pathways.."
## [1] "Downloading png files for hsa03410, 1/1 pathways.."
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa03410.43.22.multi.png
str(pv.out)
## List of 2
##  $ plot.data.gene:'data.frame':  43 obs. of  13 variables:
##   ..$ kegg.names: chr [1:43] "328" "2237" "142" "3978" ...
##   ..$ labels    : chr [1:43] "APEX1" "FEN1" "PARP1" "LIG1" ...
##   ..$ type      : chr [1:43] "gene" "gene" "gene" "gene" ...
##   ..$ x         : num [1:43] 909 955 955 863 954 955 909 863 862 681 ...
##   ..$ y         : num [1:43] 957 957 974 974 786 991 991 957 794 610 ...
##   ..$ width     : num [1:43] 46 46 46 46 46 46 46 46 46 46 ...
##   ..$ height    : num [1:43] 17 17 17 17 17 17 17 17 17 17 ...
##   ..$ MM1       : num [1:43] -0.338 0.457 -0.307 NA -0.545 ...
##   ..$ MM2       : num [1:43] -0.0473 0.3756 0.4787 NA -0.2693 ...
##   ..$ MM3       : num [1:43] 1.398 -0.894 -1.588 NA 1.501 ...
##   ..$ MM1.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM2.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM3.col   : Factor w/ 4 levels "#8FCE8F","#BEBEBE",..: 3 2 1 4 3 2 2 2 2 2 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##  $ plot.data.cpd : NULL
head(pv.out$plot.data.gene)
##   kegg.names labels type   x   y width height     MM1      MM2     MM3
## 1        328  APEX1 gene 909 957    46     17 -0.3381 -0.04732  1.3979
## 2       2237   FEN1 gene 955 957    46     17  0.4573  0.37563 -0.8937
## 3        142  PARP1 gene 955 974    46     17 -0.3071  0.47869 -1.5877
## 4       3978   LIG1 gene 863 974    46     17      NA       NA      NA
## 5       3980   LIG3 gene 954 786    46     17 -0.5451 -0.26930  1.5015
## 7      54107  POLE3 gene 955 991    46     17 -0.2837 -0.08812  0.0000
##   MM1.col MM2.col MM3.col
## 1 #BEBEBE #BEBEBE #CE8F8F
## 2 #BEBEBE #BEBEBE #BEBEBE
## 3 #BEBEBE #BEBEBE #8FCE8F
## 4 #FFFFFF #FFFFFF #FFFFFF
## 5 #BEBEBE #BEBEBE #CE8F8F
## 7 #BEBEBE #BEBEBE #BEBEBE
# test on http://www.nature.com/nri/journal/v5/n5/fig_tab/nri1604_F1.html 
pathview(hda.d,pathway.id="hsa04630",gene.idtype="UNIPROT",out.suffix="jak")
## [1] "Downloading xml files for hsa04630, 1/1 pathways.."
## [1] "Downloading png files for hsa04630, 1/1 pathways.."
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa04630.jak.multi.png
pathview(hda.d,pathway.id="hsa00020",gene.idtype="UNIPROT",out.suffix="joxphos")
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa00020.joxphos.multi.png
# Trasfection
hda=readExpData("L:/Elite/LARS/2014/mai/transfection 3rd paralell/ListPV1in20FC2.txt",row.names=1)
hda=as.matrix(hda)
# herpes simplex
pathview(hda,pathway.id="hsa03410",gene.idtype="UNIPROT")
## No of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa03410.pathview.png
pathview(hda,pathway.id="hsa05168",gene.idtype="UNIPROT")
## [1] "Downloading xml files for hsa05168, 1/1 pathways.."
## [1] "Downloading png files for hsa05168, 1/1 pathways.."
## No of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa05168.pathview.png
library("org.Hs.eg.db")
frame = toTable(org.Hs.egGO)
goframeData = data.frame(frame$go_id, frame$Evidence, frame$gene_id)
head(goframeData)
##   frame.go_id frame.Evidence frame.gene_id
## 1  GO:0008150             ND             1
## 2  GO:0001869            IDA             2
## 3  GO:0002576            TAS             2
## 4  GO:0007596            TAS             2
## 5  GO:0007264            TAS             2
## 6  GO:0007597            TAS             2
goFrame=GOFrame(goframeData,organism="Homo sapiens")
## Loading required package: GO.db
goAllFrame=GOAllFrame(goFrame)
library("GSEABase")
## Loading required package: annotate
gsc <- GeneSetCollection(goAllFrame, setType = GOCollection())
library("GOstats")
## Loading required package: Category
## Loading required package: Matrix
## 
## Attaching package: 'GOstats'
## 
## The following object is masked from 'package:AnnotationDbi':
## 
##     makeGOGraph
universe = Lkeys(org.Hs.egGO)
genes = universe[1:500]
params <- GSEAGOHyperGParams(name="My Custom GSEA based annot Params",
                              geneSetCollection=gsc,
                              geneIds = genes,
                              universeGeneIds = universe,
                              ontology = "MF",
                              pvalueCutoff = 0.05,
                              conditional = FALSE,
                              testDirection = "over")
Over <- hyperGTest(params)
head(summary(Over))
##       GOMFID    Pvalue OddsRatio ExpCount Count Size
## 1 GO:0016805 0.0001768    127.58 0.020074     2   14
## 2 GO:0019239 0.0007253     58.83 0.040149     2   28
## 3 GO:0004180 0.0012680     43.67 0.053054     2   37
## 4 GO:0033882 0.0014339       Inf 0.001434     1    1
## 5 GO:0050129 0.0014339       Inf 0.001434     1    1
## 6 GO:0004060 0.0028658    729.52 0.002868     1    2
##                                     Term
## 1                   dipeptidase activity
## 2                     deaminase activity
## 3              carboxypeptidase activity
## 4         choloyl-CoA hydrolase activity
## 5 N-formylglutamate deformylase activity
## 6 arylamine N-acetyltransferase activity
pv.out <- pathview(hda.d,pathway.id="hsa03410",gene.idtype="UNIPROT", limit = list(gene = 5, cpd = 1), out.suffix="fcnn", kegg.native = TRUE)
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa03410.fcnn.multi.png
str(pv.out)
## List of 2
##  $ plot.data.gene:'data.frame':  43 obs. of  13 variables:
##   ..$ kegg.names: chr [1:43] "328" "2237" "142" "3978" ...
##   ..$ labels    : chr [1:43] "APEX1" "FEN1" "PARP1" "LIG1" ...
##   ..$ type      : chr [1:43] "gene" "gene" "gene" "gene" ...
##   ..$ x         : num [1:43] 909 955 955 863 954 955 909 863 862 681 ...
##   ..$ y         : num [1:43] 957 957 974 974 786 991 991 957 794 610 ...
##   ..$ width     : num [1:43] 46 46 46 46 46 46 46 46 46 46 ...
##   ..$ height    : num [1:43] 17 17 17 17 17 17 17 17 17 17 ...
##   ..$ MM1       : num [1:43] -0.338 0.457 -0.307 NA -0.545 ...
##   ..$ MM2       : num [1:43] -0.0473 0.3756 0.4787 NA -0.2693 ...
##   ..$ MM3       : num [1:43] 1.398 -0.894 -1.588 NA 1.501 ...
##   ..$ MM1.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM2.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM3.col   : Factor w/ 4 levels "#8FCE8F","#BEBEBE",..: 3 2 1 4 3 2 2 2 2 2 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##  $ plot.data.cpd : NULL
head(pv.out$plot.data.gene)
##   kegg.names labels type   x   y width height     MM1      MM2     MM3
## 1        328  APEX1 gene 909 957    46     17 -0.3381 -0.04732  1.3979
## 2       2237   FEN1 gene 955 957    46     17  0.4573  0.37563 -0.8937
## 3        142  PARP1 gene 955 974    46     17 -0.3071  0.47869 -1.5877
## 4       3978   LIG1 gene 863 974    46     17      NA       NA      NA
## 5       3980   LIG3 gene 954 786    46     17 -0.5451 -0.26930  1.5015
## 7      54107  POLE3 gene 955 991    46     17 -0.2837 -0.08812  0.0000
##   MM1.col MM2.col MM3.col
## 1 #BEBEBE #BEBEBE #CE8F8F
## 2 #BEBEBE #BEBEBE #BEBEBE
## 3 #BEBEBE #BEBEBE #8FCE8F
## 4 #FFFFFF #FFFFFF #FFFFFF
## 5 #BEBEBE #BEBEBE #CE8F8F
## 7 #BEBEBE #BEBEBE #BEBEBE
## MS to SNP

source("http://bioconductor.org/biocLite.R")
## Bioconductor version 2.13 (BiocInstaller 1.12.1), ?biocLite for help
## A newer version of Bioconductor is available after installing a new
##   version of R, ?BiocUpgrade for help
biocLite("sapFinder")
## BioC_mirror: http://bioconductor.org
## Using Bioconductor version 2.13 (BiocInstaller 1.12.1), R version 3.0.3.
## Installing package(s) 'sapFinder'
## Warning: package 'sapFinder' is not available (for R version 3.0.3)
## Old packages: 'nleqslv'
library("sapFinder")
## Error: there is no package called 'sapFinder'
browseVignettes("sapFinder")
## Error: there is no package called 'sapFinder'
vcf <- system.file("extdata/sapFinder_test.vcf",
                   package="sapFinder")
annotation <- system.file("extdata/sapFinder_test_ensGene.txt",
                          package="sapFinder")
refseq <- system.file("extdata/sapFinder_test_ensGeneMrna.fa",
                      package="sapFinder")
outdir <- "db_dir"
prefix <- "sapFinder_test"
db.files <- dbCreator(vcf=vcf, annotation=annotation,
                      refseq=refseq, outdir=outdir,prefix=prefix)
## Error: could not find function "dbCreator"
#DEMO

#load data

#KEGG view: gene data only
i <- 1
pv.out <- pathview(gene.data = gse16873.d[, 1], pathway.id =
                     demo.paths$sel.paths[i], species = "hsa", out.suffix = "gse16873",
                   kegg.native = TRUE)
## Error: object 'gse16873.d' not found
str(pv.out)
## List of 2
##  $ plot.data.gene:'data.frame':  43 obs. of  13 variables:
##   ..$ kegg.names: chr [1:43] "328" "2237" "142" "3978" ...
##   ..$ labels    : chr [1:43] "APEX1" "FEN1" "PARP1" "LIG1" ...
##   ..$ type      : chr [1:43] "gene" "gene" "gene" "gene" ...
##   ..$ x         : num [1:43] 909 955 955 863 954 955 909 863 862 681 ...
##   ..$ y         : num [1:43] 957 957 974 974 786 991 991 957 794 610 ...
##   ..$ width     : num [1:43] 46 46 46 46 46 46 46 46 46 46 ...
##   ..$ height    : num [1:43] 17 17 17 17 17 17 17 17 17 17 ...
##   ..$ MM1       : num [1:43] -0.338 0.457 -0.307 NA -0.545 ...
##   ..$ MM2       : num [1:43] -0.0473 0.3756 0.4787 NA -0.2693 ...
##   ..$ MM3       : num [1:43] 1.398 -0.894 -1.588 NA 1.501 ...
##   ..$ MM1.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM2.col   : Factor w/ 2 levels "#BEBEBE","#FFFFFF": 1 1 1 2 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##   ..$ MM3.col   : Factor w/ 4 levels "#8FCE8F","#BEBEBE",..: 3 2 1 4 3 2 2 2 2 2 ...
##   .. ..- attr(*, "names")= chr [1:43] "1" "2" "3" "4" ...
##  $ plot.data.cpd : NULL
head(pv.out$plot.data.gene)
##   kegg.names labels type   x   y width height     MM1      MM2     MM3
## 1        328  APEX1 gene 909 957    46     17 -0.3381 -0.04732  1.3979
## 2       2237   FEN1 gene 955 957    46     17  0.4573  0.37563 -0.8937
## 3        142  PARP1 gene 955 974    46     17 -0.3071  0.47869 -1.5877
## 4       3978   LIG1 gene 863 974    46     17      NA       NA      NA
## 5       3980   LIG3 gene 954 786    46     17 -0.5451 -0.26930  1.5015
## 7      54107  POLE3 gene 955 991    46     17 -0.2837 -0.08812  0.0000
##   MM1.col MM2.col MM3.col
## 1 #BEBEBE #BEBEBE #CE8F8F
## 2 #BEBEBE #BEBEBE #BEBEBE
## 3 #BEBEBE #BEBEBE #8FCE8F
## 4 #FFFFFF #FFFFFF #FFFFFF
## 5 #BEBEBE #BEBEBE #CE8F8F
## 7 #BEBEBE #BEBEBE #BEBEBE
#result PNG file in current directory

#Graphviz view: gene data only
data(gse16873.d)
data(demo.paths)
pv.out <- pathview(gene.data = gse16873.d[, 1], pathway.id = demo.paths$sel.paths[i], species = "hsa", out.suffix = "gse16873",kegg.native = FALSE, sign.pos = demo.paths$spos[i])
## [1] "Downloading xml files for hsa04110, 1/1 pathways.."
## [1] "Downloading png files for hsa04110, 1/1 pathways.."
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa04110.gse16873.pdf
#result PDF file in current directory

png('test.png')
plot(sin(seq(-10,10,0.1)),type = "l", cex = .1, col = "dark red")
dev.off()
## pdf 
##   2
?plot
## starting httpd help server ... done
# from Luo
hda=as.matrix(hda)
hda.d=hda[,1:3]-hda[,4:6]
pathview(hda.d,pathway.id="hsa03410",gene.idtype="UNIPROT", limit = list(gene = 5, cpd = 1), out.suffix="fc")
## No of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Working in directory C:/ani/home/animeshs/misccb
## Writing image file hsa03410.fc.png
## PCA
data=read.csv("L:/Elite/gaute/test/SHBER.csv")
gene=as.character(data[,1])
data=data[,-1]
prco=prcomp(data)
sco <- prco$x
loa <- prco$rotation
plot(sco[,1], sco[,2], col=row.names(data))

plot of chunk check pathview for base excision repair pathway

plot(loa[,1], loa[,2], col=row.names(data))

plot of chunk check pathview for base excision repair pathway

library(org.Hs.eg.db)
egmap <- revmap(org.Hs.egUNIPROT)[as.character(gene)]

x <- org.Hs.egGO
mapped_genes <- mappedkeys(x)
xx <- as.list(x[mapped_genes])
goids <- xx[2:3]
names(goids[[1]])
## [1] "GO:0006805" "GO:0044281" "GO:0005829" "GO:0004060"
mod <- prcomp(E)
## Error: object 'E' not found
## scores (sample space)
scores <- mod$x
## Error: object 'mod' not found
## loadings (gene space)
loads <- mod$rotation
## Error: object 'mod' not found
## gene symbols
sym <- as.character(mget(colnames(E), illuminaHumanv3SYMBOL))
## Error: error in evaluating the argument 'x' in selecting a method for function 'mget': Error in colnames(E) : 
##   error in evaluating the argument 'x' in selecting a method for function 'colnames': Error: object 'E' not found
# plot
pc1 = 1 # ok, checked
pc2 = 2
# score plot (sample space)
x11()
expVar <- 100 *mod$sdev / sum(mod$sdev) 
## Error: object 'mod' not found
hist(expVar)
## Error: object 'expVar' not found
plot(scores[,pc1], scores[,pc2], col=cl, pch=19, main="Score-plot (Sample map)", xlab=paste("PC:", expVar[pc1]), ylab=paste("PC:", expVar[pc2]))
## Error: error in evaluating the argument 'x' in selecting a method for function 'plot': Error: object 'scores' not found
text(scores[,pc1], scores[,pc2], rownames(E), pos=3)
## Error: object 'scores' not found
# loadings (gene space)
plot.factor.name <- "lrT - T"
col <- factor(a[,plot.factor.name])
## Error: object 'a' not found
plot(loads[,pc1], loads[,pc2], col=col, pch='.', main="Gene weights", xlab=paste("PC:", expVar[pc1]), ylab=paste("PC:", expVar[pc2]))
## Error: error in evaluating the argument 'x' in selecting a method for function 'plot': Error: object 'loads' not found
text(loads[,pc1], loads[,pc2], sym, pos=3, cex=0.4)
## Error: object 'loads' not found
# extract a quad 
kk <- loads[,pc1] > 0 & loads[,pc2] > 0
## Error: object 'loads' not found
plot(kk)
## Error: error in evaluating the argument 'x' in selecting a method for function 'plot': Error: object 'kk' not found
# GO mapping
x <- org.Hs.egGO
mapped_genes <- mappedkeys(x)
xx <- as.list(x[mapped_genes])
goids <- xx[2:3]
names(goids[[1]])
## [1] "GO:0006805" "GO:0044281" "GO:0005829" "GO:0004060"

You can also embed plots, for example:

plot(cars)

plot of chunk unnamed-chunk-1

set.seed(1213)  # for reproducibility
x = cumsum(rnorm(100))
mean(x)  # mean of x
## [1] -1.94
## [1] -1.939758
plot(x, type = 'l')  # Brownian motion

plot of chunk generate randome graph for test

x = cumsum(rnorm(100))
mean(x)  # mean of x
## [1] -0.2687