生物信息学生R入门教程–R的数据结构

本教程版权为《糗世界》所有,任何组织或个人不得未经书面许可转载。

在阅读本节前请自测,如果全部都知道正确答案,则可跳过本节。

  • 什么是atomic vector的常用基本类型?
  • list和atomic vector有哪些不同?
  • matrix和data frame有什么不同?
  • data frame的列可以保存list吗?
  • data frame的每一行或者每一列的长度可以不一致吗?
  • 如何将factor正确地转换为数字?

[……]

Read more

生物信息学生R入门教程–输入与输出

本教程版权为《糗世界》所有,任何组织或个人不得未经书面许可转载。

在阅读本节前请自测,如果全部都知道正确答案,则可跳过本节。

  • 如何变更工作目录?
  • 文本输入时如何只读取前五行?
  • 读取xls文件需要调用什么软件包?
  • 读写数据库文件有哪些软件包?
  • 保存R数据时如何达到最大压缩比?
  • 使用Excel读取R输出的文本文件时有哪些潜在的错误可能?

[……]

Read more

生物信息学生R入门教程–绘图基础

本教程版权为《糗世界》所有,任何组织或个人不得未经书面许可转载。

在阅读本节前请自测,如果全部都知道正确答案,则可跳过本节。

  • 如何保存一个图像?
  • 如何绘制散点图?
  • 如何绘制饼图?
  • 如何绘制箱线图?
  • 如何绘制线图?
  • 如何绘制柱状图?

[……]

Read more

生物信息学生R入门教程–R编程

本教程版权为《糗世界》所有,任何组织或个人不得未经书面许可转载。

在阅读本节前请自测,如果全部都知道正确答案,则可跳过本节。

  • 常用函数是否全部掌握(见下表)?
  • 请编写一个函数并达到自己的设计要求。
  • 为什么R应该尽量避免for循环?
  • 什么是引用类(refference class)

[……]

Read more

如何得到exon numbering

以前,NCBI的genomic browser中是可以得到一个所谓的exon number的。其规则为同一个基因内所有exon按5′->3’排序,用排序序号来做为exon number。对于有重叠的exon使用同一序号,以排序先后为其加上a,b,c…的后缀。
art9-cc-figure-1(fix)
但是这其实是一个很不靠谱的表达方式,因为随着基因的转录本的变化,有可能会新增加exon。所以一些文章中使用exon number来标[……]

Read more

使用R绘制带比例的饼图 Pie Chart with Annotated Percentages

pie1 <- function (x, labels = names(x), edges = 200, radius = 0.8, clockwise = FALSE, 
    init.angle = if (clockwise) 90 else 0, density = NULL, angle = 45, 
    col = NULL, border = NULL, lty = NULL, main = NULL, percentage=T, rawNumber=F, digits=3, cutoff=0.01, 
    legend=F, legendpos="topright", legendcol=2, ...) 
{
    if (!is.numeric(x) || any(is.na(x) | x < 0)) 
        stop("'x' values must be positive.")
    if (is.null(labels)) 
        labels <- as.character(seq_along(x))
    else labels <- as.graphicsAnnot(labels)
    rawX <- x
    x <- c(0, cumsum(x)/sum(x))
    dx <- diff(x)
    nx <- length(dx)
    plot.new()
    pin <- par("pin")
    xlim <- ylim <- c(-1, 1)
    if (pin[1L] > pin[2L]) 
        xlim <- (pin[1L]/pin[2L]) * xlim
    else ylim <- (pin[2L]/pin[1L]) * ylim
    dev.hold()
    on.exit(dev.flush())
    plot.window(xlim, ylim, "", asp = 1)
    if (is.null(col)) 
        col <- if (is.null(density)) 
            c("white", "lightblue", "mistyrose", "lightcyan", 
                "lavender", "cornsilk", "pink")
        else par("fg")
    if (!is.null(col)) 
        col <- rep_len(col, nx)
    if (!is.null(border)) 
        border <- rep_len(border, nx)
    if (!is.null(lty)) 
        lty <- rep_len(lty, nx)
    angle <- rep(angle, nx)
    if (!is.null(density)) 
        density <- rep_len(density, nx)
    twopi <- if (clockwise) 
        -2 * pi
    else 2 * pi
    t2xy <- function(t) {
        t2p <- twopi * t + init.angle * pi/180
        list(x = radius * cos(t2p), y = radius * sin(t2p))
    }
    for (i in 1L:nx) {
        n <- max(2, floor(edges * dx[i]))
        P <- t2xy(seq.int(x[i], x[i + 1], length.out = n))
        polygon(c(P$x, 0), c(P$y, 0), density = density[i], angle = angle[i], 
            border = border[i], col = col[i], lty = lty[i])
        if(!legend){
        	P <- t2xy(mean(x[i + 0:1]))
	        lab <- as.character(labels[i])
	        if (!is.na(lab) && nzchar(lab)) {
	            lines(c(1, 1.05) * P$x, c(1, 1.05) * P$y)
	            text(1.1 * P$x, 1.1 * P$y, labels[i], xpd = TRUE, 
	                adj = ifelse(P$x < 0, 1, 0), ...)
	        }
        }
    }
    if (percentage) {
    	for (i in 1L:nx){
    		if(dx[i]>cutoff){
    			P <- t2xy(mean(x[i + 0:1]))
            	text(.8 * P$x, .8 * P$y, paste(formatC(dx[i]*100, digits=digits), "%", sep=""), xpd = TRUE, 
                	adj = .5, ...)
    		}
        }
    }else{
        if(rawNumber){
		for (i in 1L:nx){
    			if(dx[i]>cutoff){
    				P <- t2xy(mean(x[i + 0:1]))
            		text(.8 * P$x, .8 * P$y, rawX[i], xpd = TRUE, 
                		adj = .5, ...)
    			}
        	}
        }
    }
    if(legend) legend(legendpos, legend=labels, fill=col, border="black", bty="n", ncol = legendcol)
    title(main = main, ...)
    invisible(NULL)
}
 
