edgeR案例学习:deepSAGE分析

接到学习任务,要求搞清楚两个R包,edgeR和REDseq。这一篇就主要讲述edgeR中的案例之一,对应edgeR包中说明文件的案例12

要了解更多关于SAGE(Serial Analysis of Gene Expression),请阅读DeepSAGE—digital transcriptomics with high sensitivity, simple experimental protocol and multiplexing of samples

安装及调用edgeR

> source("http://bioconductor.org/biocLite.R")
BiocInstaller version 1.4.3, ?biocLite for help
> biocLite("edgeR")
BioC_mirror: http://bioconductor.org
Using R version 2.15, BiocInstaller version 1.4.3.
Installing package(s) 'edgeR'
trying URL 'http://www.bioconductor.org/packages/2.10/bioc/bin/macosx/leopard/contrib/2.15/edgeR_2.6.0.tgz'
Content type 'application/x-gzip' length 1558515 bytes (1.5 Mb)
opened URL
==================================================
downloaded 1.5 Mb
 
 
The downloaded binary packages are in
	/var/folders/Dj/Dj+bWjS7HxiNJ0kYFeKdTE+++TI/-Tmp-//RtmpcteSIZ/downloaded_packages
> library(edgeR)
Loading required package: limma

我们需要做准备一个名为targets.txt的文本文件,用于对样品的说明。接下来的工作就是读取数据。

> library(GEOquery)
Loading required package: Biobase
Loading required package: BiocGenerics
Welcome to Bioconductor
 
    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.
 
Setting options('download.file.method.GEOquery'='auto')
> gset<-getGEO("GSE10782",GSEMatrix = F)
File stored at: 
/var/folders/Dj/Dj+bWjS7HxiNJ0kYFeKdTE+++TI/-Tmp-//RtmpBoLGWD/GSE10782.soft.gz
Parsing....
Found 9 entities...
GPL9185 (1 of 9 entities)
GSM272105 (2 of 9 entities)
GSM272106 (3 of 9 entities)
GSM272318 (4 of 9 entities)
GSM272319 (5 of 9 entities)
GSM272320 (6 of 9 entities)
GSM272321 (7 of 9 entities)
GSM272322 (8 of 9 entities)
GSM272323 (9 of 9 entities)
There were 50 or more warnings (use warnings() to see the first 50)
> slotNames(gset)
[1] "header" "gsms"   "gpls"  
> x<-slot(gset,"gsms")
> setwd("~/Documents/shRNA/edgeR")
> for(i in 1:length(x)){
+	name<-Accession(x[[i]]);
+	data<-Table(dataTable(x[[i]]));
+	write.table(data,paste(name,".txt",sep=""),quote=F,sep="\t",row.names=F)
+ }
> targets<-read.delim("targets.txt",stringsAsFactors=F)
> targets
          files group                          description
1 GSM272105.txt  DCLK transgenic (Dclk1) mouse hippocampus
2 GSM272106.txt    WT          wild-type mouse hippocampus
3 GSM272318.txt  DCLK transgenic (Dclk1) mouse hippocampus
4 GSM272319.txt    WT          wild-type mouse hippocampus
5 GSM272320.txt  DCLK transgenic (Dclk1) mouse hippocampus
6 GSM272321.txt    WT          wild-type mouse hippocampus
7 GSM272322.txt  DCLK transgenic (Dclk1) mouse hippocampus
8 GSM272323.txt    WT          wild-type mouse hippocampus
> d<-readDGE(targets)
> colnames(d)
[1] "1" "2" "3" "4" "5" "6" "7" "8"
> colnames(d)<-c("DCLK1","WT1","DCLK2","WT2","DCLK3","WT3","DCLK4","WT4")

UTF8_E[……]

Read more

序列比对那点事儿

本来这应该是一本书,那样的话的确需要花一点心思,就写成一篇短文吧。

从字符比对开始说起吧。第一个问题最简单,如何判断两个字符串是相等的。

