标签：软件

软件与网络应用

每天更新R原码 svn R source

~ admin ~ Leave a comment

在一个特定的目录下，每天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

~ admin ~ Leave a comment

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表达差异

~ admin ~ 12 Comments

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

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

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

Read more

如何统计BAM文件中的reads数

~ admin ~ 11 Comments

当完成测序的比对工作之后，我们得到了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）

~ admin ~ Leave a comment

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

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

Read more

升级Bioconductor至2.13

~ admin ~ 1 Comment

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

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

然后在R中运行：

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

在R中肩并肩（side-by-side）绘制箱线图（box-plot）

~ admin ~ Leave a comment

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

> 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当中的circos——OmicCircos

~ admin ~ 10 Comments

R也可以画出类似circos的图了，这是最新的OmicCircos给我们带来的新功能。

OmicCircos与circos相同的是，它可以画出连线(link), 2D图象(2D Tracks)[这其中包括scatter plots, line plots, histograms, heatmaps, polygons, boxplots]。但是遗憾的是，它不能绘制ideograms，glyp[……]

Read more