R的输入输出主要有几个类型：文本文件，xls文件，二进制文件，数据库文件，R文件，其它统计软件来源的文件，
首先是最简单的，文本文件，包括csv文件。假设我们有一个文件，题为qiuworld.com.txt，内容为
#this is a sample ArrayDataFile SourceName FactorValue GSM286765.CEL t_24h_rep1 treated_24h GSM286759.CEL t_12h_rep1 treated_12h GSM286763.CEL c_24h_rep2 control_24h GSM286760.CEL t_12h_rep2 treated_12h GSM286757.CEL c_12h_rep2 control_12h GSM286766.CEL t_24h_rep2 treated_24h GSM286756.CEL c_12h_rep1 control_12h GSM286762.CEL c_24h_rep1 control_24h
我们现在需要把它读入R，可以使用read.table命令。read.table还有几个快捷的形式，比如read.delim,read.delim2,read.csv,read.csv2。这几个快捷的方式帮助我们减少参数的书写。一般的，如果是csv文件，它的分隔符是逗号，字符串的两边会加上引号，可以直接使用read.csv(“文件名”)的方式读入数据。
> read.table("qiuworld.com.txt",header=TRUE,sep="\t",quote="",comment.char="#") ArrayDataFile SourceName FactorValue 1 GSM286765.CEL t_24h_rep1 treated_24h 2 GSM286759.CEL t_12h_rep1 treated_12h 3 GSM286763.CEL c_24h_rep2 control_24h 4 GSM286760.CEL t_12h_rep2 treated_12h 5 GSM286757.CEL c_12h_rep2 control_12h 6 GSM286766.CEL t_24h_rep2 treated_24h 7 GSM286756.CEL c_12h_rep1 control_12h 8 GSM286762.CEL c_24h_rep1 control_24h > read.delim("qiuworld.com.txt",header=TRUE,sep="\t",quote="",comment.char="#") ArrayDataFile SourceName FactorValue 1 GSM286765.CEL t_24h_rep1 treated_24h 2 GSM286759.CEL t_12h_rep1 treated_12h 3 GSM286763.CEL c_24h_rep2 control_24h 4 GSM286760.CEL t_12h_rep2 treated_12h 5 GSM286757.CEL c_12h_rep2 control_12h 6 GSM286766.CEL t_24h_rep2 treated_24h 7 GSM286756.CEL c_12h_rep1 control_12h 8 GSM286762.CEL c_24h_rep1 control_24h > read.csv("qiuworld.com.txt",header=TRUE,sep="\t",quote="",comment.char="#") ArrayDataFile SourceName FactorValue 1 GSM286765.CEL t_24h_rep1 treated_24h 2 GSM286759.CEL t_12h_rep1 treated_12h 3 GSM286763.CEL c_24h_rep2 control_24h 4 GSM286760.CEL t_12h_rep2 treated_12h 5 GSM286757.CEL c_12h_rep2 control_12h 6 GSM286766.CEL t_24h_rep2 treated_24h 7 GSM286756.CEL c_12h_rep1 control_12h 8 GSM286762.CEL c_24h_rep1 control_24h
有时候我们只需要读取文件的前几行，那么可以在上面的命令当中加入nrows参数，比如
> read.delim("qiuworld.com.txt",comment.char="#",nrows=5) ArrayDataFile SourceName FactorValue 1 GSM286765.CEL t_24h_rep1 treated_24h 2 GSM286759.CEL t_12h_rep1 treated_12h 3 GSM286763.CEL c_24h_rep2 control_24h 4 GSM286760.CEL t_12h_rep2 treated_12h 5 GSM286757.CEL c_12h_rep2 control_12h
[......]
Tags: bioconductor, 生物信息学, 编程
我们在分析了差异表达数据之后，经常要生成一种直观图－－热图(heatmap)。这一节就以基因芯片数据为例，示例生成高品质的热图。
比如
钢蓝渐白配色的热图
首先还是从最简单的heatmap开始。
> library(ggplot2) > library(ALL) #可以使用biocLite("ALL")安装该数据包 > data("ALL") > library(limma) > eset<-ALL[,ALL$mol.biol %in% c("BCR/ABL","ALL1/AF4")] > f<-factor(as.character(eset$mol.biol)) > design<-model.matrix(~f) > fit<-eBayes(lmFit(eset,design)) #对基因芯片数据进行分析，得到差异表达的数据 > selected <- p.adjust(fit$p.value[, 2]) <0.001 > esetSel <- eset[selected,] #选择其中一部分绘制热图 > dim(esetSel) #从这尺度上看，数目并不多，但也不少。如果基因数过多，可以分两次做图。 Features Samples 84 47 > library(hgu95av2.db) > data<-exprs(esetSel) > probes<-rownames(data) > symbol<-mget(probes,hgu95av2SYMBOL,ifnotfound=NA) > symbol<-do.call(rbind,symbol) > symbol[is.na(symbol[,1]),1]<-rownames(symbol)[is.na(symbol[,1])] > rownames(data)<-symbol[probes,1] #给每行以基因名替换探针名命名，在绘制热图时直接显示基因名。 > heatmap(data,cexRow=0.5)
使用heatmap函数默认颜色生成的热图
这个图有三个部分，样品分枝树图和基因分枝树图，以及热图本身。之所以对样品进行聚类分析排序，是因为这次的样品本身并没有分组。如果有分组的话，那么可以关闭对样品的聚类分析。对基因进行聚类分析排序，主要是为了色块好看，其实可以选择不排序，或者使用GO聚类分析排序。上面的这种热图，方便简单，效果非常不错。
接下来我们假设样品是分好组的，那么我们想用不同的颜色来把样品组标记出来，那么我们可以使用ColSideColors参数来实现。同时，我们希望变更热图的渐变填充色，可以使用col参数来实现。
[......]
Tags: bioconductor, 生物信息学, 编程
R当中的坐标中断一般都使用plotrix库中的axis.break(), gap.plot(), gap.barplot(), gap.boxplot()等几个函数来实现，例：
gap plot
> library(plotrix) > opar<-par(mfrow=c(3,2)) > plot(sample(5:7,20,replace=T),main="Axis break test",ylim=c(2,8)) > axis.break(axis=2,breakpos=2.5,style="gap") > axis.break(axis=2,breakpos=3.5,style="slash") > axis.break(axis=2,breakpos=4.5,style="zigzag") > twogrp<-c(rnorm(5)+4,rnorm(5)+20,rnorm(5)+5,rnorm(5)+22) > gap.plot(twogrp,gap=c(8,16,25,35), + xlab="X values",ylab="Y values",xlim=c(1,30),ylim=c(0,25), + main="Test two gap plot with the lot",xtics=seq(0,30,by=5), + ytics=c(4,6,18,20,22,38,40,42), + lty=c(rep(1,10),rep(2,10)), + pch=c(rep(2,10),rep(3,10)), + col=c(rep(2,10),rep(3,10)), + type="b") > gap.plot(21:30,rnorm(10)+40,gap=c(8,16,25,35),add=TRUE, + lty=rep(3,10),col=rep(4,10),type="l") > gap.barplot(twogrp,gap=c(8,16),xlab="Index",ytics=c(3,6,17,20), + ylab="Group values",main="Barplot with gap") > gap.barplot(twogrp,gap=c(8,16),xlab="Index",ytics=c(3,6,17,20), + ylab="Group values",horiz=TRUE,main="Horizontal barplot with gap") > twovec<-list(vec1=c(rnorm(30),-6),vec2=c(sample(1:10,40,TRUE),20)) > gap.boxplot(twovec,gap=list(top=c(12,18),bottom=c(-5,-3)), + main="Show outliers separately") > gap.boxplot(twovec,gap=list(top=c(12,18),bottom=c(-5,-3)),range=0, + main="Include outliers in whiskers") > par(opar)
从图像效果上来看，这样的坐标中断只能说实现了坐标中断，但效果上是非常一般的。甚至远不如excel, openoffice当中出图效果好。为此，我们需要对plotrix库中的gap.plot做出修改，以达到满意的效果。
[......]
Tags: bioconductor, 生物信息学, 编程
在之前的一节当中，图型名称有些混乱，从这一节开始将做如下统一（不全面）：
| 英文名称 | 中文名称 |
| bar | 条形图 |
| line | 线图 |
| area | 面积图 |
| pie | 饼图 |
| high-low | 高低图 |
| pareto | 帕累托图 |
| control | 控制图 |
| boxplot | 箱线图 |
| error bar | 误差条图 |
| scatter | 散点图 |
| P-P | P-P正态概率图 |
| Q-Q | Q-Q正态概率图 |
| sequence | 序列图 |
| ROC Curve | ROC分类效果曲线图 |
| Time Series | 时间序列图 |
好了，言归正传。那么什么又是点柱图(dot histogram)呢？之前我又称之为蜂群图(beeswarm)。还有称之为抖点图(jitter plots)。总之无论如何，在糗世界里我都称之为点柱图吧。
我们先看点柱图效果：
点柱图
以下是代码
> require(beeswarm) > data(breast) > head(breast) ER ESR1 ERBB2 time_survival event_survival 100.CEL.gz neg 8.372133 13.085894 39 1 103.CEL.gz pos 10.559356 9.491683 97 0 104.CEL.gz pos 12.299905 9.599574 11 1 105.CEL.gz pos 10.776632 9.681747 99 0 106.CEL.gz pos 10.505124 9.436763 40 1 107.CEL.gz neg 10.377741 8.695576 94 0 > require(plotrix) > cluster<-cluster.overplot(breast$event_survival, breast$time_survival) > png("dothist.png",width=1000,height=1000) > opar<-par(mfrow=c(3,3)) > plot (breast$event_survival, breast$time_survival, main="Multiple points on coordinate",col=as.numeric(breast$ER),xaxt="n",xlim=c(-1,2)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > plot(jitter(breast$event_survival), breast$time_survival, main="Using Jitter on x-axis",col=as.numeric(breast$ER),xaxt="n",xlim=c(-0.5,1.5)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > plot(jitter(breast$event_survival), jitter(breast$time_survival), main="Using Jitter on x and y-axis",col=as.numeric(breast$ER),xaxt="n",xlim=c(-0.5,1.5)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > sunflowerplot(breast$event_survival, breast$time_survival, main="Using Sunflowers",xaxt="n",xlim=c(-0.5,1.5)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > plot(cluster, main="Using cluster.overplot",col=as.numeric(breast$ER),xaxt="n",xlim=c(-0.5,1.5)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > count.overplot(jitter(breast$event_survival), jitter(breast$time_survival), main="Using cout.overplot",col=as.numeric(breast$ER),xaxt="n") > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > sizeplot(breast$event_survival, breast$time_survival, main="Using sizeplot",col=as.numeric(breast$ER),xaxt="n",xlim=c(-0.5,1.5)) > axis(1,at=c(0,1),labels=c("Censored","Metastasis")) > beeswarm(time_survival ~ event_survival, data = breast, + method = 'swarm', + pch = 16, pwcol = as.numeric(ER), + xlab = '', ylab = 'Follow-up time (months)', + labels = c('Censored', 'Metastasis')) > dev.off() quartz 2 > par(opar)
以下是解释[......]
Tags: bioconductor, 生物信息学, 编程
一，布局
Ｒ绘图所占的区域，被分成两大部分，一是外围边距，一是绘图区域。
外围边距可使用par()函数中的oma来进行设置。比如oma=c(4,3,2,1)，就是指外围边距分别为下边距：4行，左边距3行，上边距2行，右边距1行。很明显这个设置顺序是从x轴开始顺时针方向。这里的行是指可以显示1行普通字体。所以当我们使用mtext中的line参数时，设置的大小就应该是[0,行数)的开区间。当我们使用mtext在外围边距上书写内容时，设置mtext中的outer=TRUE即可。
绘图区域可使用par()函数中的mfrow, mfcol来进行布局。mfrow和mfcol可以使用绘图区域被区分为多个区域。默认值为mfrow(1,1)。
比如mfrow(2,3)就是指将绘图区域分成2行3列，并按行的顺序依次绘图填充；
比如mfcol(3,2)就是指将绘图区域分成3行2列，并按列的顺序依次绘图填充；
我们将每一个细分的绘图区域分为两个部分，一是绘图边距，一是主绘图。
绘图边距需要容纳的内容有坐标轴，坐标轴标签，标题。通常来讲，我们都只需要一个x轴，一个y轴，所以在设置时，一般的下边距和左边距都会大一些。如果多个x轴或者y轴，才考虑将上边距或者右边距放大一些。绘图边距可以使用par()函数中mar来设置。比如mar=c(4,3,2,1)，与外围边距的设置类似，是指绘图边距分别为下边距：4行，左边距3行，上边距2行，右边距1行。很明显这个设置顺序是从x轴开始顺时针方向。行的概念与之前的相同。也可以使用mai来设置。mai与mar唯一不同之处在于mai不是以行为单位，而是以inch为单位。
SOUTH< -1; WEST<-2; NORTH<-3; EAST<-4; GenericFigure <- function(ID, size1, size2) { plot(0:10, 0:10, type="n", xlab="X", ylab="Y") text(5,5, ID, col="red", cex=size1) box("plot", col="red") mtext(paste("cex",size2,sep=""), SOUTH, line=3, adj=1.0, cex=size2, col="blue") title(paste("title",ID,sep="")) } MultipleFigures <- function() { GenericFigure("1", 3, 0.5) box("figure", lty="dotted", col="blue") GenericFigure("2", 3, 1) box("figure", lty="dotted", col="blue") GenericFigure("3", 3, 1.5) box("figure", lty="dotted", col="blue") GenericFigure("4", 3, 2) box("figure", lty="dotted", col="blue") } par(mfrow=c(2,2),mar=c(6,4,2,1),oma=c(4,3,2,1)) MultipleFigures() box("inner", lty="dotted", col="green") box("outer", lty="solid", col="green") mtext("Outer Margin Area (oma) of South: 6", SOUTH, line=1, cex=1, outer=TRUE) plotline<-function(n,direc){ for(i in 0:n){ mtext(paste("line",i,sep=""), direc, line=i, cex=1, col="black", adj=1, outer=TRUE) } } plotline(4,SOUTH)
Ｒ绘图布局
UTF8_E[......]
Tags: bioconductor, 生物信息学, 编程
任务：有MoEx-1_0-st-v1芯片的注释文件，MoEx-1_0-st-v1.na31.mm9.probeset.csv。现在需要为其创建一个注释包，以方便注释及使用。
MoEx-1_0-st-v1.na31.mm9.probeset.csv可以从http://www.affymetrix.com/estore/browse/products.jsp?productId=131474&a[......]
Tags: bioconductor, 生物信息学, 编程
近期评论