木有，但是看过例子
你的问题多多，不知道这个例子够不够～
> library(Biobase)
Installing Bioconductor
> source("http://bioconductor.org/biocLite.R")
> # install the core Bioconductor packages
> biocLite()
> # install one or more specific packages
> biocLite(c("limma", "GenomicRanges"))
> # upgrade all currently installed packages
> biocLite(character())
Loading Raw Data
> library(affy)
> targets <- readTargets()
> targets
> # filename dmnt1 celltype sample.label
> # 1 (Mouse430_2)_WT1.CEL.gz wt hsc HSC +/+ 1
> # 2 (Mouse430_2)_WT3.CEL.gz wt hsc HSC +/+ 2
> # 3 (Mouse430_2)_WT7.CEL.gz wt mp MP +/+ 1
> # 4 (Mouse430_2)_WT8.CEL.gz wt mp MP +/+ 2
> # 5 (Mouse430_2)_WT9.CEL.gz wt mp MP +/+ 3
> # 6 (Mouse430_2)__2.CEL.gz wt b BP +/+ 1
> # 7 (Mouse430_2)_3.CEL.gz wt b BP +/+ 2
> # 8 GSM371218.CEL.gz wt etp ETP 1
> # 9 GSM371219.CEL.gz wt etp ETP 2
> # 10 GSM371220.CEL.gz wt etp ETP 3
> raw <- ReadAffy(filenames=targets$filename,
> celfile.path="./input/raw",
> sampleNames=targets$sample.label,
> phenoData=targets,
> annotation="mouse4302.db")
Briey on Annotation in R
> library(hgu95av2.db)
> library(annotate)
> mget(c("738_at", "40840_at", "41668_r_at"), hgu95av2GENENAME)
> hgu95av2() # includes info about when the package was updated
Working with ExpressionSets
> # Load an ExpressionSet object
> library(ALL)
> data(ALL)
> pData(ALL)[1:3, 1:3]
## cod diagnosis sex
## 01005 1005 5/21/1997 M
## 01010 1010 3/29/2000 M
## 03002 3002 6/24/1998 F
> annotation(ALL)
## [1] "hgu95av2"
> sampleNames(ALL)[1:10]
## [1] "01005" "01010" "03002" "04006" "04007" "04008" "04010" "04016"
## [9] "06002" "08001"
> # Pull out the matrix of expression values
> x <- exprs(ALL)
> # Can be subset like a data frame, but keeps associated annotation
> # information
> ALL_subset <- ALL[, 1:5]
> head(pData(ALL_subset))[, c("diagnosis", "sex", "age")]
## diagnosis sex age
## 01005 5/21/1997 M 53
## 01010 3/29/2000 M 19
## 03002 6/24/1998 F 52
## 04006 7/17/1997 M 38
## 04007 7/22/1997 M 57
Quality assessment
> library(arrayQualityMetrics)
> # This will create a report in the subdirectory qa-norm
> arrayQualityMetrics(eset, outdir = "./qa-norm/")
> # Run it on both normalized and raw data
> arrayQualityMetrics(raw, outdir = "./qa-raw/")
Di erential Expression Analysis
> library(limma) > limmaUsersGuide()
Functional Group Enrichment
> library(GOstats)
> # Define the universe first, then create a params object
> params <- new("GOHyperGParams",
> geneIds=selectedEntrezIds,
> universeGeneIds=myUniverse,
> annotation="moe430v2.db",
>,
> pvalueCutoff=0.01,
> conditional=FALSE,
> testDirection="over")
> hgOver <- hyperGTest(params)
> summary(hgOver)
Beyond Microarrays: ChIP-seq
> library(GenomicRanges)
> data <- read.table("test.bed", header=FALSE)
> colnames(data) <- c('chr','start','end','id','score','strand')
> gr <- with(data, GRanges(chr, IRanges(start, end),
> strand, score, id=id))
> # Keep ranges that overlap with the ranges
> # in another GRanges object:
> overlaps <- subsetByOverlap(gr, other_gr)
Beyond Microarrays: Sequence data
> library(BSgenome.Hsapiens.UCSC.hg19)
> # Specify chromosome, start and end
> seqs <- getSeq(Hsapiens, "chr1", 100, 10000)
> # Can also take a GRanges object
> seqs2 <- getSeq(Hsapiens, gr)
> # Has precomputed upstream sequences as well
> upstream2k <- Hsapiens$upstream2000
扫码加好友，拉您进群
全部版块
我的主页
收藏
