安装HTSeq，需要Python版本在2.5以上(但是在Python 3下不行)，并且需要安装NumPy。如果已经安装了NumPy的话，安装HTSeq并不困难。但是如果没有安装的话，可能会比较麻烦。
首先，需要安装python 2.5以上的版本。因为centOS 5所带的Python版本是2.4，无法满足HTSeq的安装要求。但是，我并不建议直接升组安装python2.5以上的版本，因为yum等很多功能都由python来实现，所以我的办法是全新安装一个python的版本到一个指定的目录下面去。
下载并安装Python 2.7.2.
[ouj@qiuworld.com ~]$ tar -xzvf Python-2.7.2.tgz [ouj@qiuworld.com ~]$ cd Python-2.7.2 [ouj@qiuworld.com Python-2.7.2]$ sudo yum install tcl #需要安装tcl/tk库 [ouj@qiuworld.com Python-2.7.2]$ sudo yum install tcl-devel [ouj@qiuworld.com Python-2.7.2]$ sudo yum install tk [ouj@qiuworld.com Python-2.7.2]$ sudo yum install tk-devel [ouj@qiuworld.com Python-2.7.2]$ ./configure --prefix=/opt/python2.7 --with-threads --enable-shared [ouj@qiuworld.com Python-2.7.2]$ make [ouj@qiuworld.com Python-2.7.2]$ sudo make install [ouj@qiuworld.com Python-2.7.2]$ sudo ln -s /opt/python2.7/bin/python /usr/bin/python2.7 [ouj@qiuworld.com Python-2.7.2]$ sudo echo '/opt/python2.7/lib'>> /etc/ld.so.conf.d/opt-python2.7.conf [ouj@qiuworld.com Python-2.7.2]$ sudo /sbin/ldconfig [ouj@qiuworld.com Python-2.7.2]$ python2.7 Python 2.7.2 (default, Nov 14 2011, 17:02:46) [GCC 4.1.2 20080704 (Red Hat 4.1.2-51)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import numpy Traceback (most recent call last): File "", line 1, in ImportError: No module named numpy
如果不新建/etc/ld.so.conf.d/opt-python2.7.conf文件并在当中写入一行/opt/python2.7/lib，将会得到如下错误：[......]
引子
闲来无事，翻看CELL杂志，发现很多基因组的图都使用circos来作图，于是就去circos.ca网上看了一眼，发现果然是个基因组研究绘图强大工具，应该为生物信息员掌握。于是花了点时间，来试验它的每一个设置。上网搜索发现，它的中文资料少得可怜，于是将心得体会做一总结，形成这一系列教程。希望能对急于掌握circos又不擅长阅读英语的人有所帮助。－－糗世界之糗糗
下载与安装
下载地址：http://circos.ca/software/download/circos/
circos是基于perl的脚本程序。它的安装难度在于安装好perl以及它所需要的模块。对于windows用户，可以试着安装Strawberry Perl或者ActiveState Perl。这两者都是不错的选择。对于Unix/linux/MacOS用户，很可能你已经安装了perl，否则，你可以到http://www.perl.org/get.html去下载安装。
我们需要测试一下perl的环境， UNIX/Linux/MacOS用户
> which perl /usr/bin/perl > perl -v This is perl, v5.10.0 built for ...
Windows用户
> perl -v This is perl, v5.10.0 built for ...
接着，我们将下载下来的circos程序解压，假设它的目录是circos-x.xx
> cd circos-x.xx > bin/circos -man
你可能得到一个帮助页面，那么你安装circos已经成功。也可能得到是出错信息。[......]
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, 生物信息学, 编程
近期评论