int strcmp(const char *s1, const char *s2)
{
  int ret = 0;
  while (!(ret = *(unsigned char *) s1 - *(unsigned char *) s2) && *s2) ++s1, ++s2;
 
  if (ret < 0)
    ret = -1;
  else if (ret > 0)
    ret = 1 ;
 
  return ret;
}

这种比较方式是逐一比较字符,如果遇到某个字符不一样了,就停下来,返回结果,否则就认为两个字符串是一样的。
[……]

Read more

R中的富集分析(一)--关于代谢途径(metabolic pathway)数据库

给定一些基因,从中找出联系,并能很好的显示给人们看,是一个生物信息学中非常常见的问题。其中pathway分析与gene ontology(GO)分析是两个很重要方面。基于种种原因,人们对GO分析往往比较注重,而随之而来的分析工具也多。而对PATH的分析,却因数据库的原因,手段寥寥。为此,我收集了一些关于pathway数据库方面的资料,总结如下:

1)代谢数据库有哪些?

这个问题很难全面[……]

Read more

使用R进行Gene Ontology(GO)富集分析过程中如何导出某一特定GO词条下所有的基因

我们在进行GO分析时,时常需要得到某一特定的GO term下所有涉及到的基因,但是R当中并没有明确的工具。下面的脚本片断可以解决这个问题。

> library("goProfiles")
> GOTermsList <- as.GOTerms.list(as.character(genelist.entrez.id),"Entrez","org.Hs.eg.db",onto="BP")
> ancestorsList<-getAncestorsLst(GOTermsList,onto="BP")
> yapply <- function(X,FUN, ...) {
+   index <- seq(length.out=length(X))
+   namesX <- names(X)
+   if(is.null(namesX))
+     namesX <- rep(NA,length(X))
+ 
+   FUN <- match.fun(FUN)
+   fnames <- names(formals(FUN))
+   if( ! "INDEX" %in% fnames ){
+     formals(FUN) <- append( formals(FUN), alist(INDEX=) )
+   }
+   if( ! "NAMES" %in% fnames ){
+     formals(FUN) <- append( formals(FUN), alist(NAMES=) ) 
+   } 
+   mapply(FUN,X,INDEX=index, NAMES=namesX,MoreArgs=list(...)) 
+ }
> ancestors<-do.call(rbind,yapply(ancestorsList,function(.ele){cbind(rep(NAMES,length(.ele)),.ele)}))
> colnames(ancestors)<-c("entrez_id","go_term")
> library("ChIPpeakAnno") #使用用其中addGeneIDs可以很方便地在种类ID之间进行转换
> anno<-addGeneIDs(as.character(genelist.entrez.id),"org.Hs.eg.db",c("symbol","genename"),"entrez_id")
> go.anno<-merge(anno,ancestors,by="entrez_id")
> c.cell.death<-lapply(cell.death.GO.terms,function(.ele,x){as.character(x[x$go_term==.ele,"symbol"])},go.anno)

今天又得到高人指点,有了新的更简单的脚本

> library(org.Hs.eg.db)
> help("org.Hs.egGO2ALLEGS")
> e<-mget(cell.death.GO.terms,org.Hs.egGO2ALLEGS)
> t<-lapply(e,function(.ele,x){intersect(as.character(x),.ele)},genelist.entrez.id)

R绘图基础(五)文氏图vennDiagram

文氏图是一种非常常用的图示手段,主要用于显示组与组之间重叠的程度。

R当中可以画文氏图的包有好几个,使用起来各有特点。最原始的工具,来自于2004年的《Venn Diagrams in R》–Duncan J.Murdoch, Journal of Statistical Software. 但是,现在已经无法找到venn包了。

之后使用得较为广泛的有两个工具,一个是LIMMA内嵌的vennDiagram, 另一个是gplots当中的venn。前者最多可以对3组数据画文氏图,后者可以最多对5组数据画文氏图。首先来介绍使用LIMMA的vennDiagram.

> source("http://bioconductor.org/biocLite.R")
> biocLite("limma")
> library(limma)
> hsb2<-read.table("http://www.ats.ucla.edu/stat/R/notes/hsb2.csv", sep=',', header=T)
> attach(hsb2)
> hw<-(write>=60)
> hm<-(math >=60)
> hr<-(read >=60)
> c3<-cbind(hw, hm, hr)
> a <- vennCounts(c3); a
     hw hm hr Counts
