如何使用bioconductor做liftover

当我们有不同版本的bed文件,比如说有的bed文件map到hg18, 有的bed文件map到hg19,或者对于小鼠来说有些文件map到mm9,有些文件map到mm10。我们如何实现对比对到不同基因组版本的数据进行比较呢?在比较之前,第一个步骤就是将不同的版本的比对文件统一至一个版本,这就需要使用到UCSC提供的liftOver工具。其可执行工具在:http://hgdownload.cse.ucs[……]

Read more

如何统计BAM文件中的reads数

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

使用Bioconductor下载GEO(Gene Expression Omnibus)上的数据

本文以GSE46106数据为例,讲述如何从GEO上下载数据。分为两个方面,第一,直接从GEO上下载表达数据,第二,直接从GEO上下载CEL文件。

首先我们来看第一种。

> library(GEOquery)
> gset <- getGEO("GSE46106", GSEMatrix =TRUE)
Found 1 file(s)
GSE46106_series_matrix.txt.gz
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/matrix/GSE46106_series_matrix.txt.gz'
ftp data connection made, file length 4110183 bytes
opened URL
==================================================
downloaded 3.9 Mb
 
File stored at: 
/var/folders/n4/11sc2xz13k56hl85z_h2rgq00000gn/T//RtmpBFQdpL/GPL570.soft
> length(gset)
[1] 1
> gset <- gset[[1]]
> head(pData(gset)[,1:5])
> # load NCBI platform annotation
> gpl <- annotation(gset)
> platf <- getGEO(gpl, AnnotGPL=TRUE)
File stored at: 
/var/folders/n4/11sc2xz13k56hl85z_h2rgq00000gn/T//RtmpBFQdpL/GPL570.annot.gz
There were 30 warnings (use warnings() to see them)
> ncbifd <- data.frame(attr(dataTable(platf), "table"))
> eset <- exprs(gset)
> head(eset[,1:5])
          GSM1123782 GSM1123783 GSM1123784 GSM1123785 GSM1123786
1007_s_at    10.1689    10.5247    10.8179    10.3539    10.4964
1053_at       9.6002     7.9436     9.8653     9.9733    10.1960
117_at        5.6808     5.0301     3.7654     2.7751     2.8179
121_at        4.2268     4.6148     4.6147     4.4977     4.6147
1255_g_at     2.1869     2.1869     2.1869     2.1869     2.1869
1294_at       2.1874     2.1874     2.1874     2.1874     2.1874
> head(ncbifd[,1:5])
         ID                                  Gene.title Gene.symbol Gene.ID UniGene.title
