使用pd.genomewidesnp.6为GenomeWideSNP_6做注释

代码如下:

> head(dat)
  Composite.Element.REF Call Confidence
1         SNP_A-2131660    2     0.0126
2         SNP_A-1967418    2     0.0029
3         SNP_A-1969580    2     0.0019
4         SNP_A-4263484    2     0.0043
5         SNP_A-1978185    0     0.0082
6         SNP_A-4264431    1     0.0088
> library(BiocInstaller)
Bioconductor version 2.13 (BiocInstaller 1.12.0), ?biocLite for help
> biocLite("pd.genomewidesnp.6")
> library(pd.genomewidesnp.6)
> con <- pd.genomewidesnp.6@getdb()
> dbListTables(con)
 [1] "featureSet"        "featureSetCNV"     "fragmentLength"    "fragmentLengthCNV" "pmfeature"         "pmfeatureCNV"     
 [7] "sequence"          "sequenceCNV"       "sqlite_stat1"      "table_info"       
> dbGetQuery(con, "select * from featureSet limit 3")[,1:10]
  fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
1      1 SNP_A-2131660          NA   rs2887286     1      1156131      0   p36.33        C        T
2      2 SNP_A-1967418          NA   rs1496555     1      2234251      0   p36.33        A        G
3      3 SNP_A-1969580          NA  rs41477744     1      2329564      0   p36.32        A        G
> dim(dat)
[1] 906600      3
> anno <- dbGetQuery(con, "select * from featureSet")
> dim(anno)
[1] 906600     13
> data <- cbind(dat, anno[match(dat[,1], anno[,2]),])
> dim(data)
[1] 906600     16
> miR146a <- data[grepl("miR146a", data[,"gene_assoc"], ignore.case=T), ]
> dim(miR146a)
[1] 49 16
> head(data[,1:13])
  Composite.Element.REF Call Confidence fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
1         SNP_A-2131660    2     0.0126      1 SNP_A-2131660          NA   rs2887286     1      1156131      0   p36.33        C        T
2         SNP_A-1967418    2     0.0029      2 SNP_A-1967418          NA   rs1496555     1      2234251      0   p36.33        A        G
3         SNP_A-1969580    2     0.0019      3 SNP_A-1969580          NA  rs41477744     1      2329564      0   p36.32        A        G
4         SNP_A-4263484    2     0.0043      4 SNP_A-4263484          NA   rs3890745     1      2553624      0   p36.32        C        T
5         SNP_A-1978185    0     0.0082      5 SNP_A-1978185          NA  rs10492936     1      2936870      1   p36.32        C        T
6         SNP_A-4264431    1     0.0088      6 SNP_A-4264431          NA  rs10489588     1      2951834      1   p36.32        C        T
> miR146a[miR146a$dbsnp_rs_id=="rs2910164",1:13]
       Composite.Element.REF Call Confidence fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
283011         SNP_A-2264525    0     0.0012 283011 SNP_A-2264525          NA   rs2910164     5    159912418      1      q34        C        G

如何让R自动保存历史记录(auto save R history before crash)

也许第一眼看到这个题目,您可能并不认为这是一个问题,因为大家都知道,R会自动在当前工作目录保存名为.Rhistory或者.Rapp.history的文件。这个文件中就是R自动保存的历史记录。但是这种自动保存的历史记录有两个缺点,第一,它只在正常退出的时候保存;第二,它有行数限制。

我提出的这个问题,其实是因为我经常遇到R崩溃的情况,或者R无限耗用内存的情况。在这些情况下,R无法正常退出,也就[……]

Read more

如何获取R中代码的运行效率

对于如何测试R代码的运行效率,首先可能想到的就是proc.time()以及system.time()了。

我们来看代码:

> ptm <- proc.time()
> for (i in 1:50) mad(stats::runif(500))
> proc.time() - ptm
   user  system elapsed 
  0.022   0.002   0.038 
> system.time(for(i in 1:100) mad(runif(1000)))
   user  system elapsed 
  0.025   0.002   0.028

这两个函数都只适合测试一整块代码的运行效率。

还有可以想到的工具应该就是rbenchmark了。

> random.array = function(rows, cols, dist=rnorm) 
+                   array(dist(rows*cols), c(rows, cols))
> random.replicate = function(rows, cols, dist=rnorm)
+                       replicate(cols, dist(rows))
> benchmark(replications=rep(100, 3),
+           random.array(100, 100),
+           random.array(100, 100),
+           columns=c('test', 'elapsed', 'replications'))
                    test elapsed replications
1 random.array(100, 100)   0.093          100
2 random.array(100, 100)   0.094          100
3 random.array(100, 100)   0.093          100
4 random.array(100, 100)   0.096          100
5 random.array(100, 100)   0.092          100
6 random.array(100, 100)   0.093          100

rb[……]

Read more

在R中肩并肩(side-by-side)绘制箱线图(box-plot)