[1,]  0  0  0    113
[2,]  0  0  1     18
[3,]  0  1  0      8
[4,]  0  1  1      8
[5,]  1  0  0     12
[6,]  1  0  1      8
[7,]  1  1  0     11
[8,]  1  1  1     22
attr(,"class")
[1] "VennCounts"
> vennDiagram(a, include = "both", names = c("High Writing", "High Math", "High Reading"), cex = 1, counts.col = "red")
LIMMA绘制文氏图
LIMMA绘制文氏图

从上面的脚本我们可以知道,LIMMA在绘制文氏图的时候,先是对数据转换成bool类型的矩阵,而后进行分组计数,分组就包括所有可能的组合,得出数字后绘图。
[……]

Read more

Circos系列教程(三)突出标记Highlight

这一节的目标是画出下面的图

亮显强调
亮显强调

所谓突出标记,或者说亮显强调,多是通过大的反差明显或者符合色彩心理学的色块来将数据分组强调出来。在使用circos绘制基因组时,可以使用这一办法,将不同区域同一组内的基因亮显出来。[……]

Read more

R数据分析当中的化整为零(Split-Apply-Combine)策略

本文心得自:The Split-Apply-Combine Strategy for Data Analysis, Hadley Wickham, Journal of Statistical Software, April 2011, V.40.

引子:
我们常常会遇到这样的问题,数据量很大,并不需要依顺序来依次处理。合理分块处理,并最终整合起来是一个不错的选择。这也就是所谓的Split-Apply-Combine Strategy策略。这在速度上会有比做一个loop有优势,因为它可以并行处理数据。

什么时候我们需要使用到化整为零的策略呢?有以下三种情况:

  1. 数据需要分组处理
  2. 数据需要按照每行或者每列来处理
  3. 数据需要分级处理,和分组很类似,但是分级时需要考虑分级之间的关系。

化整为零策略有点类似于由Google推广的map-reduce策略。当然map-reduce策略的基础是网格,而这里的Split-Apply-Combine的基础完全可以是单机,甚至不支持并行处理的单机都可以。

然而,化整为零并不是一个很直观的编程过程。最直观的过程是使用Loop循环。这里使用一个例子来讲解一下如何实现化整为零策略。在plyr包中有数据ozone,它是一个三维矩阵(24X24X72),其中最后一维72是指的6年12个月每个月的结果。也就是ozone是一个包括了连续72个月24X24的三维矩阵数据。三维分别是lat,long,time。我们需要由对时间robust linear model之后的残基residuals。

> library(plyr) # need for dataset ozone
> library(MASS) # need for function rlm
> month <- ordered(rep(1:12, length=72)) #set time sequence
> #try one set
> one <- ozone[1,1,]
> model <- rlm(one ~ month - 1); model
Call:
rlm(formula = one ~ month - 1)
Converged in 9 iterations
 
Coefficients:
  month1   month2   month3   month4   month5   month6   month7   month8   month9  month10  month11  month12 
264.3964 259.2036 255.0000 252.0052 258.5089 265.3387 274.0000 276.6724 277.0000 285.0000 283.6036 273.1964 
 
Degrees of freedom: 72 total; 60 residual
Scale estimate: 4.45 
> deseas <- resid(model)

UTF8_E[……]

Read more

Circos系列教程(二)染色体示意图ideograms

本节的目标就是画出如下的图

circos绘制简单的ideogram

基础:circos作业流程

circos流程图

定义:
The symbolic representation of chromosomes are called ideograms.

circos为了能准确地画出染色体示意图,染色体的定义,位置,大小,以及显示的形式都是circos需要考虑的。这些要素需要在数据文件当中定义出来。[……]

Read more

Circos系列教程(一)安装

引子

Circular genome and data visualization with Circos (950 x 234)

闲来无事,翻看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已经成功。也可能得到是出错信息。[……]

Read more