1 1007_s_at discoidin domain receptor tyrosine kinase 1        DDR1     780          <NA>
2   1053_at replication factor C (activator 1) 2, 40kDa        RFC2    5982          <NA>
3    117_at        heat shock 70kDa protein 6 (HSP70B')       HSPA6    3310          <NA>
4    121_at                                paired box 8        PAX8    7849          <NA>
5 1255_g_at     guanylate cyclase activator 1A (retina)      GUCA1A    2978          <NA>
6   1294_at ubiquitin-like modifier activating enzyme 7        UBA7    7318          <NA>

接着是第二种。

> library(GEOquery)
> getGEOSuppFiles("GSE46106")
[1] "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl/"
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl//GSE46106_RAW.tar'
ftp data connection made, file length 399247360 bytes
opened URL
===============================================
downloaded 380.8 Mb
 
trying URL 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE46nnn/GSE46106/suppl//filelist.txt'
ftp data connection made, file length 3308 bytes
opened URL
==================================================
downloaded 3308 bytes
> setwd("GSE46106/")
> dir()
[1] "filelist.txt"     "GSE46106_RAW.tar"
> untar("GSE46106_RAW.tar")
> files <- dir(pattern="gz$")
> sapply(files, gunzip)
GSM1123782_2147TG11_U133plus2.CEL.gz GSM1123783_2147TG15_U133plus2.CEL.gz  GSM1123784_2147TG1_U133plus2.CEL.gz 
                            33380229                             32529207                             33390958 
 GSM1123785_2147TG6_U133plus2.CEL.gz  GSM1123786_2277TG1_U133plus2.CEL.gz  GSM1123787_2277TG9_U133plus2.CEL.gz 
                            33210730                             33082826                             33531908 
GSM1123788_2665TG15_U133plus2.CEL.gz  GSM1123789_2665TG1_U133plus2.CEL.gz  GSM1123790_2665TG6_U133plus2.CEL.gz 
                            32232774                             33380782                             33432878 
 GSM1123791_3104TG1_U133plus2.CEL.gz  GSM1123792_3104TG3_U133plus2.CEL.gz  GSM1123793_3107TG1_U133plus2.CEL.gz 
                            32850644                             32588606                             32584342 
 GSM1123794_3107TG5_U133plus2.CEL.gz  GSM1123795_3143TG1_U133plus2.CEL.gz  GSM1123796_3143TG5_U133plus2.CEL.gz 
                            33358932                             32369391                             32806108 
 GSM1123797_3204TG1_U133plus2.CEL.gz  GSM1123798_3204TG5_U133plus2.CEL.gz  GSM1123799_3561TG1_U133plus2.CEL.gz 
                            32557782                             33016324                             32883029 
 GSM1123800_3561TG5_U133plus2.CEL.gz  GSM1123801_3611TG1_U133plus2.CEL.gz  GSM1123802_3611TG5_U133plus2.CEL.gz 
                            33454274                             32657648                             33126603 
 GSM1123803_3613TG1_U133plus2.CEL.gz  GSM1123804_3613TG5_U133plus2.CEL.gz  GSM1123805_3807TG1_U133plus2.CEL.gz 
                            32581902                             33193826                             32791917 
 GSM1123806_3807TG5_U133plus2.CEL.gz  GSM1123807_3824TG2_U133plus2.CEL.gz  GSM1123808_3887TG1_U133plus2.CEL.gz 
                            33266865                             33078422                             32479022 
 GSM1123809_3887TG5_U133plus2.CEL.gz  GSM1123810_3904TG1_U133plus2.CEL.gz  GSM1123811_3904TG5_U133plus2.CEL.gz 
                            33291287                             32645537                             33425796 
 GSM1123812_3936TG1_U133plus2.CEL.gz  GSM1123813_3936TG5_U133plus2.CEL.gz  GSM1123814_3963TG1_U133plus2.CEL.gz 
                            32939923                             33317196                             32725410 
 GSM1123815_3963TG5_U133plus2.CEL.gz  GSM1123816_4013TG1_U133plus2.CEL.gz  GSM1123817_4169TG2_U133plus2.CEL.gz 
                            33276082                             33151406                             33167461 
 GSM1123818_4175TG1_U133plus2.CEL.gz  GSM1123819_4272TG1_U133plus2.CEL.gz  GSM1123820_4400TG1_U133plus2.CEL.gz 
                            33139829                             33111635                             33206021 
 GSM1123821_4400TG5_U133plus2.CEL.gz  GSM1123822_4664TG1_U133plus2.CEL.gz  GSM1123823_4849TG1_U133plus2.CEL.gz 
                            33026806                             33364390                             33300210 
 GSM1123824_4888TG1_U133plus2.CEL.gz  GSM1123825_4913TG1_U133plus2.CEL.gz 
                            33299396                             33207528 
> filelist <- dir(pattern="CEL$")
> library(affy)
> library(annotate)
> data <- ReadAffy(filenames=filelist)
> affydb<-annPkgName(data@annotation,type="db")
> require(affydb, character.only=TRUE)
> eset<-rma(data,verbose=FALSE)
> eset.e <- exprs(eset)
> library(annaffy)
> symbols<-as.character(aafSymbol(as.character(rownames(eset)),affydb))
> genes<-as.character(aafUniGene(as.character(rownames(eset)),affydb))

视频链接:
youtube

土豆网[……]

Read more

Bioconductor并行(parallel)运算

当我们进入大数据时代之后,并行运算就十分重要了。什么是大数据呢?引用一句流行语吧:

Big data is hype. Everybody talks about big data; nobody knows exactly what it is.

对于生物学来说,大数据对应的应该是第二代测序技术的兴起。所以我们谈到并行运算,主要是为了处理测序数据。在R中,有许多并行运算的包,Hi[……]

Read more

使用pd.genomewidesnp.6为GenomeWideSNP_6做注释

代码如下:

> head(dat)
  Composite.Element.REF Call Confidence
1         SNP_A-2131660    2     0.0126
2         SNP_A-1967418    2     0.0029
3         SNP_A-1969580    2     0.0019
4         SNP_A-4263484    2     0.0043
5         SNP_A-1978185    0     0.0082
6         SNP_A-4264431    1     0.0088
> library(BiocInstaller)
Bioconductor version 2.13 (BiocInstaller 1.12.0), ?biocLite for help
> biocLite("pd.genomewidesnp.6")
> library(pd.genomewidesnp.6)
> con <- pd.genomewidesnp.6@getdb()
> dbListTables(con)
 [1] "featureSet"        "featureSetCNV"     "fragmentLength"    "fragmentLengthCNV" "pmfeature"         "pmfeatureCNV"     
 [7] "sequence"          "sequenceCNV"       "sqlite_stat1"      "table_info"       
> dbGetQuery(con, "select * from featureSet limit 3")[,1:10]
  fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
1      1 SNP_A-2131660          NA   rs2887286     1      1156131      0   p36.33        C        T
2      2 SNP_A-1967418          NA   rs1496555     1      2234251      0   p36.33        A        G
3      3 SNP_A-1969580          NA  rs41477744     1      2329564      0   p36.32        A        G
> dim(dat)
[1] 906600      3
> anno <- dbGetQuery(con, "select * from featureSet")
> dim(anno)
[1] 906600     13
> data <- cbind(dat, anno[match(dat[,1], anno[,2]),])
> dim(data)
[1] 906600     16
> miR146a <- data[grepl("miR146a", data[,"gene_assoc"], ignore.case=T), ]
> dim(miR146a)
[1] 49 16
> head(data[,1:13])
  Composite.Element.REF Call Confidence fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
1         SNP_A-2131660    2     0.0126      1 SNP_A-2131660          NA   rs2887286     1      1156131      0   p36.33        C        T
2         SNP_A-1967418    2     0.0029      2 SNP_A-1967418          NA   rs1496555     1      2234251      0   p36.33        A        G
3         SNP_A-1969580    2     0.0019      3 SNP_A-1969580          NA  rs41477744     1      2329564      0   p36.32        A        G
4         SNP_A-4263484    2     0.0043      4 SNP_A-4263484          NA   rs3890745     1      2553624      0   p36.32        C        T
5         SNP_A-1978185    0     0.0082      5 SNP_A-1978185          NA  rs10492936     1      2936870      1   p36.32        C        T
6         SNP_A-4264431    1     0.0088      6 SNP_A-4264431          NA  rs10489588     1      2951834      1   p36.32        C        T
> miR146a[miR146a$dbsnp_rs_id=="rs2910164",1:13]
       Composite.Element.REF Call Confidence fsetid    man_fsetid affy_snp_id dbsnp_rs_id chrom physical_pos strand cytoband allele_a allele_b
283011         SNP_A-2264525    0     0.0012 283011 SNP_A-2264525          NA   rs2910164     5    159912418      1      q34        C        G

如何在R中使用超几何分布计算p-value

问题,假设有两组ChIP-seq的数据A和B,我想知道它们之间是否相关。A组一共有50个数据,B组一共有60个数据,它们之间重叠的数据有30个,总计的binding位点如果是200个,那么如果它们相关的p-value如何计算呢?这个问题其实就是转变成抽到重叠数比30或者比30大的可能性是多大。所以我们需要计算抽到重叠数为29的可能性,然后用1减去它就可以了。
公式:p.value = phype[……]

Read more

如何让R自动保存历史记录(auto save R history before crash)

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

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

Read more

如何获取R中代码的运行效率

对于如何测试R代码的运行效率,首先可能想到的就是proc.time()以及system.time()了。

我们来看代码:

> ptm <- proc.time()
> for (i in 1:50) mad(stats::runif(500))
> proc.time() - ptm
   user  system elapsed 
  0.022   0.002   0.038 
> system.time(for(i in 1:100) mad(runif(1000)))
   user  system elapsed 
  0.025   0.002   0.028

这两个函数都只适合测试一整块代码的运行效率。

还有可以想到的工具应该就是rbenchmark了。

> random.array = function(rows, cols, dist=rnorm) 
+                   array(dist(rows*cols), c(rows, cols))
> random.replicate = function(rows, cols, dist=rnorm)
+                       replicate(cols, dist(rows))
> benchmark(replications=rep(100, 3),
+           random.array(100, 100),
+           random.array(100, 100),
+           columns=c('test', 'elapsed', 'replications'))
                    test elapsed replications
1 random.array(100, 100)   0.093          100
2 random.array(100, 100)   0.094          100
3 random.array(100, 100)   0.093          100
4 random.array(100, 100)   0.096          100
5 random.array(100, 100)   0.092          100
6 random.array(100, 100)   0.093          100

rb[……]

Read more