这里分为两部分,第一,肩并肩绘制图型,第二,对于有显著差异的相邻的数据标记星号。

> a<-matrix(nrow=100,ncol=3,data=runif(300,max=2))
> b<-matrix(nrow=100,ncol=3,data=runif(300,max=1))
> colnames(a)<-c("case 1","case 2","case 3")
> colnames(b)<-c("case 1","case 2","case 3")
> n <- 3
> xpos <- 0:(n-1)*3+1.5
> ypos.a <- apply(a, 2, max)
> ypos.b <- apply(b, 2, max)
> pvalue <- c(0.5, 0.05, 0.001)
> mark <- symnum(pvalue, cutpoints=c(0, 0.05, 1), symbols=c("*", NA))
> dist <- max(range(a,b))/20
> ylim <- range(a, b)
> ylim[2] <- ylim[2]+dist
> boxplot(a, at=0:(n-1)*3 + 1, xlim=c(0,n*3), ylim=ylim, xaxt="n", col="yellow")
> boxplot(b, at=0:(n-1)*3+2, xaxt="n", add=TRUE, col="red")
> axis(1, at = 0:(n-1)*3 + 1.5, labels = colnames(a), tick = TRUE)
> for(i in 1:length(mark)){
+ 	if(!is.na(mark[i])){
+ 		segments(xpos[i]-.5, ypos.a[i]+dist/2, xpos[i]-.5, max(ypos.a[i], ypos.b[i])+dist)
+ 		segments(xpos[i]+.5, ypos.b[i]+dist/2, xpos[i]+.5, max(ypos.a[i], ypos.b[i])+dist)
+ 		segments(xpos[i]-.5, max(ypos.a[i], ypos.b[i])+dist, xpos[i]-0.1, max(ypos.a[i], ypos.b[i])+dist)
+ 		segments(xpos[i]+.5, max(ypos.a[i], ypos.b[i])+dist, xpos[i]+0.1, max(ypos.a[i], ypos.b[i])+dist)
+ 		text(x=xpos[i], y=max(ypos.a[i], ypos.b[i])+dist, label=mark[i], col="red")
+ 	}
+ }

boxplot.sidebyside

使用R从网页中提取表格

有时候我们在数据挖掘的过程中需要提取网页中的表格。我们可以使用下面的代码:

> library(XML)
> url <- "http://www.bioguo.org/AnimalTFDB/BrowseAllTF.php?spe=Mus_musculus"
> tables <- readHTMLTable(url)
> length(tables)
[1] 6
> names(tables)
[1] "NULL"   "NULL"   "NULL"   "NULL"   "table1" "NULL"  
> dim(tables$table1)
[1] 1458    5
> head(tables$table1)
  No.         Ensembl ID Gene ID  Symbol            Family
1   1 ENSMUSG00000029313   17355    Aff1              AF-4
2   2 ENSMUSG00000031189   14266    Aff2              AF-4
3   3 ENSMUSG00000037138   16764    Aff3              AF-4
4   4 ENSMUSG00000049470   93736    Aff4              AF-4
5   5 ENSMUSG00000046532   11835      Ar Androgen receptor
6   6 ENSMUSG00000021359   21418 Tcfap2a              AP-2

从org.Xx.eg.db包中提取基因起止位置(本文已被历史淘汰,请使用TxDb包)

基因往往会有多个转录子版本,所以基因的起止位置并不一定。这里给出的示例代码用于演示如何从org.db包中提取出所有的基因范围。

library(org.Mm.eg.db)
##随机生成一些基因名
keys <- sample(keys(org.Mm.eg.db, keytype="ENTREZID"), 100)
genelist <- select(org.Mm.eg.db, keys=keys, cols="SYMBOL", keytype="ENTREZID")$SYMBOL
##查看数据库中的列名,以及key值名称
cols(org.Mm.eg.db)
keytypes(org.Mm.eg.db)
##获取基因转录子的起止点
location <- select(org.Mm.eg.db, keys=genelist, cols=c("CHR", "CHRLOC", "CHRLOCEND"), keytype="SYMBOL")
##或者间接方式
entrez <- select(org.Mm.eg.db, keys=genelist, cols=c("ENTREZID"), keytype="SYMBOL")$ENTREZID
location <- select(org.Mm.eg.db, keys=genelist, cols=c("SYMBOL", "CHR", "CHRLOC", "CHRLOCEND"), keytype="ENTREZID")

我们可以通过biomaRt来进行相类的操作

> library(biomaRt)
> mart <- useMart("ensembl", "mmusculus_gene_ensembl")
> loc <- getBM(c("mgi_symbol", "chromosome_name", "start_position", "end_position", "strand"), "mgi_symbol", genelist, mart)
> loc2 <- getBM(c("entrezgene", "chromosome_name", "start_position", "end_position", "strand"), "entrezgene", keys, mart)