|
Glancing at the code for heatmap.2 I’m fairly sure that the default is to use dist, and it’s default is in turn to use euclidean distances.
The reason your attempt at passing distfun = dist(method = ‘euclidean’) didn’t work is thatdistfun (and hclustfun) are supposed to simply be name of functions. So if you want […]
## load the libraries required library(made4) require(graphics) library(cluster) library(stats)
#### read data in from a comma-delimited file with header array_data<-read.csv(“Overall_Outsidetreat_sort_top100.csv”,header=T, row.names=1, check.names=T) ## array_data is a data frame type ## rows are taxa_strings and columns are samples
# now transpose rows and columns # array_data<-t(array_data) ## now rows are chips and columns are taxa
dim(array_data) […]
我们在分析了差异表达数据之后,经常要生成一种直观图--热图(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])] > […]
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 […]
在之前的一节当中,图型名称有些混乱,从这一节开始将做如下统一(不全面):
英文名称 中文名称 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 […]
一,布局
R绘图所占的区域,被分成两大部分,一是外围边距,一是绘图区域。
外围边距可使用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″, […]
def stdDeviation(a): l=len(a) m=sum(a)/l d=0 for i in a:
d+=(i-m)**2 return (d*(1/l))**0.5
a=[5,6,8,9] print(stdDeviation(a)) ======== 1.5811388300841898
|
|
Recent Comments