使用R绘制带比例的饼图 Pie Chart with Annotated Percentages

pie1 <- function (x, labels = names(x), edges = 200, radius = 0.8, clockwise = FALSE, 
    init.angle = if (clockwise) 90 else 0, density = NULL, angle = 45, 
    col = NULL, border = NULL, lty = NULL, main = NULL, percentage=T, rawNumber=F, digits=3, cutoff=0.01, 
    legend=F, legendpos="topright", legendcol=2, ...) 
{
    if (!is.numeric(x) || any(is.na(x) | x < 0)) 
        stop("'x' values must be positive.")
    if (is.null(labels)) 
        labels <- as.character(seq_along(x))
    else labels <- as.graphicsAnnot(labels)
    rawX <- x
    x <- c(0, cumsum(x)/sum(x))
    dx <- diff(x)
    nx <- length(dx)
    plot.new()
    pin <- par("pin")
    xlim <- ylim <- c(-1, 1)
    if (pin[1L] > pin[2L]) 
        xlim <- (pin[1L]/pin[2L]) * xlim
    else ylim <- (pin[2L]/pin[1L]) * ylim
    dev.hold()
    on.exit(dev.flush())
    plot.window(xlim, ylim, "", asp = 1)
    if (is.null(col)) 
        col <- if (is.null(density)) 
            c("white", "lightblue", "mistyrose", "lightcyan", 
                "lavender", "cornsilk", "pink")
        else par("fg")
    if (!is.null(col)) 
        col <- rep_len(col, nx)
    if (!is.null(border)) 
        border <- rep_len(border, nx)
    if (!is.null(lty)) 
        lty <- rep_len(lty, nx)
    angle <- rep(angle, nx)
    if (!is.null(density)) 
        density <- rep_len(density, nx)
    twopi <- if (clockwise) 
        -2 * pi
    else 2 * pi
    t2xy <- function(t) {
        t2p <- twopi * t + init.angle * pi/180
        list(x = radius * cos(t2p), y = radius * sin(t2p))
    }
    for (i in 1L:nx) {
        n <- max(2, floor(edges * dx[i]))
        P <- t2xy(seq.int(x[i], x[i + 1], length.out = n))
        polygon(c(P$x, 0), c(P$y, 0), density = density[i], angle = angle[i], 
            border = border[i], col = col[i], lty = lty[i])
        if(!legend){
        	P <- t2xy(mean(x[i + 0:1]))
	        lab <- as.character(labels[i])
	        if (!is.na(lab) && nzchar(lab)) {
	            lines(c(1, 1.05) * P$x, c(1, 1.05) * P$y)
	            text(1.1 * P$x, 1.1 * P$y, labels[i], xpd = TRUE, 
	                adj = ifelse(P$x < 0, 1, 0), ...)
	        }
        }
    }
    if (percentage) {
    	for (i in 1L:nx){
    		if(dx[i]>cutoff){
    			P <- t2xy(mean(x[i + 0:1]))
            	text(.8 * P$x, .8 * P$y, paste(formatC(dx[i]*100, digits=digits), "%", sep=""), xpd = TRUE, 
                	adj = .5, ...)
    		}
        }
    }else{
        if(rawNumber){
		for (i in 1L:nx){
    			if(dx[i]>cutoff){
    				P <- t2xy(mean(x[i + 0:1]))
            		text(.8 * P$x, .8 * P$y, rawX[i], xpd = TRUE, 
                		adj = .5, ...)
    			}
        	}
        }
    }
    if(legend) legend(legendpos, legend=labels, fill=col, border="black", bty="n", ncol = legendcol)
    title(main = main, ...)
    invisible(NULL)
}
 
pie1(1:5)

pie

每天更新R原码 svn R source

在一个特定的目录下,每天checkout开发版的R,并且完成安装,应该怎么做呢?

cd path/to/put/R
svn checkout http://svn.r-project.org/R/trunk/ R-devel
cd R-devel
./configure --without-recommended-packages
make
make install
svn cleanup .

可以把上面的命令写成一个shell文件,然后通过daemons来自动调用。
比如在mac下,shell文件(autoUpdateRdev.sh)可以写成

#!/bin/bash
cd /qiuworld/R
svn checkout http://svn.r-project.org/R/trunk/ R-devel
cd R-devel
./configure --without-recommended-packages
make
make install
svn cleanup .

[……]

Read more

理解edgeR

from: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data

edgeR假设RNA-seq的reads计数对于每个基因来说是符合负二项式分布(Negative Binomial distribution)的。其公式可以表示为Pr(K[......]

Read more

使用DEXSeq分析NGS数据中的exon表达差异

对于RNA-seq,除了gene水平的差异分析外,还可以进行exon水平的差异分析。这时可以使用Bioconductor的DEXSeq软件包。

使用DEXSeq软件包,其输入为Count table。要生成这个Count table,与DESeq, edgeR类似的,都需要使用到htseq。DEXSeq提供了两个python角本,来调用htseq生成计数文件。

关于HTSeq的安装,请[……]

Read more

如何统计BAM文件中的reads数

当完成测序的比对工作之后,我们得到了bam/sam文件。那么,如何得到reads的统计数据呢?
这有很多途径:
1.读取日志文件。对于bowtie的日志,其中会包括如下的描述:
31991083 reads; of these:
31991083 (100.00%) were unpaired; of these:
6844445 (21.39%) aligned 0 times
18[……]

Read more

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

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

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

Read more

升级Bioconductor至2.13

Bioconductor已经升级至2.13了,我们如何升级自己的bioconductor呢?

首先最好是先升级自己R至3.02版本。至少也得是3.0版本。

然后在R中运行:

source("http://bioconductor.org/biocLite.R")
biocLite("BiocUpgrade")

在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