pie1(1:5)

pie

使用Rcurl提交表单

有时候需要向在线工具批量提交申请。手动操作非常繁琐,使用curl工具提交会非常方便。

Rcurl已经整合好了相关代码,使用简便。下面就给出一个示例:

library(RCurl) #调用curl
curl <- getCurlHandle() #虚拟一个浏览器
curlSetOpt(cookiejar=tempfile(), curl=curl) #生成cookie
getURL("http://some.webpage.com/index.php", curl=curl) #填充cookie
response <- postForm("http://some.webpage.com/form.handle.page.php", style="HTTPPOST", data=data, file=fileUpload(filename="tobeuploadedfile.txt", contentType="text/plain"), curl=curl) #提交表单,表单对应的网址写在第一个参数中,具体帮助请使用?postForm查看。

使用R按照指定宽度计算覆盖率(calculate coverage by bins)

如何使用R来按照bin(比如说每200碱基宽)来计数覆盖度并生成一个bedgraph文件呢?
这个其实就是所谓的resample问题。可以使用?tileGenome来解决。

library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg19) ## human
tiles <- tileGenome(seqinfo(Hsapiens), tilewidth=200, cut.last.tile.in.chrom=TRUE) ##
tail(tiles[seqnames(tiles)=="chr1"]) ##
head(tiles[seqnames(tiles)=="chr2"]) ##
library(Rsamtools)
library(GenomicAlignments)
binnedAverage <- function(bins, numvar, mcolname)
{
    stopifnot(is(bins, "GRanges"))
    stopifnot(is(numvar, "RleList"))
    stopifnot(identical(seqlevels(bins), names(numvar)))
    bins_per_chrom <- split(ranges(bins), seqnames(bins))
    means_list <- lapply(names(numvar),
        function(seqname) {
            views <- Views(numvar[[seqname]],
                           bins_per_chrom[[seqname]])
            viewMeans(views)
        })
    new_mcol <- unsplit(means_list, as.factor(seqnames(bins)))
    mcols(bins)[[mcolname]] <- new_mcol
    bins
}
 
library(rtracklayer)
bams <- dir(pattern="*.bam$")
for(bam in bams){
	dat <- readGAlignmentsFromBam(bam)
	dat <- granges(dat)
	score <- coverage(dat)
	tiles.s <- keepSeqlevels(tiles, names(score))
	bins <- binnedAverage(tiles.s, score, "score")
	export(bins, gsub(".bam", ".bedgraph", bam), "BedGraph")
}

注意到viewMeans是和rowMeans类似的一种函数,相应的,它也有viewSums, viewMins, viewMaxs等。
如[……]

Read more

每天更新R原码 svn R source

在一个特定的目录下,每天checkout开发版的R,并且完成安装,应该怎么做呢?

cd path/to/put/R
svn checkout http://svn.r-project.org/R/trunk/ R-devel
cd R-devel
./configure --without-recommended-packages
make
make install
svn cleanup .

可以把上面的命令写成一个shell文件,然后通过daemons来自动调用。
比如在mac下,shell文件(autoUpdateRdev.sh)可以写成

#!/bin/bash
cd /qiuworld/R
svn checkout http://svn.r-project.org/R/trunk/ R-devel
cd R-devel
./configure --without-recommended-packages
make
make install
svn cleanup .

[……]

Read more